chi-sqaures

import%20marimo%0A%0A__generated_with%20%3D%20%220.17.0%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20sys%0A%0A%20%20%20%20%23%20---%20WASM%20%2F%20Pyodide%20Environment%20---%0A%0A%20%20%20%20%23%201.%20Handle%20Plotting%20Libraries%0A%20%20%20%20%23%20Install%20matplotlib%20and%20seaborn%20%20%20%20%20%20%20%20%0A%20%20%20%20import%20matplotlib%0A%20%20%20%20import%20seaborn%20as%20sns%0A%20%20%20%20sns.set_theme(style%3D%22whitegrid%22)%20%23%20Apply%20a%20nice%20default%20style%0A%0A%20%20%20%20%23%202.%20Handle%20Scipy%0A%20%20%20%20%23%20This%20is%20the%20%22safe%22%20pattern%20for%20WASM.%0A%20%20%20%20import%20scipy%0A%20%20%20%20plt%20%3D%20matplotlib.pyplot%20%0A%20%20%20%20chi2_contingency%20%3D%20scipy.stats.chi2_contingency%0A%20%20%20%20fisher_exact%20%3D%20scipy.stats.fisher_exact%0A%20%20%20%20gaussian_kde%20%3D%20scipy.stats.gaussian_kde%0A%20%20%20%20chi2%20%3D%20scipy.stats.chi2%0A%0A%20%20%20%20%23%20---%20Your%20code%20below%20this%20cell%20can%20now%20safely%20use%20all%20functions%20---%0A%20%20%20%20return%20chi2%2C%20chi2_contingency%2C%20gaussian_kde%2C%20mo%2C%20np%2C%20pd%2C%20plt%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%23%20Bivariate%20Analysis%20of%20Categorical%20Variables%20Illustration%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22%23%23%20Set%20Up%20Illustration%20Baseline%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20%23%20---%20Dimensions%20---%0A%20%20%20%20n_rows%20%3D%20mo.ui.number(start%3D2%2C%20step%3D1%2C%20label%3D%22Number%20of%20Categories%20in%20First%20Variable%20(%23%20Rows)%22)%0A%20%20%20%20n_cols%20%3D%20mo.ui.number(%20start%3D2%2C%20step%3D1%2C%20label%3D%22Number%20of%20Categories%20in%20Second%20Variable%20(%23%20Columns)%22)%0A%0A%20%20%20%20%23%20---%20Sample%20size%20---%0A%20%20%20%20%23%20n%20%3D%20mo.ui.number(value%3D100%2C%20start%3D5%2C%20step%3D1%2C%20label%3D%22Sample%20size%20(n)%22)%0A%0A%20%20%20%20%23%20Random%20Seed%0A%20%20%20%20seed%20%3D%20mo.ui.number(start%3D42%2Cstep%3D1%2C%20label%3D%22RNG%20Seed%22)%0A%20%20%20%20mo.vstack(%5Bn_rows%2C%20n_cols%2C%20seed%5D)%0A%20%20%20%20return%20n_cols%2C%20n_rows%2C%20seed%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22Enter%20the%20contingency%20table%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo%2C%20n_cols%2C%20n_rows)%3A%0A%0A%20%20%20%20R%2C%20C%20%3D%20int(n_rows.value)%2C%20int(n_cols.value)%0A%0A%20%20%20%20%23%201)%20Row-label%20%22column%22%3A%20%5B%20corner%20%2C%20row1_name%20%2C%20...%20%2C%20rowR_name%20%5D%0A%20%20%20%20corner%20%3D%20mo.ui.text(%22%22%2C%20disabled%3DTrue)%0A%20%20%20%20row_label_cells%20%3D%20%5Bmo.ui.text(f%22Row%20%7Bi%2B1%7D%22)%20for%20i%20in%20range(R)%5D%0A%20%20%20%20row_col%20%3D%20mo.ui.array(%5Bcorner%2C%20*row_label_cells%5D%2C%20label%3D%22%20%22)%20%20%23%20blank%20label%20to%20keep%20UI%20minimal%0A%20%20%20%20row_col_arr%20%3D%20mo.vstack(row_col)%0A%0A%20%20%20%20%23%202)%20Data%20columns%3A%20each%20array%20is%20%5B%20col_header%20%2C%20x_11%20%2C%20x_21%20%2C%20...%20%2C%20x_R1%20%5D%20etc.%0A%20%20%20%20col_arrays%20%3D%20%5B%5D%0A%20%20%20%20col_vals%20%3D%20%5B%5D%0A%20%20%20%20for%20j%20in%20range(C)%3A%0A%20%20%20%20%20%20%20%20header%20%3D%20mo.ui.text(f%22Col%20%7Bj%2B1%7D%22)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20editable%20column%20header%0A%20%20%20%20%20%20%20%20col_cells%20%3D%20%5Bmo.ui.number(value%3D0%2C%20start%3D0%2C%20step%3D1)%20%20%20%20%23%20editable%20counts%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20_%20in%20range(R)%5D%0A%20%20%20%20%20%20%20%20col_cells_arr%20%3D%20mo.ui.array(%5Bheader%2C%20*col_cells%5D%2C%20label%3D%22%20%22)%0A%20%20%20%20%20%20%20%20col_vals.append(col_cells_arr)%0A%20%20%20%20%20%20%20%20col_arrays.append(mo.vstack(col_cells_arr))%0A%0A%20%20%20%20%23%203)%20Render%20as%20a%20table%20by%20stacking%20the%20column%20arrays%20horizontally%0A%20%20%20%20ui%20%3D%20mo.hstack(%5Brow_col_arr%2C%20*col_arrays%5D%2C%20gap%3D%22xs%22%2C%20align%3D%22center%22)%0A%20%20%20%20table_button%20%3D%20mo.ui.run_button(label%3D%22Generate%20Table%22)%0A%20%20%20%20mo.vstack(%5Bui%2C%20table_button%5D)%0A%20%20%20%20return%20C%2C%20R%2C%20col_vals%2C%20row_col%2C%20table_button%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20R%2C%20col_vals%2C%20mo%2C%20pd%2C%20row_col%2C%20table_button)%3A%0A%20%20%20%20if%20not%20table_button.value%3A%0A%20%20%20%20%20%20%20%20gate1%3DFalse%0A%20%20%20%20%20%20%20%20dis_ui%3Dmo.stop(True%2C%20mo.callout(mo.md(%22Click%20**Generate%20Table**%20to%20initialize%20data%20table%22)%2C%20kind%3D%22neutral%22))%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%23%20Labels%0A%20%20%20%20%20%20%20%20row_labels%20%3D%20list(row_col.value%5B1%3A%5D)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20skip%20the%20disabled%20corner%0A%20%20%20%20%20%20%20%20col_labels%20%3D%20%5Barr.value%5B0%5D%20for%20arr%20in%20col_vals%5D%20%20%20%23%20first%20entry%20in%20each%20column%20array%0A%0A%0A%20%20%20%20%20%20%20%20%23%20Data%20(gather%20row-wise%20from%20each%20column%20array)%0A%20%20%20%20%20%20%20%20data%20%3D%20%5B%5Bcol_vals%5Bj%5D.value%5B1%3A%5D%5Bi%5D%20for%20j%20in%20range(C)%5D%20for%20i%20in%20range(R)%5D%0A%0A%20%20%20%20%20%20%20%20obs_df%20%3D%20pd.DataFrame(data%2C%20index%3Drow_labels%2C%20columns%3Dcol_labels)%0A%20%20%20%20%20%20%20%20dis_ui%3Dmo.ui.table(obs_df)%0A%20%20%20%20%20%20%20%20gate1%3DTrue%0A%20%20%20%20dis_ui%0A%20%20%20%20return%20col_labels%2C%20gate1%2C%20obs_df%2C%20row_labels%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20%20R%2C%20C%20%3D%20int(n_rows.value)%2C%20int(n_cols.value)%0A%0A%20%20%20%20%23%20%23%20Editable%20category%20names%0A%20%20%20%20%23%20row_names%20%3D%20mo.ui.array(%5Bmo.ui.text(f%22Row%20%7Bi%2B1%7D%20name%22%2C%20label%3D%22%22)%20for%20i%20in%20range(R)%5D)%0A%20%20%20%20%23%20col_names%20%3D%20mo.ui.array(%5Bmo.ui.text(f%22Col%20%7Bj%2B1%7D%20name%22%2C%20label%3D%22%22)%20for%20j%20in%20range(C)%5D)%0A%0A%20%20%20%20%23%20%23%20Marginal%20probabilities%20(%25)%2C%20editable%0A%20%20%20%20%23%20row_probs%20%3D%20mo.ui.array(%5Bmo.ui.number(value%3D100.0%2FR%2C%20start%3D0.0%2C%20stop%3D100.0%2C%20step%3D0.01)%20for%20_%20in%20range(R)%5D)%0A%20%20%20%20%23%20col_probs%20%3D%20mo.ui.array(%5Bmo.ui.number(value%3D100.0%2FC%2C%20start%3D0.0%2C%20stop%3D100.0%2C%20step%3D0.01)%20for%20_%20in%20range(C)%5D)%0A%0A%20%20%20%20%23%20%23%20Expose%20these%20as%20outputs%20so%20downstream%20cells%20auto-run%20when%20any%20value%20changes%3A%0A%20%20%20%20%23%20ctrl_title%3Dmo.hstack(%5Bmo.md(%22**Row%20Names**%22)%2C%20mo.md(%22**Row%20Marginal%20Probabilities**%22)%2Cmo.md(%22**Column%20Names**%22)%2C%20mo.md(%22**Column%20Marginal%20Probabilities**%22)%20%5D)%0A%20%20%20%20%23%20ctrl_ui%3Dmo.hstack(%5Brow_names%2C%20row_probs%2Ccol_names%2C%20%20col_probs%5D)%0A%0A%20%20%20%20%23%20mo.vstack(%5Bctrl_title%2C%20ctrl_ui%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(obs_df)%3A%0A%20%20%20%20row_probs%3D(obs_df.sum(axis%3D1)%2Fobs_df.values.sum()).values*100%0A%20%20%20%20col_probs%3D(obs_df.sum(axis%3D0)%2Fobs_df.values.sum()).values*100%0A%20%20%20%20return%20col_probs%2C%20row_probs%0A%0A%0A%40app.cell%0Adef%20_(col_probs%2C%20mo%2C%20row_probs)%3A%0A%20%20%20%20row_sum%20%3D%20sum(row_probs)%0A%20%20%20%20col_sum%20%3D%20sum(col_probs)%0A%0A%20%20%20%20ok%20%3D%20abs(row_sum%20-%20100.0)%20%3C%201e-6%20and%20abs(col_sum%20-%20100.0)%20%3C%201e-6%0A%0A%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20f%22Row%20sum%20%3D%20%7Brow_sum%3A.2f%7D%25%2C%20Col%20sum%20%3D%20%7Bcol_sum%3A.2f%7D%25%20%E2%80%94%20%22%0A%20%20%20%20%20%20%20%20%2B%20(%22OK%20%E2%9C%85%22%20if%20ok%20else%20%22must%20sum%20to%20100%25%20%E2%9D%97%22)%2C%0A%20%20%20%20%20%20%20%20kind%3D%22success%22%20if%20ok%20else%20%22warn%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gate1%2C%20mo)%3A%0A%20%20%20%20gate1%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%23%20Random%20Null%20Tables%20%0A%20%20%20%20Let's%20draw%20some%20random%20tables%20under%20the%20Null%20Hypothesis.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(col_probs%2C%20mo%2C%20np%2C%20obs_df%2C%20row_probs)%3A%0A%20%20%20%20N%3Dint(obs_df.values.sum())%0A%20%20%20%20p_r1%20%3D%20np.array(row_probs%2C%20dtype%3Dfloat)%20%2F%20100.0%0A%20%20%20%20p_c1%20%3D%20np.array(col_probs%2C%20dtype%3Dfloat)%20%2F%20100.0%0A%0A%20%20%20%20%23%20normalize%20softly%20in%20case%20the%20sliders%2Finputs%20aren't%20exactly%20100%0A%20%20%20%20p_r%20%3D%20p_r1%20%2F%20p_r1.sum()%0A%20%20%20%20p_c%20%3D%20p_c1%20%2F%20p_c1.sum()%0A%0A%20%20%20%20%23%20joint%20cell%20probabilities%20under%20independence%20(flattened%20row-major)%0A%20%20%20%20p_orig%20%3D%20np.outer(p_r%2C%20p_c)%0A%20%20%20%20p%20%3D%20p_orig.reshape(-1)%0A%0A%20%20%20%20%23%20controls%0A%20%20%20%20M%20%3D%20mo.ui.number(value%3D1000%2C%20start%3D1%2C%20step%3D1%2C%20label%3D%22Number%20of%20tables%20(M)%22)%0A%20%20%20%20draw_btn%20%3D%20mo.ui.run_button(label%3D%22Draw%20tables%22)%0A%20%20%20%20mo.vstack(%5BM%2C%20draw_btn%5D%2C%20gap%3D%22md%22%2C%20align%3D%22start%22)%0A%20%20%20%20return%20M%2C%20N%2C%20draw_btn%2C%20p%2C%20p_orig%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20M%2C%20N%2C%20R%2C%20draw_btn%2C%20mo%2C%20np%2C%20p%2C%20seed)%3A%0A%20%20%20%20%23%20do%20the%20draws%0A%20%20%20%20if%20not%20draw_btn.value%3A%0A%20%20%20%20%20%20%20%20mo.stop(True%2C%20mo.callout(mo.md(%22Set%20M%20and%20click%20**Draw%20tables**.%22)%2C%20kind%3D%22neutral%22))%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20M_int%20%3D%20int(M.value)%0A%20%20%20%20%20%20%20%20rng%20%3D%20np.random.default_rng(None%20if%20int(seed.value)%20%3D%3D%200%20else%20int(seed.value))%0A%0A%20%20%20%20%20%20%20%20%23%20vectorized%20multinomial%3A%20shape%20(M%2C%20R*C)%20-%3E%20reshape%20to%20(M%2C%20R%2C%20C)%0A%20%20%20%20%20%20%20%20flat%20%3D%20rng.multinomial(N%2C%20p%2C%20size%3DM_int)%0A%20%20%20%20%20%20%20%20samples%20%3D%20flat.reshape(M_int%2C%20R%2C%20C)%0A%20%20%20%20return%20(samples%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20obs_df%2C%20samples)%3A%0A%20%20%20%20%23Expected%20cell%20values%20%0A%20%20%20%20%23%20Per-draw%20totals%0A%20%20%20%20Ns%20%20%20%20%20%20%3D%20samples.sum(axis%3D(1%2C%202))%5B%3A%2C%20None%2C%20None%5D%20%20%20%23%20(M%2C1%2C1)%0A%20%20%20%20rowtot%20%3D%20samples.sum(axis%3D2)%5B%3A%2C%20%3A%2C%20None%5D%20%20%20%20%20%20%20%20%20%20%20%23%20(M%2CR%2C1)%0A%20%20%20%20coltot%20%3D%20samples.sum(axis%3D1)%5B%3A%2C%20None%2C%20%3A%5D%20%20%20%20%20%20%20%20%20%20%20%23%20(M%2C1%2CC)%0A%0A%20%20%20%20%23%20Expected%20counts%20per%20draw%20under%20independence%20(MLE%20from%20that%20draw%E2%80%99s%20margins)%0A%20%20%20%20E_hat%20%3D%20(rowtot%20*%20coltot)%20%2F%20Ns%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20(M%2CR%2CC)%0A%0A%20%20%20%20%23%20Pearson%20chi-square%20per%20draw%0A%20%20%20%20contrib_numerator%3D(samples%20-%20E_hat)%0A%20%20%20%20contrib%20%3D%20%20contrib_numerator**2%2F%20np.where(E_hat%20%3E%200%2C%20E_hat%2C%20np.nan)%0A%20%20%20%20chi2_stats%20%3D%20np.nansum(contrib%2C%20axis%3D(1%2C%202))%20%20%20%20%20%20%20%23%20(M%2C)%0A%0A%20%20%20%20%23%20Observed%20Table%0A%20%20%20%20obs_row_tot%20%3D%20obs_df.values.sum(axis%3D1)%0A%20%20%20%20obs_col_tot%20%3D%20obs_df.values.sum(axis%3D0)%0A%20%20%20%20obs_expected%20%3D%20np.outer(obs_row_tot%2C%20obs_col_tot)%20%2F%20obs_df.values.sum()%0A%20%20%20%20observed_chi2%20%3D%20np.nansum((obs_df.values%20-%20obs_expected)**2%2F%20np.where(obs_expected%20%3E%200%2C%20obs_expected%2C%20np.nan))%0A%20%20%20%20return%20chi2_stats%2C%20obs_expected%2C%20observed_chi2%0A%0A%0A%40app.cell%0Adef%20_(col_labels%2C%20mo%2C%20obs_expected%2C%20pd%2C%20row_labels)%3A%0A%20%20%20%20obs_expected_df%20%3D%20pd.DataFrame(obs_expected)%0A%20%20%20%20obs_expected_df.index%3Drow_labels%20%0A%20%20%20%20obs_expected_df.columns%3Dcol_labels%20%0A%0A%20%20%20%20mo.vstack(%5Bmo.md(%22%23%23%20Expected%20Cell%20Values%20under%20Null%22)%2C%20mo.md(%22%24E_%7Bij%7D%3DN*p_%7Bij%7D%3DN*p_i*p_j%24%22)%2C%20mo.ui.table(obs_expected_df)%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(col_labels%2C%20mo%2C%20obs_expected%2C%20p_orig%2C%20pd%2C%20row_labels)%3A%0A%20%20%20%20obs_expected_var_df%20%3D%20pd.DataFrame(obs_expected*(1-p_orig))%0A%20%20%20%20obs_expected_var_df.index%3Drow_labels%20%0A%20%20%20%20obs_expected_var_df.columns%3Dcol_labels%20%0A%0A%20%20%20%20mo.vstack(%5Bmo.md(%22%23%23%20Expected%20Cell%20Values%20Variance%20under%20Null%22)%2C%20mo.md(%22%24Var_%7Bij%7D%3DN*p_%7Bij%7D(1-p_%7Bij%7D)%3DE_%7Bij%7D*(1-p_%7Bij%7D)%24%22)%2C%20mo.ui.table(obs_expected_var_df)%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gate1%2C%20mo)%3A%0A%20%20%20%20gate1%20%0A%20%20%20%20mo.md(%22%22%22%23%23%20Distribution%20of%20Each%20Cell%20Value%20-%20Expected%20Value%20from%20Random%20Draws%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(col_labels%2C%20np%2C%20plt%2C%20row_labels%2C%20samples)%3A%0A%20%20%20%20def%20fd_bin_width(x)%3A%0A%20%20%20%20%20%20%20%20x%20%3D%20np.asarray(x)%3B%20x%20%3D%20x%5Bnp.isfinite(x)%5D%0A%20%20%20%20%20%20%20%20if%20x.size%20%3C%202%3A%20return%201.0%0A%20%20%20%20%20%20%20%20q75%2C%20q25%20%3D%20np.percentile(x%2C%20%5B75%2C%2025%5D)%3B%20iqr%20%3D%20q75%20-%20q25%0A%20%20%20%20%20%20%20%20h%20%3D%202%20*%20iqr%20*%20(x.size%20**%20(-1%2F3))%0A%20%20%20%20%20%20%20%20return%20h%20if%20h%20%3E%200%20else%201.0%0A%0A%20%20%20%20def%20scott_bin_width(x)%3A%0A%20%20%20%20%20%20%20%20x%20%3D%20np.asarray(x)%3B%20x%20%3D%20x%5Bnp.isfinite(x)%5D%0A%20%20%20%20%20%20%20%20if%20x.size%20%3C%202%3A%20return%201.0%0A%20%20%20%20%20%20%20%20sigma%20%3D%20np.std(x%2C%20ddof%3D1)%0A%20%20%20%20%20%20%20%20h%20%3D%203.5%20*%20sigma%20*%20(x.size%20**%20(-1%2F3))%0A%20%20%20%20%20%20%20%20return%20h%20if%20h%20%3E%200%20else%201.0%0A%0A%20%20%20%20def%20compute_cell_figs(A%2C%20row_labels%2C%20col_labels%2C%20use%3D%22scott%22)%3A%0A%20%20%20%20%20%20%20%20A%20%3D%20np.asarray(A%2C%20float)%20%20%23%20shape%20(M%2CR%2CC)%0A%20%20%20%20%20%20%20%20M%2C%20R%2C%20C%20%3D%20A.shape%0A%20%20%20%20%20%20%20%20bw_rule%20%3D%20scott_bin_width%20if%20use%20%3D%3D%20%22scott%22%20else%20fd_bin_width%0A%0A%20%20%20%20%20%20%20%20means%20%3D%20np.empty((R%2C%20C))%3B%20vars_%20%3D%20np.empty((R%2C%20C))%0A%0A%20%20%20%20%20%20%20%20%23%20Create%20matplotlib%20subplot%20grid%0A%20%20%20%20%20%20%20%20fig%2C%20axes%20%3D%20plt.subplots(%0A%20%20%20%20%20%20%20%20%20%20%20%20R%2C%20C%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20figsize%3D(max(8%2C%203.5%20*%20C)%2C%20max(6%2C%203%20*%20R))%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20squeeze%3DFalse%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20sharex%3DFalse%2C%20%23%20Decouple%20axes%0A%20%20%20%20%20%20%20%20%20%20%20%20sharey%3DFalse%0A%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(R)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(C)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax%20%3D%20axes%5Bi%2C%20j%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20A%5B%3A%2C%20i%2C%20j%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20x%5Bnp.isfinite(x)%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20m%20%3D%20np.mean(x)%20if%20x.size%20else%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20v%20%3D%20np.var(x%2C%20ddof%3D1)%20if%20x.size%20%3E%201%20else%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20means%5Bi%2C%20j%5D%2C%20vars_%5Bi%2C%20j%5D%20%3D%20m%2C%20v%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20per-cell%20bins%3A%20data-driven%20width%2C%20full%20data%20range%20%2B%20padding%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20x.size%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bw%20%3D%20bw_rule(x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pad%20%3D%200.5%20*%20bw%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20xmin%2C%20xmax%20%3D%20np.min(x)%2C%20np.max(x)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Calculate%20bin%20edges%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20start%20%3D%20np.floor((xmin%20-%20pad)%20%2F%20bw)%20*%20bw%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20end%20%3D%20np.ceil((xmax%20%2B%20pad)%20%2F%20bw)%20*%20bw%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Ensure%20at%20least%20one%20bin%20if%20bw%20is%200%20or%20data%20is%20constant%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20end%20%3C%3D%20start%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20end%20%3D%20start%20%2B%201%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bins%20%3D%20np.arange(start%2C%20end%20%2B%20bw%2C%20bw)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20len(bins)%20%3C%202%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bins%20%3D%20np.linspace(xmin%20-%200.5%2C%20xmax%20%2B%200.5%2C%2010)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.hist(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bins%3Dbins%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20density%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.80%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20edgecolor%3D'black'%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20linewidth%3D0.5%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.set_xlim(-1%2C%201)%20%23%20Set%20dummy%20limits%20if%20no%20data%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20mean%20line%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20x.size%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.axvline(x%3Dm%2C%20linestyle%3D%22--%22%2C%20color%3D'k'%2C%20linewidth%3D1)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20Add%20annotation%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.text(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%200.98%2C%200.95%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22mean%3D%7Bm%3A.3g%7D%5Cnvar%3D%7Bv%3A.3g%7D%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ha%3D'right'%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20va%3D'top'%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20transform%3Dax.transAxes%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D9%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.set_title(f%22%7Brow_labels%5Bi%5D%7D%20%C3%97%20%7Bcol_labels%5Bj%5D%7D%22%2C%20fontsize%3D10)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.tick_params(axis%3D'both'%2C%20which%3D'major'%2C%20labelsize%3D8)%0A%0A%20%20%20%20%20%20%20%20fig.suptitle(%22Per-cell%20distributions%20(M%20draws)%20with%20mean%20%C2%B1%20var%22%2C%20fontsize%3D16)%0A%20%20%20%20%20%20%20%20fig.tight_layout(rect%3D%5B0%2C%200.03%2C%201%2C%200.95%5D)%20%23%20Adjust%20for%20suptitle%0A%0A%20%20%20%20%20%20%20%20return%20fig%2C%20means%2C%20vars_%0A%0A%0A%20%20%20%20%23%20---%20call%20it%20---%0A%20%20%20%20%23%20A%20%3D%20contrib_numerator%20%20%23%20(M%2CR%2CC)%0A%20%20%20%20fig_cells%2C%20means%2C%20sds%20%3D%20compute_cell_figs(samples%2C%20row_labels%2C%20col_labels)%0A%20%20%20%20fig_cells%20%23%20Marimo%20will%20render%20the%20matplotlib%20figure%0A%20%20%20%20return%20means%2C%20sds%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20R%2C%20col_labels%2C%20means%2C%20mo%2C%20pd%2C%20row_labels%2C%20sds)%3A%0A%20%20%20%20%23%20----%20tidy%20summary%20table%20----%0A%20%20%20%20summary%20%3D%20pd.DataFrame(%0A%20%20%20%20%20%20%20%20%5B(row_labels%5Bi%5D%2C%20col_labels%5Bj%5D%2C%20means%5Bi%2C%20j%5D%2C%20sds%5Bi%2C%20j%5D)%0A%20%20%20%20%20%20%20%20%20for%20i%20in%20range(R)%20for%20j%20in%20range(C)%5D%2C%0A%20%20%20%20%20%20%20%20columns%3D%5B%22row%22%2C%20%22col%22%2C%20%22mean%22%2C%20%22var%22%5D%0A%20%20%20%20)%0A%20%20%20%20summary%20%3D%20summary.round(6)%0A%20%20%20%20mo.ui.table(summary)%0A%20%20%20%20summary%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(C%2C%20R%2C%20gate1%2C%20mo)%3A%0A%20%20%20%20gate1%0A%20%20%20%20mo.md(f%22%22%22KDE%20from%20Random%20Draws%20vs%20%24%7B%7B%5C%5Cchi%7D%7D%5E2%20%20(dof%3D%7B(R-1)*(C-1)%7D)%24%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(chi2%2C%20chi2_stats%2C%20gaussian_kde%2C%20mo%2C%20n_cols%2C%20n_rows%2C%20np%2C%20plt)%3A%0A%0A%20%20%20%20stats%20%3D%20np.asarray(chi2_stats%2C%20dtype%3Dfloat)%0A%20%20%20%20stats%20%3D%20stats%5Bnp.isfinite(stats)%20%26%20(stats%20%3E%3D%200)%5D%0A%0A%20%20%20%20if%20stats.size%20%3C%202%3A%0A%20%20%20%20%20%20%20%20mo_fig%20%3D%20mo.callout(%22Need%20at%20least%202%20samples%20to%20estimate%20a%20KDE.%22%2C%20kind%3D%22warn%22)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20df%20%3D%20max(1%2C%20(int(n_rows.value)-1)*(int(n_cols.value)-1))%0A%0A%20%20%20%20%20%20%20%20%23%20x-grid%0A%20%20%20%20%20%20%20%20xmax%20%3D%20max(np.quantile(stats%2C%200.995)*1.2%2C%20float(chi2.ppf(0.999%2C%20df)))%0A%20%20%20%20%20%20%20%20x%20%3D%20np.linspace(0.0%2C%20xmax%2C%20400)%0A%0A%20%20%20%20%20%20%20%20%23%20----%20Decide%20whether%20to%20reflect%20based%20on%20where%20the%20plain%20KDE%20peaks%20----%0A%20%20%20%20%20%20%20%20threshold%20%3D%201.0%20%20%20%20%20%20%20%20%20%20%23%20mode%20%3C%20threshold%20%3D%3E%20use%20reflection%20(tune%20as%20you%20like)%0A%20%20%20%20%20%20%20%20bw%20%3D%20None%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20or%20e.g.%200.9%20for%20smoother%20%2F%200.7%20for%20sharper%0A%0A%20%20%20%20%20%20%20%20kde_plain%20%3D%20gaussian_kde(stats%2C%20bw_method%3Dbw)%0A%20%20%20%20%20%20%20%20y_plain%20%3D%20kde_plain(x)%0A%0A%20%20%20%20%20%20%20%20mode_x%20%3D%20x%5Bnp.argmax(y_plain)%5D%0A%20%20%20%20%20%20%20%20use_reflection%20%3D%20mode_x%20%3C%20threshold%0A%0A%20%20%20%20%20%20%20%20if%20use_reflection%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20data_reflect%20%3D%20np.r_%5Bstats%2C%20-stats%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20kde_ref%20%3D%20gaussian_kde(data_reflect%2C%20bw_method%3Dbw)%0A%20%20%20%20%20%20%20%20%20%20%20%20y_kde%20%3D%20kde_ref(x)%20%2B%20kde_ref(-x)%0A%20%20%20%20%20%20%20%20%20%20%20%20kde_name%20%3D%20f%22KDE%20(reflection%2C%20mode%E2%89%88%7Bmode_x%3A.2g%7D)%22%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20y_kde%20%3D%20y_plain%0A%20%20%20%20%20%20%20%20%20%20%20%20kde_name%20%3D%20f%22KDE%20(plain%2C%20mode%E2%89%88%7Bmode_x%3A.2g%7D)%22%0A%0A%20%20%20%20%20%20%20%20%23%20Reference%20chi-square%20PDF%0A%20%20%20%20%20%20%20%20y_chi%20%3D%20chi2.pdf(x%2C%20df)%0A%0A%20%20%20%20%20%20%20%20%23%20Create%20Matplotlib%20Figure%0A%20%20%20%20%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots()%0A%20%20%20%20%20%20%20%20ax.hist(%0A%20%20%20%20%20%20%20%20%20%20%20%20stats%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20bins%3D40%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20density%3DTrue%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.25%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22Samples%20(hist)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20edgecolor%3D'black'%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20linewidth%3D0.5%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20ax.plot(x%2C%20y_kde%2C%20label%3Dkde_name%2C%20linewidth%3D2)%0A%20%20%20%20%20%20%20%20ax.plot(x%2C%20y_chi%2C%20label%3Df%22%CF%87%C2%B2%20PDF%20(df%3D%7Bdf%7D)%22%2C%20linestyle%3D%22--%22%2C%20linewidth%3D2)%0A%0A%20%20%20%20%20%20%20%20ax.set_title(%22KDE%20vs%20%CF%87%C2%B2%20reference%22)%0A%20%20%20%20%20%20%20%20ax.set_xlabel(%22%CF%87%C2%B2%20statistic%22)%0A%20%20%20%20%20%20%20%20ax.set_ylabel(%22Density%22)%0A%20%20%20%20%20%20%20%20ax.legend()%0A%20%20%20%20%20%20%20%20ax.set_ylim(bottom%3D0)%0A%0A%20%20%20%20%20%20%20%20mo_fig%20%3D%20fig%20%23%20Marimo%20renders%20the%20figure%0A%0A%20%20%20%20mo_fig%0A%20%20%20%20return%20(stats%2C)%0A%0A%0A%40app.cell%0Adef%20_(chi2%2C%20mo%2C%20n_cols%2C%20n_rows%2C%20np%2C%20observed_chi2%2C%20plt%2C%20stats)%3A%0A%20%20%20%20def%20draw_ecdf()%3A%0A%20%20%20%20%20%20%20%20if%20stats.size%20%3C%202%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20fig%3Dmo.callout(%22Need%20at%20least%202%20samples%20to%20build%20an%20ECDF.%22%2C%20kind%3D%22warn%22)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20df%20%3D%20max(1%2C%20(int(n_rows.value)-1)*(int(n_cols.value)-1))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20x-range%20for%20smooth%20reference%20CDF%0A%20%20%20%20%20%20%20%20%20%20%20%20xmax%20%3D%20max(np.quantile(stats%2C%200.995)*1.2%2C%20float(chi2.ppf(0.999%2C%20df)))%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20np.linspace(0.0%2C%20xmax%2C%20600)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Reference%20CDF%0A%20%20%20%20%20%20%20%20%20%20%20%20ref_cdf%20%3D%20chi2.cdf(x%2C%20df)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Empirical%20CDF%0A%20%20%20%20%20%20%20%20%20%20%20%20s%20%3D%20np.sort(stats)%0A%20%20%20%20%20%20%20%20%20%20%20%20ecdf_y%20%3D%20np.arange(1%2C%20s.size%20%2B%201)%20%2F%20s.size%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Key%20verticals%0A%20%20%20%20%20%20%20%20%20%20%20%20x_ref95%20%3D%20float(chi2.ppf(0.95%2C%20df))%0A%20%20%20%20%20%20%20%20%20%20%20%20x_emp95%20%3D%20float(np.quantile(stats%2C%200.95))%0A%20%20%20%20%20%20%20%20%20%20%20%20x_obs%20%20%20%3D%20float(observed_chi2)%20if%20'observed_chi2'%20in%20locals()%20or%20'observed_chi2'%20in%20globals()%20else%20None%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Build%20Matplotlib%20figure%0A%20%20%20%20%20%20%20%20%20%20%20%20fig_cdf%2C%20ax%20%3D%20plt.subplots()%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.plot(x%2C%20ref_cdf%2C%20label%3Df%22%CF%87%C2%B2%20CDF%20(df%3D%7Bdf%7D)%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.plot(s%2C%20ecdf_y%2C%20label%3D%22Empirical%20CDF%20(samples)%22)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Vertical%20lines%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.axvline(x%3Dx_ref95%2C%20linestyle%3D%22--%22%2C%20color%3D'k'%2C%20alpha%3D0.7)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.text(x_ref95%2C%200.02%2C%20%22ref%2095%25%22%2C%20transform%3Dax.get_xaxis_transform()%2C%20ha%3D'right'%2C%20va%3D'bottom'%2C%20fontsize%3D9%2C%20rotation%3D90)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23ax.axvline(x%3Dx_emp95%2C%20linestyle%3D%22%3A%22%2C%20color%3D'k'%2C%20alpha%3D0.7)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23ax.text(x_emp95%2C%200.98%2C%20%22emp%2095%25%22%2C%20transform%3Dax.get_xaxis_transform()%2C%20ha%3D'left'%2C%20va%3D'top'%2C%20fontsize%3D9%2C%20rotation%3D90)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20x_obs%20is%20not%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.axvline(x%3Dx_obs%2C%20linestyle%3D%22-%22%2C%20color%3D'k'%2C%20alpha%3D0.9%2C%20linewidth%3D1.5)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ax.text(x_obs%2C%200.98%2C%20%22observed%20%CF%87%C2%B2%22%2C%20transform%3Dax.get_xaxis_transform()%2C%20ha%3D'right'%2C%20va%3D'top'%2C%20fontsize%3D9%2C%20rotation%3D90)%0A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Optionally%2C%20add%20horizontal%20guide%20at%200.95%20for%20clarity%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.axhline(y%3D0.95%2C%20linestyle%3D%22%3A%22%2C%20color%3D%22rgba(0%2C0%2C0%2C0.3)%22%2C%20linewidth%3D1)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.set_title(%22Empirical%20CDF%20vs%20reference%20%CF%87%C2%B2%20CDF%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.set_xlabel(%22%CF%87%C2%B2%20statistic%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.set_ylabel(%22CDF%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.legend(loc%3D%22lower%20right%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.set_ylim(0%2C%201.2)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20fig%20%3D%20fig_cdf%0A%0A%20%20%20%20%20%20%20%20return%20fig%0A%0A%20%20%20%20mo_fig2%3Ddraw_ecdf()%0A%20%20%20%20mo_fig2%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gate1%2C%20mo)%3A%0A%20%20%20%20gate1%0A%20%20%20%20mo.md(%22%22%22%23%20Fisher's%20Exact%20(Monte%20Carlo%20Version)%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(chi2_contingency%2C%20np)%3A%0A%0A%20%20%20%20def%20sample_fixed_margins(row_tot%2C%20col_tot%2C%20rng)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Sequential%20hypergeometric%20(Patefield-style)%20sampler.%0A%20%20%20%20%20%20%20%20Returns%20an%20(R%2C%20C)%20integer%20table%20with%20the%20exact%20given%20margins.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20r%20%3D%20row_tot.astype(int).copy()%0A%20%20%20%20%20%20%20%20c%20%3D%20col_tot.astype(int).copy()%0A%20%20%20%20%20%20%20%20R%2C%20C%20%3D%20len(r)%2C%20len(c)%0A%20%20%20%20%20%20%20%20T%20%3D%20np.zeros((R%2C%20C)%2C%20dtype%3Dint)%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(R%20-%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(C%20-%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20N_rem%20%3D%20c%5Bj%3A%5D.sum()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ngood%20%3D%20c%5Bj%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20draws%20%3D%20r%5Bi%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20rng.hypergeometric(ngood%2C%20N_rem%20-%20ngood%2C%20draws)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20T%5Bi%2C%20j%5D%20%3D%20x%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20r%5Bi%5D%20-%3D%20x%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%5Bj%5D%20-%3D%20x%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20last%20column%20in%20row%20i%0A%20%20%20%20%20%20%20%20%20%20%20%20T%5Bi%2C%20C%20-%201%5D%20%3D%20r%5Bi%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20c%5BC%20-%201%5D%20-%3D%20r%5Bi%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20r%5Bi%5D%20%3D%200%0A%20%20%20%20%20%20%20%20%23%20last%20row%20determined%20by%20remaining%20column%20totals%0A%20%20%20%20%20%20%20%20T%5BR%20-%201%2C%20%3A%5D%20%3D%20c%0A%20%20%20%20%20%20%20%20return%20T%0A%0A%20%20%20%20def%20stat_value(X%2C%20kind%3D%22pearson%22)%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Compute%20the%20test%20statistic%20for%20a%20table%20X%20(2D%20ndarray).%0A%20%20%20%20%20%20%20%20kind%3A%20%22pearson%22%20(%CF%87%C2%B2)%20or%20%22g2%22%20(likelihood-ratio%20%2F%20deviance).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20X%20%3D%20np.asarray(X%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20if%20kind%20%3D%3D%20%22pearson%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20chi2%2C%20_%2C%20_%2C%20_%20%3D%20chi2_contingency(X%2C%20correction%3DFalse)%20%20%23%20recompute%20E%20from%20X's%20margins%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20chi2%0A%20%20%20%20%20%20%20%20elif%20kind%20%3D%3D%20%22g2%22%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20G%5E2%20%3D%202%20*%20sum%20O_ij%20*%20log(O_ij%20%2F%20Ehat_ij)%2C%20with%200*log(0%2F.)%20defined%20as%200%0A%20%20%20%20%20%20%20%20%20%20%20%20n%20%3D%20X.sum()%0A%20%20%20%20%20%20%20%20%20%20%20%20r%20%3D%20X.sum(axis%3D1%2C%20keepdims%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20c%20%3D%20X.sum(axis%3D0%2C%20keepdims%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20E%20%3D%20r%20%40%20c%20%2F%20n%0A%20%20%20%20%20%20%20%20%20%20%20%20mask%20%3D%20(X%20%3E%200)%20%26%20(E%20%3E%200)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%202.0%20*%20np.sum(X%5Bmask%5D%20*%20np.log(X%5Bmask%5D%20%2F%20E%5Bmask%5D))%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError(%22Unknown%20kind%3B%20use%20'pearson'%20or%20'g2'%22)%0A%20%20%20%20return%20sample_fixed_margins%2C%20stat_value%0A%0A%0A%40app.cell%0Adef%20_(mo%2C%20np%2C%20obs_df)%3A%0A%20%20%20%20%23%20observed%20table%20(any%20R%20x%20C)%0A%20%20%20%20X_obs%20%3D%20np.asarray(obs_df.values%2C%20dtype%3Dint)%0A%20%20%20%20row_tot%20%3D%20X_obs.sum(axis%3D1)%0A%20%20%20%20col_tot%20%3D%20X_obs.sum(axis%3D0)%0A%0A%20%20%20%20%23%20controls%0A%20%20%20%20stat_radio%20%3D%20mo.ui.radio(%0A%20%20%20%20%20%20%20%20options%3D%5B%22pearson%22%2C%20%22g2%22%5D%2C%0A%20%20%20%20%20%20%20%20value%3D%22pearson%22%2C%0A%20%20%20%20%20%20%20%20label%3D%22Statistic%20for%20extremeness%22%0A%20%20%20%20)%0A%20%20%20%20M1%20%3D%20mo.ui.number(value%3D5000%2C%20start%3D100%2C%20step%3D100%2C%20label%3D%22Monte%20Carlo%20draws%20(M)%22)%0A%20%20%20%20run%20%3D%20mo.ui.run_button(label%3D%22Run%20MC%20Fisher%20(fixed%20margins)%22)%0A%0A%20%20%20%20mo.vstack(%5Bstat_radio%2C%20M1%2C%20run%5D%2C%20gap%3D%22md%22%2C%20align%3D%22start%22)%0A%20%20%20%20return%20X_obs%2C%20col_tot%2C%20row_tot%2C%20run%2C%20stat_radio%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20C%2C%0A%20%20%20%20M%2C%0A%20%20%20%20R%2C%0A%20%20%20%20X_obs%2C%0A%20%20%20%20chi2%2C%0A%20%20%20%20col_tot%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20row_tot%2C%0A%20%20%20%20run%2C%0A%20%20%20%20sample_fixed_margins%2C%0A%20%20%20%20seed%2C%0A%20%20%20%20stat_radio%2C%0A%20%20%20%20stat_value%2C%0A)%3A%0A%20%20%20%20def%20fisher_exact_show()%3A%0A%20%20%20%20%20%20%20%20rng%20%3D%20np.random.default_rng(None%20if%20int(seed.value)%20%3D%3D%200%20else%20int(seed.value))%0A%0A%20%20%20%20%20%20%20%20%23%20observed%20stat%0A%20%20%20%20%20%20%20%20kind%20%3D%20stat_radio.value%0A%20%20%20%20%20%20%20%20s_obs%20%3D%20stat_value(X_obs%2C%20kind%3Dkind)%0A%0A%20%20%20%20%20%20%20%20%23%20MC%20loop%0A%20%20%20%20%20%20%20%20M_int%20%3D%20int(M.value)%0A%20%20%20%20%20%20%20%20stats%20%3D%20np.empty(M_int%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20for%20m%20in%20range(M_int)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20T%20%3D%20sample_fixed_margins(row_tot%2C%20col_tot%2C%20rng)%0A%20%20%20%20%20%20%20%20%20%20%20%20stats%5Bm%5D%20%3D%20stat_value(T%2C%20kind%3Dkind)%0A%0A%20%20%20%20%20%20%20%20%23%20two-sided%20(by%20statistic)%3A%20%22as-or-more-extreme%22%0A%20%20%20%20%20%20%20%20%23%20add-one%20correction%20to%20avoid%20zero%20p%20when%20M%20small%0A%20%20%20%20%20%20%20%20p_mc%20%3D%20(1%20%2B%20np.sum(stats%20%3E%3D%20s_obs))%20%2F%20(M_int%20%2B%201)%0A%0A%20%20%20%20%20%20%20%20%23%20asymptotic%20chi-square%20ref%20(for%20context)%0A%20%20%20%20%20%20%20%20df%20%3D%20(R%20-%201)%20*%20(C%20-%201)%0A%20%20%20%20%20%20%20%20p_asym%20%3D%20float(1%20-%20chi2.cdf(s_obs%2C%20df))%0A%0A%20%20%20%20%20%20%20%20%23%20----%20display%20----%0A%20%20%20%20%20%20%20%20dis1%3Dmo.callout(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22**MC%20Fisher%20(fixed%20margins)**%20using%20*%7Bkind%7D*%20statistic%20%20%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22Observed%20stat%20%3D%20%7Bs_obs%3A.4g%7D%20%20%20%7C%20%20%20MC%20p%20%E2%89%88%20**%7Bp_mc%3A.4g%7D**%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22(asymptotic%20%CF%87%C2%B2%20df%3D%7Bdf%7D%3A%20p%20%E2%89%88%20%7Bp_asym%3A.4g%7D)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20kind%3D%22info%22%0A%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20%23%20Histogram%20%2B%20vline%20at%20observed%0A%20%20%20%20%20%20%20%20nbins%20%3D%20max(20%2C%20int(np.sqrt(M_int)))%0A%20%20%20%20%20%20%20%20fig_h%2C%20ax_h%20%3D%20plt.subplots()%0A%20%20%20%20%20%20%20%20ax_h.hist(%0A%20%20%20%20%20%20%20%20%20%20%20%20stats%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20bins%3Dnbins%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20density%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22MC%20null%22%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20alpha%3D0.6%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20edgecolor%3D'black'%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20linewidth%3D0.5%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20ax_h.axvline(x%3Ds_obs%2C%20color%3D'k'%2C%20linewidth%3D2)%0A%20%20%20%20%20%20%20%20ax_h.text(s_obs%2C%200.98%2C%20%22observed%22%2C%20transform%3Dax_h.get_xaxis_transform()%2C%20ha%3D'right'%2C%20va%3D'top'%2C%20fontsize%3D9%2C%20rotation%3D90)%0A%20%20%20%20%20%20%20%20ax_h.set_title(f%22MC%20Fisher%20null%20of%20%7Bkind.upper()%7D%20with%20observed%20stat%22)%0A%20%20%20%20%20%20%20%20ax_h.set_xlabel(f%22%7Bkind.upper()%7D%20statistic%22)%0A%20%20%20%20%20%20%20%20ax_h.set_ylabel(%22Density%22)%0A%20%20%20%20%20%20%20%20ax_h.legend()%0A%20%20%20%20%20%20%20%20dis2%20%3D%20fig_h%20%23%20Marimo%20renders%20the%20figure%0A%0A%20%20%20%20%20%20%20%20%23%20Empirical%20CDF%20%2B%20reference%20%CF%87%C2%B2%20CDF%20(context%20only)%0A%20%20%20%20%20%20%20%20xs%20%3D%20np.sort(stats)%3B%20ecdf%20%3D%20np.arange(1%2C%20M_int%2B1)%20%2F%20M_int%0A%20%20%20%20%20%20%20%20xgrid%20%3D%20np.linspace(0%2C%20max(xs.max()%2C%20s_obs)*1.1%2C%20600)%0A%0A%20%20%20%20%20%20%20%20fig_cdf%2C%20ax_cdf%20%3D%20plt.subplots()%0A%20%20%20%20%20%20%20%20ax_cdf.plot(xs%2C%20ecdf%2C%20label%3D%22Empirical%20CDF%20(MC)%22)%0A%20%20%20%20%20%20%20%20ax_cdf.plot(xgrid%2C%20chi2.cdf(xgrid%2C%20df)%2C%20linestyle%3D%22--%22%2C%20label%3Df%22%CF%87%C2%B2%20CDF%20(df%3D%7Bdf%7D)%22)%0A%0A%20%20%20%20%20%20%20%20ax_cdf.axvline(x%3Ds_obs%2C%20color%3D'k'%2C%20linewidth%3D2)%0A%20%20%20%20%20%20%20%20ax_cdf.text(s_obs%2C%200.98%2C%20%22observed%22%2C%20transform%3Dax_cdf.get_xaxis_transform()%2C%20ha%3D'right'%2C%20va%3D'top'%2C%20fontsize%3D9%2C%20rotation%3D90)%0A%0A%20%20%20%20%20%20%20%20ax_cdf.axhline(y%3D1%20-%20p_mc%2C%20linestyle%3D%22%3A%22%2C%20color%3D'k'%2C%20alpha%3D0.7)%0A%20%20%20%20%20%20%20%20ax_cdf.text(0.98%2C%201%20-%20p_mc%2C%20%221%20-%20MC%20p%22%2C%20transform%3Dax_cdf.transAxes%2C%20ha%3D'right'%2C%20va%3D'bottom'%2C%20fontsize%3D9)%0A%0A%20%20%20%20%20%20%20%20ax_cdf.set_title(%22Empirical%20CDF%20vs%20%CF%87%C2%B2%20CDF%20(reference%20shown%20for%20context)%22)%0A%20%20%20%20%20%20%20%20ax_cdf.set_xlabel(f%22%7Bkind.upper()%7D%20statistic%22)%0A%20%20%20%20%20%20%20%20ax_cdf.set_ylabel(%22CDF%22)%0A%20%20%20%20%20%20%20%20ax_cdf.set_ylim(0%2C%201)%0A%20%20%20%20%20%20%20%20ax_cdf.legend(loc%3D'lower%20right')%0A%0A%20%20%20%20%20%20%20%20dis3%20%3D%20fig_cdf%20%23%20Marimo%20renders%20the%20figure%0A%0A%20%20%20%20%20%20%20%20return%20mo.vstack(%5Bdis1%2C%20dis2%2C%20dis3%5D)%0A%0A%20%20%20%20if%20not%20run.value%3A%0A%20%20%20%20%20%20%20%20fisher_dis%3Dmo.callout(%22Set%20parameters%20and%20click%20**Run%20MC%20Fisher**.%22%2C%20kind%3D%22neutral%22)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20fisher_dis%3Dfisher_exact_show()%0A%20%20%20%20fisher_dis%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A