import%20marimo%0A%0A__generated_with%20%3D%20%220.23.9%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Bursty%20genes%20and%20the%20two-state%20promoter%20model%20in%20%60scribe%60%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20*This%20tutorial%20is%20a%20follow-up%20to*%20%5BModeling%20Assumptions%20for%20Single-Cell%20RNA-seq%20with%20%60scribe%60%5D(jurkat_cells.md).%20*It%20assumes%20you%20are%20comfortable%20with%20the%20ideas%20built%20up%20there%20%E2%80%94%20negative%20binomial%20counts%2C%20variable%20capture%2C%20parameterizations%2C%20hierarchical%20gene-specific%20%24p_g%24%20%E2%80%94%20and%20want%20to%20ask%20a%20question%20those%20models%20structurally%20cannot%20answer.*%0A%0A%20%20%20%20Every%20model%20in%20the%20previous%20tutorial%20was%2C%20at%20heart%2C%20a%20**negative%20binomial**.%20We%20made%20it%20more%20flexible%20in%20many%20ways%3A%20per-cell%20capture%20%24%5Cnu_c%24%2C%20friendlier%20coordinates%20%24(%5Cmu_g%2C%20%5Cphi)%24%2C%20joint%20low-rank%20guides%2C%20even%20a%20gene-specific%20success%20probability%20%24p_g%24.%20But%20all%20of%20those%20moves%20change%20the%20*parameters*%20of%20a%20negative%20binomial%20%E2%80%94%20none%20of%20them%20changes%20its%20*shape*.%20And%20a%20negative%20binomial%2C%20for%20any%20choice%20of%20%24(r%2C%20p)%24%2C%20is%20always%20either%20monotonically%20decreasing%20from%20zero%20or%20unimodal%20away%20from%20the%20boundary.%20**It%20can%20never%20be%20genuinely%20bimodal.**%0A%0A%20%20%20%20Some%20mammalian%20genes%2C%20however%2C%20have%20promoters%20that%20switch%20slowly%20between%20an%20OFF%20and%20an%20ON%20state.%20When%20that%20switching%20is%20slow%20relative%20to%20mRNA%20degradation%2C%20a%20population%20of%20cells%20splits%20into%20two%20groups%20%E2%80%94%20cells%20caught%20in%20a%20long%20OFF%20stretch%20(counts%20near%20zero)%20and%20cells%20caught%20in%20a%20long%20ON%20stretch%20(a%20second%20mode%20at%20moderate%20counts)%20%E2%80%94%20producing%20a%20two-mode%20count%20distribution.%20No%20negative%20binomial%20can%20fit%20that.%0A%0A%20%20%20%20This%20tutorial%20has%20**two%20phases**%3A%0A%0A%20%20%20%201.%20**Validate%20on%20synthetic%20data.**%20We%20generate%20counts%20from%20a%20*known*%20two-state%20process%20%E2%80%94%20some%20genes%20deeply%20bursty%2C%20some%20negative-binomial%20%E2%80%94%20fit%20the%20model%2C%20and%20check%20two%20things%3A%20(a)%20that%20it%20recovers%20the%20bimodality%20a%20negative%20binomial%20structurally%20cannot%2C%20and%20(b)%20*which*%20biophysical%20parameters%20are%20actually%20identifiable%20from%20snapshot%20counts%20and%20which%20are%20not.%20This%20is%20where%20the%20model's%20power%20is%20unambiguous%2C%20because%20we%20control%20the%20ground%20truth.%0A%20%20%20%202.%20**Apply%20to%20a%20real%20monoculture.**%20We%20then%20fit%20a%20%5Bdeeply-sequenced%20**K562**%20dataset%5D(https%3A%2F%2Fwww.10xgenomics.com%2Fdatasets%2F10k-human-k562-r-cells-singleplex-sample-1-standard).%20The%20honest%20spoiler%3A%20genuine%20promoter-bursting%20bimodality%20is%20*hard*%20to%20see%20in%20droplet%20total%20counts%20%E2%80%94%20for%20reasons%20we%20will%20make%20precise%20%E2%80%94%20so%20most%20genes%20sit%20in%20the%20negative-binomial%20limit%2C%20which%20the%20two-state%20model%20correctly%20reports.%20We%20use%20the%20fit%20to%20find%20the%20genes%20that%20deviate%20and%20ask%20whether%20they%20make%20biological%20sense.%0A%0A%20%20%20%20%3E%20**A%20note%20on%20dataset%20choice.**%20It%20is%20tempting%20to%20reach%20for%20a%20heterogeneous%20tissue%20(we%20tried%20%5BPBMC%5D(https%3A%2F%2Fwww.10xgenomics.com%2Fdatasets%2F10-k-peripheral-blood-mononuclear-cells-pbm-cs-from-a-healthy-donor-single-indexed-3-1-standard-4-0-0))%2C%20where%20marker%20genes%20look%20strikingly%20bimodal.%20But%20that%20conflates%20two%20very%20different%20mechanisms%3A%20*promoter%20bursting%20within%20one%20cell%20type*%20versus%20*a%20gene%20being%20on%20in%20one%20cell%20type%20and%20off%20in%20another*.%20The%20two-state%20model%20is%20built%20for%20the%20first%3B%20the%20second%20is%20a%20job%20for%20a%20mixture%20model.%20To%20study%20bursting%20cleanly%20we%20need%20a%20**monoculture**%20%E2%80%94%20Jurkat%20from%20the%20previous%20tutorial%20was%20one%2C%20but%20too%20shallow%3B%20K562%20is%20a%20deeper%20one.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20Import%20basic%20packages%0A%20%20%20%20from%20pathlib%20import%20Path%0A%20%20%20%20import%20pickle%0A%0A%20%20%20%20%23%20Import%20functionality%20to%20clear%20GPU%20memory%0A%20%20%20%20import%20gc%0A%20%20%20%20from%20jax%20import%20clear_caches%0A%0A%20%20%20%20%23%20Import%20our%20main%20package%0A%20%20%20%20import%20scribe%0A%0A%20%20%20%20%23%20Import%20useful%20tools%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20%23%20Import%20plotting%20packages%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20seaborn%20as%20sns%0A%0A%20%20%20%20%23%20Set%20our%20plotting%20style%20(totally%20optional)%0A%20%20%20%20scribe.viz.matplotlib_style()%0A%20%20%20%20return%20Path%2C%20clear_caches%2C%20gc%2C%20np%2C%20pickle%2C%20plt%2C%20scribe%2C%20sns%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20One%20backbone%2C%20different%20priors%20on%20the%20latent%20rate%0A%0A%20%20%20%20It%20is%20worth%20pausing%20to%20see%20how%20cleanly%20the%20three%20models%20in%20the%20paper%20relate.%20Every%20count%20model%20%60scribe%60%20implements%20shares%20the%20same%20backbone%3A%20an%20observed%20UMI%20count%20is%20a%20**Poisson%20draw%20from%20a%20latent%20per-gene%20rate**%2C%20and%20the%20per-cell%20capture%20efficiency%20%24%5Cnu_c%24%20thins%20that%20rate.%20What%20distinguishes%20the%20models%20is%20*the%20prior%20placed%20on%20the%20latent%20rate*%20%E2%80%94%20and%20that%20single%20choice%20controls%20two%20things%20at%20once%3A%20the%20**shape**%20the%20count%20distribution%20can%20take%2C%20and%20whether%20the%20clean%20Dirichlet%E2%80%93Multinomial%20composition%20(the%20backbone%20of%20%60scribe%60's%20differential-expression%20machinery)%20survives.%0A%0A%20%20%20%20**The%20negative%20binomial%20is%20the%20Poisson%E2%80%93Gamma%20compound.**%20This%20is%20the%20model%20from%20the%20previous%20tutorial%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Clambda_g%20%5Csim%20%5Cmathrm%7BGamma%7D(r_g%2C%20%5Ctheta)%2C%0A%20%20%20%20%5Cqquad%0A%20%20%20%20m_g%20%5Cmid%20%5Clambda_g%20%5Csim%20%5Cmathrm%7BPoisson%7D(%5Clambda_g)%0A%20%20%20%20%5C%3B%5C%3B%5CLongrightarrow%5C%3B%5C%3B%0A%20%20%20%20m_g%20%5Csim%20%5Cmathrm%7BNB%7D(r_g%2C%20p).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Biophysically%20this%20is%20the%20steady%20state%20of%20a%20**one-state%2C%20constitutively%20bursty%20promoter**%3A%20the%20gene%20is%20always%20able%20to%20fire%2C%20transcription%20arrives%20in%20geometric%20bursts%2C%20and%20the%20*unbounded*%20support%20of%20the%20Gamma%20reflects%20the%20fact%20that%20bursts%20can%2C%20in%20principle%2C%20pile%20up%20arbitrarily%20high.%0A%0A%20%20%20%20**The%20two-state%20promoter%20model%20replaces%20the%20Gamma%20with%20a%20Beta.**%20A%20promoter%20that%20switches%20%24%5Cmathrm%7BOFF%7D%5Cxrightarrow%7Bk%5E%2B%7D%5Cmathrm%7BON%7D%24%20and%20%24%5Cmathrm%7BON%7D%5Cxrightarrow%7Bk%5E-%7D%5Cmathrm%7BOFF%7D%24%2C%20producing%20mRNA%20at%20rate%20%24r%24%20only%20while%20ON%20and%20degrading%20it%20at%20rate%20%24%5Cgamma%24%20(we%20non-dimensionalize%20all%20rates%20by%20%24%5Cgamma%24%2C%20labeling%20them%20as%20%24%5Chat%20x%24)%2C%20has%20a%20steady-state%20count%20distribution%20that%20is%20*exactly*%20a%20**Poisson%E2%80%93Beta**%20compound%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20p%20%5Csim%20%5Cmathrm%7BBeta%7D(%5Chat%7Bk%7D%5E%2B_g%2C%20%5Chat%7Bk%7D%5E-_g)%2C%0A%20%20%20%20%5Cqquad%0A%20%20%20%20m_g%20%5Cmid%20p%20%5Csim%20%5Cmathrm%7BPoisson%7D(%5Chat%20r_g%5C%2C%20p).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Here%20%24p%20%5Cin%20%5B0%2C%201%5D%24%20is%20the%20**memory-weighted%20fraction%20of%20recent%20history%20the%20promoter%20spent%20ON**%20(weighted%20by%20the%20mRNA-lifetime%20kernel)%2C%20and%20%24%5Chat%20k%5E%2B%2C%20%5Chat%20k%5E-%24%20are%20the%20switching%20rates%20in%20units%20of%20the%20mRNA%20lifetime.%20The%20structural%20parallel%20to%20the%20negative%20binomial%20is%20exact%20%E2%80%94%20*same%20Poisson%2C%20different%20mixing%20distribution%20for%20the%20rate*%20%E2%80%94%20but%20the%20consequence%20is%20dramatic%3A%0A%0A%20%20%20%20-%20The%20Gamma%20has%20**unbounded%20support**%20%24%5B0%2C%20%5Cinfty)%24%3A%20it%20can%20be%20peaked%20but%20never%20has%20a%20second%20mode%20at%20the%20top.%0A%20%20%20%20-%20The%20Beta%20has%20**bounded%20support**%20%24%5B0%2C%201%5D%24%2C%20and%20when%20both%20shape%20parameters%20are%20below%201%20it%20is%20**U-shaped**%20%E2%80%94%20mass%20piled%20at%20*both*%20ends.%20Pushed%20through%20the%20Poisson%2C%20that%20produces%20a%20genuinely%20**bimodal**%20count%20distribution%3A%20a%20spike%20at%20zero%20(cells%20caught%20OFF)%20and%20a%20second%20mode%20near%20%24%5Cmathrm%7BPoisson%7D(%5Chat%20r)%24%20(cells%20caught%20ON).%20The%20bounded%20support%20is%20the%20biophysical%20signature%3A%20even%20a%20permanently-ON%20promoter%20caps%20out%20at%20%24%5Cmathrm%7BPoisson%7D(%5Chat%20r)%24%2C%20because%20production%20and%20degradation%20are%20finite.%0A%0A%20%20%20%20The%20next%20cell%20makes%20this%20concrete%20with%20synthetic%20draws%20%E2%80%94%20same%20Poisson%2C%20the%20only%20difference%20is%20Gamma%20vs.%20Beta%20for%20the%20rate.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20%23%20Same%20Poisson%20backbone%2C%20two%20different%20priors%20on%20the%20latent%20rate.%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(0)%0A%20%20%20%20_n%20%3D%2050_000%0A%20%20%20%20_target_mean%20%3D%2020.0%0A%0A%20%20%20%20%23%20---%20Poisson-Gamma%20(negative%20binomial)%3A%20unbounded%20Gamma%20rate%20---------%0A%20%20%20%20_r%20%3D%202.0%0A%20%20%20%20_lam%20%3D%20_rng.gamma(shape%3D_r%2C%20scale%3D_target_mean%20%2F%20_r%2C%20size%3D_n)%0A%20%20%20%20_counts_nb%20%3D%20_rng.poisson(_lam)%0A%0A%20%20%20%20%23%20---%20Poisson-Beta%20(two-state)%3A%20bounded%2C%20U-shaped%20Beta%20rate%20-----------%0A%20%20%20%20%23%20k%2B%20%3D%20k-%20%3D%200.3%20%3C%201%20makes%20the%20Beta%20U-shaped%20(promoter%20mostly%20fully%20ON%0A%20%20%20%20%23%20or%20fully%20OFF)%2C%20and%20r_hat%20sets%20where%20the%20%22ON%22%20mode%20lands.%0A%20%20%20%20_k_on%2C%20_k_off%20%3D%200.3%2C%200.3%0A%20%20%20%20_r_hat%20%3D%20_target_mean%20%2F%20(_k_on%20%2F%20(_k_on%20%2B%20_k_off))%20%20%23%20so%20E%5Bcount%5D%20%3D%20target%0A%20%20%20%20_p%20%3D%20_rng.beta(_k_on%2C%20_k_off%2C%20size%3D_n)%0A%20%20%20%20_counts_ts%20%3D%20_rng.poisson(_r_hat%20*%20_p)%0A%0A%20%20%20%20%23%20Plot%20distributions%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(9%2C%203.5)%2C%20sharey%3DTrue)%0A%20%20%20%20_bins%20%3D%20np.arange(0%2C%2090)%0A%0A%20%20%20%20_axes%5B0%5D.hist(_counts_nb%2C%20bins%3D_bins%2C%20color%3D%22C0%22%2C%20density%3DTrue)%0A%20%20%20%20_axes%5B0%5D.set_title(r%22Poisson%E2%80%93Gamma%20(negative%20binomial)%22)%0A%20%20%20%20_axes%5B0%5D.set_xlabel(%22count%20%24m%24%22)%0A%20%20%20%20_axes%5B0%5D.set_ylabel(%22density%22)%0A%0A%20%20%20%20_axes%5B1%5D.hist(_counts_ts%2C%20bins%3D_bins%2C%20color%3D%22C1%22%2C%20density%3DTrue)%0A%20%20%20%20_axes%5B1%5D.set_title(r%22Poisson%E2%80%93Beta%20(two-state%20promoter)%22)%0A%20%20%20%20_axes%5B1%5D.set_xlabel(%22count%20%24m%24%22)%0A%0A%20%20%20%20_fig.suptitle(%0A%20%20%20%20%20%20%20%20f%22Same%20mean%20(%24%5C%5Capprox%7B_target_mean%3A.0f%7D%24)%2C%20same%20Poisson%20%E2%80%94%20%22%0A%20%20%20%20%20%20%20%20f%22only%20the%20rate%20prior%20differs%22%0A%20%20%20%20)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Both%20panels%20have%20the%20*same%20mean*%20and%20the%20*same%20Poisson%20observation%20step*.%20The%20negative%20binomial%20(left)%20is%20unimodal%20%E2%80%94%20a%20peak%20near%20zero%20with%20a%20tail.%20The%20Poisson%E2%80%93Beta%20(right)%20is%20unmistakably%20**bimodal**%3A%20a%20population%20of%20cells%20with%20the%20promoter%20caught%20OFF%20(counts%20near%20zero)%20and%20a%20separate%20population%20caught%20ON%20(a%20second%20mode%20near%20%24%5Chat%20r%24).%20This%20is%20the%20qualitative%20behaviour%20the%20negative%20binomial%20structurally%20cannot%20produce%2C%20and%20it%20is%20the%20entire%20reason%20the%20two-state%20model%20exists.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20When%20does%20the%20two-state%20model%20collapse%20back%20to%20the%20negative%20binomial%3F%0A%0A%20%20%20%20Crucially%2C%20the%20two-state%20model%20**contains**%20the%20negative%20binomial%20as%20a%20special%20case.%20As%20the%20OFF%20rate%20grows%2C%20%24%5Chat%20k%5E-_g%20%5Cgg%201%24%20(the%20promoter%20spends%20only%20brief%2C%20frequent%20ON%20visits%2C%20and%20the%20long%20mRNA%20lifetime%20averages%20over%20many%20of%20them)%2C%20the%20Poisson%E2%80%93Beta%20limits%20*exactly*%20to%20a%20negative%20binomial%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cpi_%7B%5Ctext%7Bsteady-state%7D%7D(m)%20%5C%3B%5Cxrightarrow%7B%5C%3B%20%5Chat%20k%5E-_g%20%5Cgg%201%5C%3B%7D%5C%3B%0A%20%20%20%20%5Cmathrm%7BNegBinom%7D%5C!%5Cleft(m%20%5C%2C%5Cmiddle%7C%5C%2C%20%5Chat%20k%5E%2B_g%2C%5C%3B%20%5Cfrac%7B1%7D%7B1%20%2B%20%5Chat%20r_g%20%2F%20%5Chat%20k%5E-_g%7D%5Cright).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20This%20gives%20a%20precise%20biophysical%20dictionary%20for%20the%20negative%20binomial%20parameters%20we%20have%20been%20fitting%20all%20along%3A%0A%0A%20%20%20%20%7C%20Negative%20binomial%20parameter%20%7C%20Two-state%20biophysical%20meaning%20%7C%0A%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%7C%20shape%20%24r_%7B%5Cmathrm%7BNB%7D%7D%20%3D%20%5Chat%20k%5E%2B_g%24%20%7C%20dimensionless%20ON%20rate%20%24%3D%24%20**burst%20frequency**%20%7C%0A%20%20%20%20%7C%20burst%20size%20%24b_g%20%3D%20%5Chat%20r_g%20%2F%20%5Chat%20k%5E-_g%24%20%7C%20mean%20mRNA%20yield%20per%20ON%20visit%20%24%3D%24%20**burst%20size**%20%7C%0A%0A%20%20%20%20The%20burst%20size%20has%20a%20clean%20derivation%3A%20each%20ON%20visit%20lasts%2C%20on%20average%2C%20%241%2F%5Chat%20k%5E-_g%24%20mRNA%20lifetimes%20and%20produces%20mRNA%20at%20rate%20%24%5Chat%20r_g%24%2C%20so%20the%20expected%20yield%20per%20visit%20is%20%24%5Chat%20r_g%20%2F%20%5Chat%20k%5E-_g%24.%0A%0A%20%20%20%20The%20**contrapositive%20is%20the%20whole%20point%20of%20the%20model**%3A%20the%20genes%20for%20which%20the%20negative%20binomial%20fits%20*poorly*%20are%20exactly%20the%20genes%20whose%20data%20*reject*%20the%20fast-switching%20limit.%20A%20gene%20with%20%24%5Chat%20k%5E-_g%20%5Clesssim%201%24%20has%20ON%2FOFF%20dwell%20times%20comparable%20to%20the%20mRNA%20lifetime%3B%20the%20lifetime%20cannot%20smooth%20over%20the%20switching%2C%20and%20the%20steady-state%20distribution%20is%20genuinely%20bimodal.%20**That%20is%20what%20%22bursty%22%20means%20here%2C%20and%20%24%5Chat%20k%5E-_g%24%20is%20the%20knob%20that%20measures%20it.**%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20%60moment_delta%60%20parameterization%20(and%20the%20fit%20configuration)%0A%0A%20%20%20%20There%20is%20a%20practical%20wrinkle.%20The%20natural%20parameters%20%24(%5Chat%20r_g%2C%20%5Chat%20k%5E%2B_g%2C%20%5Chat%20k%5E-_g)%24%20are%20**not%20jointly%20identifiable**%3A%20you%20can%20scale%20%24%5Chat%20r_g%20%5Cto%20%5Clambda%20%5Chat%20r_g%24%20and%20%24%5Chat%20k%5E%2B_g%20%5Cto%20%5Clambda%20%5Chat%20k%5E%2B_g%24%20and%20leave%20both%20the%20mean%20and%20(to%20leading%20order)%20the%20Fano%20factor%20unchanged%20%E2%80%94%20a%20flat%2C%20banana-shaped%20ridge%20in%20parameter%20space%20that%20mean-field%20variational%20inference%20handles%20badly%2C%20exactly%20like%20the%20canonical%20%24(r%2C%20p)%24%20ridge%20in%20the%20first%20tutorial.%0A%0A%20%20%20%20%60scribe%60%20fixes%20this%20the%20same%20way%20it%20did%20before%3A%20reparameterize%20into%20data-identifiable%20combinations.%20The%20two-state%20model%20offers%20four%20mean-preserving%20parameterizations%3B%20we%20use%20**%60moment_delta%60**%2C%20whose%20sampled%20coordinates%20are%0A%0A%20%20%20%20-%20%24%5Cmu_g%24%20%E2%80%94%20the%20gene%20mean%20expression%20(a%20level%20knob%2C%20exactly%20the%20quantity%20most%20researchers%20think%20about)%3B%0A%20%20%20%20-%20%24F_g%20%3D%20%5Cmathrm%7BVar%7D%2F%5Cmu_g%20-%201%24%20%E2%80%94%20the%20**excess%20Fano%20factor**%2C%20the%20overdispersion-above-Poisson%20knob%3B%0A%20%20%20%20-%20%24%5Cdelta_g%20%3D%20%5Cdfrac%7B1%7D%7B%5Ckappa_g%20%2B%201%7D%20%5Cin%20(0%2C%201)%24%2C%20where%20%24%5Ckappa_g%20%3D%20%5Chat%20k%5E%2B_g%20%2B%20%5Chat%20k%5E-_g%24%20%E2%80%94%20the%20**regime%20knob**.%0A%0A%20%20%20%20The%20first%20two%20moments%20are%20preserved%20*by%20construction*%2C%20so%20%24%5Cmu_g%24%20and%20%24F_g%24%20are%20pinned%20down%20directly%20by%20the%20data.%20The%20regime%20coordinate%20%24%5Cdelta_g%24%20carries%20the%20part%20the%20data%20constrain%20most%20weakly%2C%20and%20it%20has%20a%20clean%20interpretation%3A%0A%0A%20%20%20%20-%20%24%5Cdelta_g%20%5Cto%200%24%20(%24%5Ckappa_g%20%5Cto%20%5Cinfty%24%2C%20fast%20switching)%20is%20the%20**negative%20binomial%20limit**%3B%0A%20%20%20%20-%20%24%5Cdelta_g%20%5Cto%201%24%20is%20the%20**extreme%20bursty**%20regime%20(the%20Beta%20degenerates%20to%20mass%20at%20%24%5C%7B0%2C%201%5C%7D%24%20%E2%80%94%20the%20promoter%20is%20fully%20ON%20or%20fully%20OFF).%0A%0A%20%20%20%20Mapping%20it%20to%20a%20bounded%20interval%20%24(0%2C%201)%24%20keeps%20the%20mean-field%20guide%20from%20wasting%20probability%20mass%20over%20an%20arbitrarily%20long%20ridge.%20From%20the%20fitted%20model%20we%20can%20read%20off%20any%20of%20the%20biophysical%20quantities%20%E2%80%94%20%24%5Chat%20k%5E%2B_g%24%2C%20%24%5Chat%20k%5E-_g%24%2C%20%24%5Chat%20r_g%24%2C%20the%20burst%20size%20%E2%80%94%20as%20deterministic%20functions%20of%20%24(%5Cmu_g%2C%20F_g%2C%20%5Cdelta_g)%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20exact%20%60scribe.fit%60%20configuration%20we%20use%20for%20the%20two-state%20fit%20is%20annotated%20below%3A%0A%0A%20%20%20%20-%20%60model%3D%22twostatevcp%22%60%20%E2%80%94%20the%20two-state%20promoter%20likelihood%20with%20per-cell%20**v**ariable%20**c**apture%20**p**robability.%20Because%20the%20Poisson%E2%80%93Beta%20is%20*closed%20under%20binomial%20thinning*%20(a%20thinned%20Poisson%20is%20just%20a%20Poisson%20with%20a%20scaled%20rate)%2C%20capture%20enters%20as%20%24%5Chat%20r_g%20%5Cmapsto%20%5Cnu_c%5C%2C%5Chat%20r_g%24%2C%20exactly%20the%20clean%20story%20we%20had%20for%20the%20negative%20binomial.%20(In%20**Phase%201**%20the%20synthetic%20data%20are%20generated%20with%20no%20capture%2C%20so%20we%20fit%20the%20plain%20%60model%3D%22twostate%22%60%20%E2%80%94%20model%20and%20data-generating%20process%20match%20exactly%3B%20in%20**Phase%202**%20the%20real%20K562%20data%20have%20capture%2C%20so%20we%20use%20%60twostatevcp%60.)%0A%20%20%20%20-%20%60parameterization%3D%22moment_delta%22%60%20%E2%80%94%20the%20%24(%5Cmu_g%2C%20F_g%2C%20%5Cdelta_g)%24%20coordinates%20above.%0A%20%20%20%20-%20%60unconstrained%3DTrue%60%20%E2%80%94%20fit%20in%20an%20unconstrained%20real%20space%20(positive%20params%20via%20a%20transform%2C%20the%20bounded%20%24%5Cdelta_g%24%20via%20a%20sigmoid)%2C%20which%20is%20numerically%20better%20behaved.%0A%20%20%20%20-%20%60priors%3D%7B%22inv_concentration%22%3A%20(0.0%2C%202.0)%7D%60%20%E2%80%94%20a%20**neutral%20prior%20on%20the%20regime%20coordinate**%20%24%5Cdelta%24.%20%60scribe%60's%20built-in%20default%20is%20deliberately%20NB-leaning%20(%60logit(%CE%B4)%20~%20Normal(-2%2C%202)%60)%3A%20it%20will%20not%20declare%20a%20gene%20bursty%20without%20strong%20evidence%2C%20which%20is%20the%20right%20conservative%20default%20for%20routine%20fits.%20But%20for%20a%20study%20that%20is%20*about*%20burstiness%2C%20that%20prior%20actively%20fights%20the%20signal%20%E2%80%94%20so%20we%20relax%20it%20to%20a%20symmetric%20%60Normal(0%2C%202)%60%20and%20let%20the%20data%2C%20not%20the%20prior%2C%20decide%20the%20regime.%20(We%20will%20see%20in%20Phase%201%20exactly%20how%20much%20this%20prior%20matters.)%0A%20%20%20%20-%20%60n_quad_nodes%3D256%60%20%E2%80%94%20the%20number%20of%20Gauss%E2%80%93Legendre%20nodes%20used%20to%20integrate%20the%20Poisson%E2%80%93Beta%20likelihood%20over%20the%20latent%20%24p%24.%20This%20is%20**the%20%60scribe%60%20default**%2C%20and%20it%20matters%20specifically%20in%20the%20bursty%20regime%3A%20a%20U-shaped%20Beta%20(%24%5Chat%20k%5E%2B_g%2C%20%5Chat%20k%5E-_g%20%3C%201%24)%20has%20near-singular%20mass%20at%20the%20endpoints%20%24p%20%5Cto%200%2C%201%24%2C%20and%20a%20coarse%20grid%20under-resolves%20it%20%E2%80%94%20systematically%20biasing%20the%20inferred%20regime%20coordinate%20toward%20the%20negative-binomial%20limit%2C%20which%20flattens%20the%20predicted%20ON%20mode.%0A%20%20%20%20-%20%60positive_transform%3D%22exp%22%60%20and%20%60optimizer_config%3D%7B...%20step_size%3A%201e-3%20...%7D%60%20%E2%80%94%20put%20every%20positive%20parameter%20on%20a%20log%20scale%20and%20use%20a%20large%20enough%20step.%20The%20regime%20coordinate%20converges%20*slowly*%3B%20an%20%60exp%60%20transform%20and%20a%20healthy%20step%20size%20are%20what%20let%20it%20actually%20reach%20the%20bursty%20solution%20rather%20than%20stalling%20NB-ward%20(see%20the%20negative%20binomial%20foil's%20%22one%20transform%20for%20everything%22%20note%20below).%0A%20%20%20%20-%20%60n_steps%60%2C%20%60batch_size%60%2C%20%60early_stopping%60%20%E2%80%94%20standard%20SVI%20controls%20with%20%60clipped_adam%60%20(gradient%20clipping%20keeps%20the%20bursty-regime%20gradients%20well%20behaved).%20We%20fit%20on%20the%20full%20gene%20set%20(no%20%60gene_coverage%60%20pre-filter)%2C%20so%20the%20two-state%20and%20negative%20binomial%20panels%20share%20exactly%20the%20same%20genes.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Phase%201%20%E2%80%94%20Validate%20on%20synthetic%20data%0A%0A%20%20%20%20Before%20trusting%20the%20model%20on%20real%20data%2C%20we%20ask%20the%20most%20basic%20question%3A%20*if%20the%20data%20really%20were%20generated%20by%20a%20two-state%20promoter%2C%20would%20the%20fit%20recover%20it%3F*%20We%20can%20answer%20that%20exactly%2C%20because%20we%20generate%20the%20data%20ourselves.%0A%0A%20%20%20%20We%20simulate%20%24G%20%3D%2020%24%20genes%20in%20%24C%20%3D%205%7B%2C%7D000%24%20cells.%20Half%20of%20the%20genes%20are%20**deeply%20bursty**%20(U-shaped%20Beta%2C%20%24%5Chat%20k%5E%2B_g%2C%20%5Chat%20k%5E-_g%20%3C%201%24%20%E2%80%94%20slow%20switching%2C%20genuinely%20bimodal)%3B%20half%20are%20**negative-binomial-like**%20controls%20(%24%5Chat%20k%5E-_g%20%5Cgg%201%24%20%E2%80%94%20fast%20switching).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20%23%20Ground-truth%20two-state%20parameters%3A%2010%20bursty%20(U-shaped)%20%2B%2010%20NB-like.%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(0)%0A%20%20%20%20n_cells_sim%20%3D%205_000%0A%20%20%20%20n_bursty%20%3D%2010%0A%20%20%20%20n_nb%20%3D%2010%0A%20%20%20%20_n_genes%20%3D%20n_bursty%20%2B%20n_nb%0A%0A%20%20%20%20k_on_true%20%3D%20np.concatenate(%5B%0A%20%20%20%20%20%20%20%20_rng.uniform(0.25%2C%200.45%2C%20n_bursty)%2C%20%20%20%23%20bursty%3A%20U-shaped%0A%20%20%20%20%20%20%20%20_rng.uniform(0.50%2C%202.00%2C%20n_nb)%2C%20%20%20%20%20%20%20%23%20NB-like%0A%20%20%20%20%5D)%0A%20%20%20%20k_off_true%20%3D%20np.concatenate(%5B%0A%20%20%20%20%20%20%20%20_rng.uniform(0.10%2C%200.25%2C%20n_bursty)%2C%20%20%20%23%20bursty%3A%20U-shaped%20(deep%20valley)%0A%20%20%20%20%20%20%20%20_rng.uniform(15.0%2C%2050.0%2C%20n_nb)%2C%20%20%20%20%20%20%20%23%20NB%20limit%20(fast%20OFF)%0A%20%20%20%20%5D)%0A%20%20%20%20r_hat_true%20%3D%20_rng.uniform(10.0%2C%20120.0%2C%20_n_genes)%20%20%20%23%20ON-mode%20location%0A%0A%20%20%20%20%23%20p_gc%20~%20Beta(k_on%2C%20k_off)%3B%20%20u_gc%20~%20Poisson(r_hat%20*%20p)%0A%20%20%20%20_p%20%3D%20_rng.beta(k_on_true%5BNone%2C%20%3A%5D%2C%20k_off_true%5BNone%2C%20%3A%5D%2C%20size%3D(n_cells_sim%2C%20_n_genes))%0A%20%20%20%20counts_sim%20%3D%20_rng.poisson(r_hat_true%5BNone%2C%20%3A%5D%20*%20_p).astype(np.int64)%0A%0A%20%20%20%20counts_sim.shape%0A%20%20%20%20return%20counts_sim%2C%20k_off_true%2C%20k_on_true%2C%20n_bursty%2C%20n_nb%2C%20r_hat_true%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20First%2C%20confirm%20the%20synthetic%20data%20actually%20look%20the%20way%20we%20intend%3A%20the%20bursty%20genes%20(top%20row)%20should%20show%20a%20deep-valley%2C%20two-mode%20shape%3B%20the%20NB-like%20controls%20(bottom%20row)%20should%20be%20unimodal.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(counts_sim%2C%20k_off_true%2C%20k_on_true%2C%20n_bursty%2C%20np%2C%20plt)%3A%0A%20%20%20%20_idx%20%3D%20%5B0%2C%201%2C%202%2C%203%2C%2010%2C%2011%2C%2012%2C%2013%5D%20%20%20%23%204%20bursty%2C%204%20NB-like%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(2%2C%204%2C%20figsize%3D(11%2C%205))%0A%20%20%20%20for%20_ax%2C%20_i%20in%20zip(_axes.ravel()%2C%20_idx)%3A%0A%20%20%20%20%20%20%20%20_ax.hist(%0A%20%20%20%20%20%20%20%20%20%20%20%20counts_sim%5B%3A%2C%20_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20bins%3Dnp.arange(0%2C%20counts_sim%5B%3A%2C%20_i%5D.max()%20%2B%202)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20density%3DTrue%2C%20histtype%3D%22step%22%2C%20color%3D%22black%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_grp%20%3D%20%22bursty%22%20if%20_i%20%3C%20n_bursty%20else%20%22NB-like%22%0A%20%20%20%20%20%20%20%20_ax.set_title(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22gene%20%7B_i%7D%20(%7B_grp%7D)%5Cn%24k%5E%2B%24%3D%7Bk_on_true%5B_i%5D%3A.2f%7D%2C%20%24k%5E-%24%3D%7Bk_off_true%5B_i%5D%3A.2f%7D%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D8%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22counts%22)%0A%20%20%20%20_fig.suptitle(%22Synthetic%20data%20%E2%80%94%20bursty%20genes%20(top)%20are%20bimodal%3B%20NB-like%20(bottom)%20are%20not%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Now%20fit%20both%20models%20to%20the%20synthetic%20counts.%20Because%20the%20synthetic%20data%20carry%20**no%20capture**%2C%20we%20fit%20the%20plain%20%60model%3D%22twostate%22%60%20%E2%80%94%20the%20fitted%20model%20is%20then%20*exactly*%20the%20data-generating%20process%2C%20which%20is%20what%20makes%20this%20a%20clean%20recovery%20test.%20We%20use%20the%20neutral%20regime%20prior%2C%20%60positive_transform%3D%22exp%22%60%2C%20and%20%60n_quad_nodes%3D256%60%20(the%20settings%20discussed%20above%3B%20the%20node%20count%20is%20what%20stops%20the%20U-shaped%20Beta%20from%20being%20mis-integrated%20and%20pulled%20NB-ward).%20The%20**negative%20binomial%20foil**%20is%20fit%20for%20the%20head-to-head.%20Both%20are%20cached%20to%20disk%20so%20the%20GPU%20fit%20runs%20once.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Path%2C%20clear_caches%2C%20counts_sim%2C%20gc%2C%20pickle%2C%20scribe)%3A%0A%20%20%20%20clear_caches()%0A%20%20%20%20gc.collect()%0A%0A%20%20%20%20%23%20Cache%20synthetic%20fits%20under%20the%20configured%20data%20root%20(out%20of%20the%20repo).%0A%20%20%20%20import%20json%20as%20_json%0A%20%20%20%20with%20Path(__file__).with_name(%22tutorial_paths.local.json%22).open(%0A%20%20%20%20%20%20%20%20%22r%22%2C%20encoding%3D%22utf-8%22%0A%20%20%20%20)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20_root%20%3D%20_json.load(_f)%5B%22SCRIBE_TUTORIAL_DATA_ROOT%22%5D%0A%20%20%20%20synth_dir%20%3D%20Path(_root).expanduser()%20%2F%20%22twostate_synthetic%22%0A%20%20%20%20synth_dir.mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20_out_path%20%3D%20synth_dir%20%2F%20%22synthetic_twostatevcp_moment_delta.pkl%22%0A%20%20%20%20if%20_out_path.exists()%3A%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22rb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20results_ts_sim%20%3D%20pickle.load(_f)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20results_ts_sim%20%3D%20scribe.fit(%0A%20%20%20%20%20%20%20%20%20%20%20%20counts_sim%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20model%3D%22twostate%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20parameterization%3D%22moment_delta%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20inference_method%3D%22svi%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20unconstrained%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_steps%3D200_000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20positive_transform%3D%22exp%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20priors%3D%7B%22inv_concentration%22%3A%20(0.0%2C%202.0)%7D%2C%20%20%20%23%20neutral%20regime%20prior%0A%20%20%20%20%20%20%20%20%20%20%20%20optimizer_config%3D%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22name%22%3A%20%22clipped_adam%22%2C%20%22step_size%22%3A%201e-3%2C%20%22grad_clip_norm%22%3A%2010.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_quad_nodes%3D256%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22wb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pickle.dump(results_ts_sim%2C%20_f)%0A%0A%20%20%20%20results_ts_sim%0A%20%20%20%20return%20results_ts_sim%2C%20synth_dir%0A%0A%0A%40app.cell%0Adef%20_(clear_caches%2C%20counts_sim%2C%20gc%2C%20pickle%2C%20scribe%2C%20synth_dir)%3A%0A%20%20%20%20clear_caches()%0A%20%20%20%20gc.collect()%0A%0A%20%20%20%20_out_path%20%3D%20synth_dir%20%2F%20%22synthetic_nbvcp_mean_odds.pkl%22%0A%20%20%20%20if%20_out_path.exists()%3A%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22rb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20results_nb_sim%20%3D%20pickle.load(_f)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20results_nb_sim%20%3D%20scribe.fit(%0A%20%20%20%20%20%20%20%20%20%20%20%20counts_sim%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20variable_capture%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20parameterization%3D%22mean_odds%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20prob_prior%3D%22gaussian%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20unconstrained%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_steps%3D100_000%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22wb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pickle.dump(results_nb_sim%2C%20_f)%0A%0A%20%20%20%20results_nb_sim%0A%20%20%20%20return%20(results_nb_sim%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20showcase%3A%20two-state%20recovers%20what%20the%20negative%20binomial%20cannot%0A%0A%20%20%20%20Posterior%20predictive%20checks%20for%20the%20bursty%20synthetic%20genes%20(integer%20indices%20%600%E2%80%937%60)%2C%20same%20genes%20in%20both%20panels.%20The%20negative%20binomial%20foil%20first%20%E2%80%94%20watch%20its%20band%20straddle%20the%20empty%20middle%20of%20the%20bimodal%20data.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(counts_sim%2C%20results_nb_sim%2C%20scribe)%3A%0A%20%20%20%20_pr%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_nb_sim%2C%20counts_sim%2C%20genes%3Dlist(range(8))%2C%0A%20%20%20%20%20%20%20%20n_rows%3D2%2C%20figsize%3D(10%2C%205)%2C%20n_samples%3D512%2C%0A%20%20%20%20)%0A%20%20%20%20_pr.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Now%20the%20two-state%20model%20on%20the%20same%20genes.%20Because%20its%20latent%20rate%20is%20Beta-distributed%2C%20it%20can%20place%20mass%20at%20zero%20*and*%20at%20a%20second%20mode%20at%20once%20%E2%80%94%20so%20its%20predictive%20band%20should%20follow%20the%20deep-valley%20shape%20the%20negative%20binomial%20could%20only%20average%20over.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(counts_sim%2C%20results_ts_sim%2C%20scribe)%3A%0A%20%20%20%20_pr%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_ts_sim%2C%20counts_sim%2C%20genes%3Dlist(range(8))%2C%0A%20%20%20%20%20%20%20%20n_rows%3D2%2C%20figsize%3D(10%2C%205)%2C%20n_samples%3D512%2C%0A%20%20%20%20)%0A%20%20%20%20_pr.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20That%20is%20the%20whole%20point%20of%20the%20model%2C%20demonstrated%20on%20data%20where%20we%20*know*%20the%20answer%3A%20the%20negative%20binomial%20cannot%20represent%20the%20bimodality%20at%20any%20setting%20of%20its%20parameters%3B%20the%20two-state%20model%20recovers%20it.%20Now%20we%20ask%20the%20harder%2C%20more%20honest%20question%20%E2%80%94%20given%20a%20good%20PPC%2C%20*which%20biophysical%20parameters%20did%20we%20actually%20learn%3F*%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(k_off_true%2C%20k_on_true%2C%20np%2C%20r_hat_true%2C%20results_ts_sim)%3A%0A%20%20%20%20%23%20Posterior%20medians%20for%20the%20fit%2C%20in%20BOTH%20coordinate%20systems.%20moment_delta%0A%20%20%20%20%23%20*samples*%20(mu%2C%20excess_fano%2C%20delta)%20and%20*derives*%20(k_on%2C%20k_off%2C%20r_hat)%20as%0A%20%20%20%20%23%20deterministic%20functions%20of%20them%20%E2%80%94%20so%20the%20same%20posterior%20can%20look%20pinned%20in%0A%20%20%20%20%23%20one%20chart%20and%20degenerate%20in%20the%20other.%20We%20extract%20both%20to%20show%20the%20contrast.%0A%20%20%20%20_post%20%3D%20results_ts_sim.get_posterior_samples(%0A%20%20%20%20%20%20%20%20n_samples%3D400%2C%20store_samples%3DFalse%2C%20convert_to_numpy%3DTrue%0A%20%20%20%20)%0A%20%20%20%20%23%20Derived%20biophysical%20rates%20(top%20row%20of%20the%20recovery%20plot).%0A%20%20%20%20kon_fit%20%3D%20np.median(_post%5B%22k_on%22%5D%2C%20axis%3D0)%0A%20%20%20%20koff_fit%20%3D%20np.median(_post%5B%22k_off%22%5D%2C%20axis%3D0)%0A%20%20%20%20rhat_fit%20%3D%20np.median(_post%5B%22r_hat%22%5D%2C%20axis%3D0)%0A%20%20%20%20%23%20Sampled%20moment_delta%20coordinates%20(bottom%20row)%20%E2%80%94%20what%20the%20data%20constrains%3A%0A%20%20%20%20%23%20mu%20%3D%20E%5Bcount%5D%2C%20excess_fano%20%3D%20Var%2FMean%20-%201%2C%20delta%20%3D%201%2F(kappa%2B1)%20in%20(0%2C%201).%0A%20%20%20%20mu_fit%20%3D%20np.median(_post%5B%22mu%22%5D%2C%20axis%3D0)%0A%20%20%20%20fano_fit%20%3D%20np.median(_post%5B%22excess_fano%22%5D%2C%20axis%3D0)%0A%20%20%20%20delta_fit%20%3D%20np.median(_post%5B%22inv_concentration%22%5D%2C%20axis%3D0)%0A%20%20%20%20%23%20log10%20inter-quantile%20width%20of%20the%20k_off%20posterior%20(%3E1%20spans%20%3E10x%20%3D%20non-identified)%0A%20%20%20%20_lo%2C%20_hi%20%3D%20np.percentile(_post%5B%22k_off%22%5D%2C%20%5B16%2C%2084%5D%2C%20axis%3D0)%0A%20%20%20%20koff_logwidth%20%3D%20np.log10(_hi)%20-%20np.log10(_lo)%0A%0A%20%20%20%20%23%20Ground-truth%20moment%20coordinates%2C%20mapped%20from%20the%20generative%0A%20%20%20%20%23%20(k_on%2C%20k_off%2C%20r_hat)%20with%20the%20SAME%20definitions%20scribe%20samples%20in%20(see%0A%20%20%20%20%23%20TwoStateMomentDeltaParameterization).%20For%20p%20~%20Beta(k_on%2C%20k_off)%20and%0A%20%20%20%20%23%20count%20~%20Poisson(r_hat%20*%20p)%2C%20with%20kappa%20%3D%20k_on%20%2B%20k_off%3A%0A%20%20%20%20%23%20%20%20mu%20%20%20%20%20%20%20%20%20%20%3D%20r_hat%20*%20k_on%20%2F%20kappa%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(%3D%20E%5Bcount%5D)%0A%20%20%20%20%23%20%20%20excess_fano%20%3D%20r_hat%20*%20k_off%20%2F%20(kappa*(kappa%2B1))%20%20%20%20(%3D%20Var%2FMean%20-%201)%0A%20%20%20%20%23%20%20%20delta%20%20%20%20%20%20%20%3D%201%20%2F%20(kappa%20%2B%201)%0A%20%20%20%20_kappa%20%3D%20k_on_true%20%2B%20k_off_true%0A%20%20%20%20mu_true%20%3D%20r_hat_true%20*%20k_on_true%20%2F%20_kappa%0A%20%20%20%20fano_true%20%3D%20r_hat_true%20*%20k_off_true%20%2F%20(_kappa%20*%20(_kappa%20%2B%201.0))%0A%20%20%20%20delta_true%20%3D%201.0%20%2F%20(_kappa%20%2B%201.0)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20delta_fit%2C%0A%20%20%20%20%20%20%20%20delta_true%2C%0A%20%20%20%20%20%20%20%20fano_fit%2C%0A%20%20%20%20%20%20%20%20fano_true%2C%0A%20%20%20%20%20%20%20%20koff_fit%2C%0A%20%20%20%20%20%20%20%20koff_logwidth%2C%0A%20%20%20%20%20%20%20%20kon_fit%2C%0A%20%20%20%20%20%20%20%20mu_fit%2C%0A%20%20%20%20%20%20%20%20mu_true%2C%0A%20%20%20%20%20%20%20%20rhat_fit%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20delta_fit%2C%0A%20%20%20%20delta_true%2C%0A%20%20%20%20fano_fit%2C%0A%20%20%20%20fano_true%2C%0A%20%20%20%20k_off_true%2C%0A%20%20%20%20k_on_true%2C%0A%20%20%20%20koff_fit%2C%0A%20%20%20%20kon_fit%2C%0A%20%20%20%20mu_fit%2C%0A%20%20%20%20mu_true%2C%0A%20%20%20%20n_bursty%2C%0A%20%20%20%20n_nb%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20r_hat_true%2C%0A%20%20%20%20rhat_fit%2C%0A)%3A%0A%20%20%20%20_colors%20%3D%20%5B%22C1%22%5D%20*%20n_bursty%20%2B%20%5B%22C0%22%5D%20*%20n_nb%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(2%2C%203%2C%20figsize%3D(12%2C%208))%0A%20%20%20%20_panels%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%23%20Top%20row%20%E2%80%94%20DERIVED%20biophysical%20rates.%20Read%20these%20for%20honesty%3A%20k%5E-%20and%0A%20%20%20%20%20%20%20%20%23%20r_hat%20are%20singular%20in%20the%20NB%20limit%20and%20will%20scatter%20for%20NB-like%20genes.%0A%20%20%20%20%20%20%20%20(k_on_true%2C%20kon_fit%2C%20%22%24k%5E%2B%24%20(burst%20frequency)%22)%2C%0A%20%20%20%20%20%20%20%20(k_off_true%2C%20koff_fit%2C%20%22%24k%5E-%24%20(OFF%20rate)%22)%2C%0A%20%20%20%20%20%20%20%20(r_hat_true%2C%20rhat_fit%2C%20r%22%24%5Chat%20r%24%20(ON-mode%20rate)%22)%2C%0A%20%20%20%20%20%20%20%20%23%20Bottom%20row%20%E2%80%94%20SAMPLED%20moment%20coordinates.%20These%20are%20what%20the%20data%0A%20%20%20%20%20%20%20%20%23%20directly%20constrains%3B%20all%20three%20land%20on%20the%20diagonal%20for%20every%20gene.%0A%20%20%20%20%20%20%20%20(mu_true%2C%20mu_fit%2C%20r%22%24%5Cmu%24%20(mean%20expression)%22)%2C%0A%20%20%20%20%20%20%20%20(fano_true%2C%20fano_fit%2C%20%22excess%20Fano%20%24F-1%24%22)%2C%0A%20%20%20%20%20%20%20%20(delta_true%2C%20delta_fit%2C%20r%22%24%5Cdelta%24%20(regime)%22)%2C%0A%20%20%20%20%5D%0A%20%20%20%20for%20_ax%2C%20(_t%2C%20_f%2C%20_name)%20in%20zip(_axes.ravel()%2C%20_panels)%3A%0A%20%20%20%20%20%20%20%20_ax.scatter(_t%2C%20_f%2C%20c%3D_colors%2C%20s%3D25)%0A%20%20%20%20%20%20%20%20_lo%20%3D%20min(_t.min()%2C%20_f.min())%0A%20%20%20%20%20%20%20%20_hi%20%3D%20max(_t.max()%2C%20_f.max())%0A%20%20%20%20%20%20%20%20_ax.plot(%5B_lo%2C%20_hi%5D%2C%20%5B_lo%2C%20_hi%5D%2C%20%22k--%22%2C%20lw%3D1)%0A%20%20%20%20%20%20%20%20_ax.set_xscale(%22log%22)%0A%20%20%20%20%20%20%20%20_ax.set_yscale(%22log%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(f%22true%20%7B_name%7D%22)%0A%20%20%20%20%20%20%20%20_ax.set_ylabel(f%22fitted%20%7B_name%7D%22)%0A%20%20%20%20_fig.suptitle(%0A%20%20%20%20%20%20%20%20%22Parameter%20recovery%20(orange%20%3D%20bursty%2C%20blue%20%3D%20NB-like)%5Cn%22%0A%20%20%20%20%20%20%20%20%22top%3A%20derived%20biophysical%20rates%20%20%20%20bottom%3A%20sampled%20moment%20coordinates%22%0A%20%20%20%20)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(delta_fit%2C%20k_off_true%2C%20koff_fit%2C%20koff_logwidth%2C%20n_bursty%2C%20np)%3A%0A%20%20%20%20%23%20Quantify%20identifiability%20per%20group%3A%20how%20well%20is%20k_off%20recovered%2C%20and%0A%20%20%20%20%23%20how%20wide%20is%20its%20posterior%3F%0A%20%20%20%20for%20_lbl%2C%20_sl%20in%20%5B%0A%20%20%20%20%20%20%20%20(%22BURSTY%20%20(truth%20k_off%20%3C%201)%22%2C%20slice(0%2C%20n_bursty))%2C%0A%20%20%20%20%20%20%20%20(%22NB-LIKE%20(truth%20k_off%20%3E%3E%201)%22%2C%20slice(n_bursty%2C%20None))%2C%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20print(f%22---%20%7B_lbl%7D%20---%22)%0A%20%20%20%20%20%20%20%20print(f%22%20%20k_off%20recovery%20%20%20median(fit%2Ftrue)%20%3D%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%7Bnp.median(koff_fit%5B_sl%5D%20%2F%20k_off_true%5B_sl%5D)%3A.2f%7D%22)%0A%20%20%20%20%20%20%20%20print(f%22%20%20k_off%20posterior%20%20median%20log10-IQR%20%3D%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%7Bnp.median(koff_logwidth%5B_sl%5D)%3A.2f%7D%20%20%20(%3E1%20spans%20%3E10x%20%3D%3E%20non-identified)%22)%0A%20%20%20%20%20%20%20%20print(f%22%20%20delta%20(sampled%20regime)%20%20median%20%3D%20%7Bnp.median(delta_fit%5B_sl%5D)%3A.3f%7D%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20What%20is%20identifiable%2C%20and%20what%20is%20not%0A%0A%20%20%20%20The%20recovery%20plot%20has%20**two%20rows**%2C%20and%20reading%20them%20together%20is%20the%20most%20important%20methodological%20lesson%20of%20this%20model.%20The%20**top%20row**%20shows%20the%20*derived*%20biophysical%20rates%20%24(k%5E%2B%2C%20k%5E-%2C%20%5Chat%20r)%24%3B%20the%20**bottom%20row**%20shows%20the%20coordinates%20the%20model%20actually%20*samples*%20%24(%5Cmu%2C%5C%2C%20F%5C!-%5C!1%2C%5C%2C%20%5Cdelta)%24.%20They%20are%20smooth%20functions%20of%20one%20another%20%E2%80%94%20the%20*same*%20posterior%20%E2%80%94%20yet%20the%20top%20row%20can%20look%20badly%20degenerate%20while%20the%20bottom%20row%20is%20far%20better%20behaved%3A%20two%20of%20its%20three%20panels%20are%20pinned%20for%20every%20gene%2C%20and%20the%20third%20fails%20*gracefully*%20(a%20bounded%20bias)%20rather%20than%20blowing%20up.%20That%20gap%20is%20the%20whole%20point%20%E2%80%94%20identifiability%20is%20a%20property%20of%20the%20*coordinate*%20you%20read%2C%20not%20of%20the%20fit.%0A%0A%20%20%20%20**Pinned%20for%20every%20gene%20%E2%80%94%20the%20observable%20moments%20%24%5Cmu%24%20and%20%24F-1%24%20(and%20the%20burst%20frequency%20%24k%5E%2B%24).**%20The%20mean%20%24%5Cmu%24%20(bottom-left)%20and%20the%20excess%20Fano%20%24F-1%24%20(bottom-middle)%20sit%20on%20the%20diagonal%20across%20the%20*entire*%20range%20%E2%80%94%20bursty%20and%20NB-like%20genes%20alike%20%E2%80%94%20because%20they%20are%20pure%20functions%20of%20the%20observed%20count%20distribution%2C%20so%20the%20data%20nail%20them%20directly.%20The%20burst%20frequency%20%24k%5E%2B%24%20(top-left)%20is%20recovered%20too.%20These%20are%20the%20directions%20you%20can%20trust%20per%20gene%2C%20and%20they%20are%20exactly%20the%20scale-free%2C%20observable-moment%20directions.%0A%0A%20%20%20%20**Non-identifiable%20in%20the%20NB%20limit%20%E2%80%94%20the%20OFF%20rate%20%24k%5E-%24%20and%20the%20exact%20regime%20%24%5Cdelta%24.**%20For%20the%20bursty%20genes%20both%20are%20pinned%3A%20%24k%5E-%24%20lands%20on%20the%20diagonal%20(top-middle)%20and%20so%20does%20%24%5Cdelta%24%20(bottom-right).%20For%20the%20NB-like%20genes%20*neither*%20does%20%E2%80%94%20and%20that%20is%20*correct*%2C%20not%20a%20fit%20failure.%20Once%20switching%20is%20fast%20the%20counts%20are%20negative-binomial%2C%20with%20only%20two%20degrees%20of%20freedom%20(burst%20frequency%20and%20burst%20size)%2C%20so%20the%20data%20carry%20no%20information%20about%20how%20*deep*%20into%20the%20NB%20limit%20a%20gene%20sits.%20The%20two%20coordinates%20express%20that%20same%20blind%20spot%20differently.%20The%20*unbounded*%20%24k%5E-%24%20shows%20it%20as%20a%20**blow-up**%3A%20it%20scatters%20over%20a%20decade%20toward%20%24%5Cinfty%24%20(the%20%60log10-IQR%60%20printed%20above%20is%20large%20for%20the%20NB-like%20group%2C%20tiny%20for%20the%20bursty%20group).%20The%20*bounded*%20%24%5Cdelta%20%3D%201%2F(%5Ckappa%2B1)%5Cin(0%2C1)%24%20cannot%20run%20to%20%24%5Cinfty%24%2C%20so%20the%20NB-like%20genes%20instead%20pile%20up%20just%20*above*%20the%20diagonal%20at%20small%20%24%5Cdelta%24%20%E2%80%94%20biased%20toward%20%22not-quite-as-deep-in-the-NB-limit-as-truth%2C%22%20held%20off%20the%20boundary%20by%20the%20prior.%20The%20payoff%20for%20practice%20is%20decisive%3A%20even%20un-pinned%2C%20%24%5Cdelta%24%20keeps%20the%20two%20regimes%20**cleanly%20separated**%20(here%20every%20NB-like%20gene%20has%20%24%5Cdelta%3C0.15%24%2C%20every%20bursty%20gene%20%24%5Cdelta%3E0.4%24)%2C%20so%20the%20bursty-vs-NB%20*classification*%20is%20robust%20even%20where%20the%20precise%20rate%20is%20not.%0A%0A%20%20%20%20**Built%20from%20the%20identified%20moments%2C%20but%20singular%20in%20the%20NB%20limit%20%E2%80%94%20the%20absolute%20rate%20%24%5Chat%20r%24%20(and%20burst%20size).**%20Compare%20the%20top-right%20(%24%5Chat%20r%24)%20and%20bottom-left%20(%24%5Cmu%24)%20panels%20for%20the%20*same*%20genes%3A%20%24%5Cmu%24%20sits%20on%20the%20diagonal%20everywhere%2C%20yet%20%24%5Chat%20r%24%20scatters%20for%20the%20NB-like%20genes.%20There%20is%20no%20contradiction%20%E2%80%94%20%24%5Chat%20r%20%3D%20%5Cmu%20%2B%20(F-1)%2F%5Cdelta%24%2C%20so%20as%20%24%5Cdelta%5Cto0%24%20(the%20NB%20limit)%20%24%5Chat%20r%5Cpropto%201%2F%5Cdelta%24%20diverges.%20The%20data%20pin%20%24%5Cmu%24%2C%20%24F-1%24%2C%20and%20the%20regime%2C%20but%20not%20the%20*exact*%20tiny%20%24%5Cdelta%24%2C%20and%20%24%5Chat%20r%24%20amplifies%20precisely%20that%20weakly-constrained%20direction%20(the%20same%20ridge%20as%20%24k%5E-%24).%20Only%20the%20finite%20combination%20%E2%80%94%20the%20NB%20shape%20%24r_%7B%5Cmathrm%7BNB%7D%7D%20%3D%20%5Cmu%2F(F-1)%24%20%E2%80%94%20is%20identified%20for%20those%20genes.%20This%20is%20the%20bottom-right%20bias%20seen%20from%20the%20other%20side%3A%20the%20NB-like%20genes%20whose%20%24%5Cdelta%24%20sits%20*above*%20its%20diagonal%20have%20%24%5Chat%20r%24%20*below*%20its%20diagonal%2C%20because%20%24%5Chat%20r%24%20falls%20as%20%24%5Cdelta%24%20rises%20%E2%80%94%20one%20fact%2C%20two%20coordinates.%20A%20*second*%2C%20independent%20scale%20ambiguity%20enters%20on%20**real**%20data%3A%20with%20variable%20capture%20(Phase%202)%20counts%20constrain%20only%20the%20*product*%20%24%5Chat%20r_g%5C%2C%5Cnu_c%24%20%E2%80%94%20rescale%20every%20gene's%20%24%5Chat%20r%24%20and%20inversely%20rescale%20every%20cell's%20%24%5Cnu_c%24%20and%20the%20likelihood%20is%20unchanged%20%E2%80%94%20so%20the%20absolute%20transcript%20rate%20rides%20on%20the%20capture%20*prior*%2C%20not%20the%20data.%20If%20you%20need%20it%20absolutely%2C%20anchor%20capture%20with%20an%20informative%20prior%20(the%20biology-informed%20capture%20anchor%20from%20the%20previous%20tutorial).%0A%0A%20%20%20%20**The%20one%20tension%20no%20parameterization%20escapes%20%E2%80%94%20valley%20depth%20vs.%20ON-population%20breadth.**%20A%20single%20Beta%20couples%20*how%20deep%20the%20valley%20is*%20and%20*how%20broad%20the%20ON%20population%20is*.%20A%20deep%20valley%20needs%20an%20extreme%20U-shaped%20Beta%2C%20which%20forces%20the%20ON%20cells%20to%20a%20narrow%20Poisson%3B%20a%20broad%20ON%20needs%20a%20milder%20Beta%2C%20which%20fills%20the%20valley.%20A%20real%20**marker%20gene**%20%E2%80%94%20bimodal%20because%20of%20discrete%20cell%20types%2C%20with%20an%20*over-dispersed*%20expressing%20population%20%E2%80%94%20needs%20both%20at%20once%2C%20and%20no%20single%20Beta%20delivers.%20That%20is%20a%20job%20for%20a%20**mixture%20model**%2C%20and%20it%20is%20exactly%20why%20we%20insist%20on%20a%20monoculture%3A%20to%20study%20bursting%20we%20must%20remove%20the%20cell-type%20axis%20the%20two-state%20model%20was%20never%20built%20to%20represent.%0A%0A%20%20%20%20With%20that%20understanding%20in%20hand%20%E2%80%94%20the%20observed%20moments%20%24%5Cmu%24%20and%20%24F-1%24%2C%20the%20burst%20frequency%20%24k%5E%2B%24%2C%20and%20the%20bursty-vs-NB%20*regime*%20are%20what%20we%20trust%20per%20gene%3B%20the%20exact%20switching%20rates%20(%24k%5E-%24%2C%20and%20%24%5Cdelta%24%20deep%20in%20the%20NB%20limit)%20and%20the%20absolute%20rate%20%24%5Chat%20r%24%20are%20*not*%20pinned%20by%20total-count%20data%3B%20and%20cell-type%20bimodality%20is%20out%20of%20scope%20%E2%80%94%20we%20turn%20to%20real%20data.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Phase%202%20%E2%80%94%20A%20real%20deeply-sequenced%20monoculture%20(K562)%0A%0A%20%20%20%20K562%20is%20a%20human%20chronic-myelogenous-leukemia%20line%20%E2%80%94%20a%20**monoculture**%2C%20so%20any%20bimodality%20we%20find%20is%20*within*%20one%20cell%20type%2C%20not%20a%20mix%20of%20types.%20We%20use%20a%20deeply-sequenced%2010x%20dataset%20(%5B10k%20K562-r%2C%20singleplex%5D(https%3A%2F%2Fwww.10xgenomics.com%2Fdatasets%2F10k-human-k562-r-cells-singleplex-sample-1-standard)).%0A%0A%20%20%20%20Set%20expectations%20honestly%20before%20we%20look.%20Genuine%20telegraph-bursting%20bimodality%20is%20**hard%20to%20see%20in%20droplet%20total%20counts**%2C%20for%20a%20fundamental%20reason%3A%20each%20diploid%20cell's%20count%20is%20the%20*sum%20of%20two%20independently-bursting%20alleles*%2C%20and%20the%20sum%20of%20two%20ON%2FOFF%20processes%20is%20far%20less%20bimodal%20than%20either%20allele%20(the%20mixed%20ON%2FOFF%20cells%20fill%20the%20valley).%20Capture%20noise%20compounds%20it%20by%20smearing%20the%20OFF%20state%20into%20the%20dropout%20floor.%20This%20is%20exactly%20why%20the%20bursting-inference%20literature%20works%20at%20the%20**allele%20level**%20%E2%80%94%20total%20counts%20are%20a%20blunt%20observable%20for%20switching.%20So%20we%20*expect*%20most%20genes%20here%20to%20sit%20in%20the%20negative-binomial%20limit%2C%20which%20the%20two-state%20model%20should%20report%20as%20such.%20The%20interesting%20question%20is%20whether%20the%20handful%20of%20genes%20it%20flags%20as%20deviating%20make%20biological%20sense.%0A%0A%20%20%20%20We%20load%20the%20data%20through%20%60scribe.data_loader%60%2C%20as%20before.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Path%2C%20scribe)%3A%0A%20%20%20%20%23%20Load%20the%20local%20tutorial%20path%20configuration.%0A%20%20%20%20import%20json%0A%20%20%20%20with%20Path(__file__).with_name(%22tutorial_paths.local.json%22).open(%0A%20%20%20%20%20%20%20%20%22r%22%2C%20encoding%3D%22utf-8%22%0A%20%20%20%20)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20_data_root%20%3D%20json.load(_f)%5B%22SCRIBE_TUTORIAL_DATA_ROOT%22%5D%0A%0A%20%20%20%20%23%20Build%20the%20dataset-specific%20path%20relative%20to%20the%20configured%20root.%0A%20%20%20%20data_dir%20%3D%20Path(_data_root).expanduser()%20%2F%20%2210k_k562r%22%0A%0A%20%20%20%20%23%20Load%20the%20data%0A%20%20%20%20adata_total%20%3D%20scribe.data_loader.load_and_preprocess_anndata(%0A%20%20%20%20%20%20%20%20data_dir%2C%20return_jax%3DFalse%0A%20%20%20%20)%0A%0A%20%20%20%20adata_total%0A%20%20%20%20return%20adata_total%2C%20data_dir%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Before%20fitting%20anything%2C%20let%20us%20look%20at%20the%20per-cell%20library%20size%20(total%20UMI%20count).%20This%20is%20the%20depth%20axis%20from%20the%20discussion%20above%3A%20it%20tells%20us%20how%20much%20of%20each%20cell's%20transcriptome%20actually%20made%20it%20into%20the%20data%2C%20and%20therefore%20how%20much%20room%20we%20have%20to%20resolve%20a%20true%20OFF%20state%20from%20a%20technical%20zero.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_total%2C%20np%2C%20plt%2C%20sns)%3A%0A%20%20%20%20%23%20Total%20UMI%20per%20cell%3B%20.A1%20flattens%20a%20sparse-matrix%20row%20sum%20to%20a%201-D%20array%0A%20%20%20%20_library_sizes%20%3D%20np.asarray(adata_total.X.sum(axis%3D1)).ravel()%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(5%2C%204))%0A%20%20%20%20sns.ecdfplot(_library_sizes%2C%20ax%3D_ax)%0A%20%20%20%20_ax.set_xlabel('total%20UMI%20count%20per%20cell')%0A%20%20%20%20_ax.set_ylabel('ECDF')%0A%20%20%20%20_ax.set_ylim(-0.01%2C%201.01)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20library-size%20distribution%20spans%20a%20wide%20range%2C%20and%20many%20barcodes%20are%20shallow.%20Since%20depth%20is%20precisely%20what%20lets%20bimodality%20survive%20into%20the%20counts%2C%20we%20keep%20only%20the%20**deeply-sequenced%20cells**%20%E2%80%94%20those%20with%20at%20least%20%60umi_thresh%60%20total%20UMIs.%20This%20is%20not%20the%20usual%20quality-control%20reflex%20of%20discarding%20%22bad%22%20cells%3B%20it%20is%20a%20deliberate%20enrichment%20for%20the%20regime%20where%20the%20two-state%20signal%20is%20detectable.%20The%20trade-off%20is%20fewer%20cells%2C%20but%20each%20carries%20more%20information%20about%20the%20*shape*%20of%20every%20gene's%20distribution%20%E2%80%94%20which%20is%20exactly%20what%20we%20are%20trying%20to%20read.%20The%20ECDF%20below%20shows%20the%20library-size%20distribution%20after%20the%20cut.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_total%2C%20np%2C%20plt%2C%20sns)%3A%0A%20%20%20%20%23%20Define%20UMI%20count%20threshold%0A%20%20%20%20umi_thresh%20%3D%2010_000%0A%0A%20%20%20%20%23%20Filter%20data%0A%20%20%20%20adata%20%3D%20adata_total%5Badata_total.X.sum(axis%3D1)%20%3E%3D%20umi_thresh%5D%0A%0A%20%20%20%20%23%20Total%20UMI%20per%20cell%3B%20.A1%20flattens%20a%20sparse-matrix%20row%20sum%20to%20a%201-D%20array%0A%20%20%20%20_library_sizes%20%3D%20np.asarray(adata.X.sum(axis%3D1)).ravel()%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(5%2C%204))%0A%20%20%20%20sns.ecdfplot(_library_sizes%2C%20ax%3D_ax)%0A%20%20%20%20_ax.set_xlabel('total%20UMI%20count%20per%20cell')%0A%20%20%20%20_ax.set_ylabel('ECDF')%0A%20%20%20%20_ax.set_ylim(-0.01%2C%201.01)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20adata%2C%20umi_thresh%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Fitting%20the%20negative%20binomial%20foil%0A%0A%20%20%20%20Our%20point%20of%20comparison%20is%20a%20**gene-specific%20%24p_g%24**%20negative%20binomial%20%E2%80%94%20the%20most%20flexible%20negative%20binomial%20from%20the%20previous%20tutorial%2C%20with%20each%20gene%20drawing%20its%20own%20success%20probability%20(%60prob_prior%3D%22gaussian%22%60%2C%20i.e.%20a%20learned%20%24%5Cmathrm%7Blogit%7D(p_g)%5Csim%5Cmathcal%7BN%7D(%5Cmu_p%2C%5Csigma_p)%24%20hierarchy)%20on%20top%20of%20per-cell%20variable%20capture%2C%20fit%20in%20the%20%60mean_odds%60%20parameterization.%20Because%20we%20are%20on%20a%20new%20dataset%2C%20we%20refit%20it%20from%20scratch%20rather%20than%20reusing%20the%20previous%20tutorial's%20pickle.%0A%0A%20%20%20%20Why%20this%20particular%20foil%3F%20Because%20it%20makes%20the%20punchline%20sharp.%20Per-gene%20%24p_g%24%20is%20the%20most%20a%20negative%20binomial%20can%20flex%20within%20its%20own%20family%2C%20and%20yet%20%E2%80%94%20being%20a%20negative%20binomial%20%E2%80%94%20it%20is%20*still%20structurally%20unimodal*.%20If%20the%20two-state%20model%20captures%20bimodality%20that%20this%20model%20misses%2C%20the%20gap%20cannot%20be%20blamed%20on%20%22the%20negative%20binomial%20just%20needed%20more%20parameters.%22%20It%20is%20a%20difference%20in%20*shape*%2C%20not%20in%20*flexibility*.%0A%0A%20%20%20%20%23%23%23%20One%20transform%20for%20everything%3A%20%60positive_transform%3D%22exp%22%60%0A%0A%20%20%20%20There%20is%20one%20change%20from%20the%20previous%20tutorial%20worth%20dwelling%20on%2C%20because%20it%20bit%20**both**%20models%20here.%20%60scribe%60%20fits%20in%20an%20unconstrained%20space%20and%20maps%20back%20to%20the%20positive%20parameters%20through%20a%20transform%3B%20the%20default%20is%20%60softplus%60.%20On%20the%20shallower%20Jurkat%20data%2C%20%60softplus%60%20was%20fine%20for%20everything%20except%20the%20mean%20(where%20we%20used%20%60exp%60).%20On%20a%20deeply-sequenced%20dataset%20like%20K562%2C%20*every*%20positive%20parameter%20needed%20%60exp%60%20%E2%80%94%20hence%20%60positive_transform%3D%22exp%22%60%20applied%20globally%20to%20both%20the%20negative%20binomial%20and%20the%20two-state%20fit%20(and%20to%20the%20synthetic%20fits%20in%20Phase%201).%0A%0A%20%20%20%20The%20reason%20is%20the%20depth%20itself.%20A%20deeply-sequenced%20library%20stretches%20the%20true%20dynamic%20range%20of%20the%20parameters%20across%20many%20orders%20of%20magnitude%3A%20gene%20means%2C%20dispersions%2C%20burst%20sizes%20and%20switching%20rates%20now%20genuinely%20span%20from%20%24%5Csim%5C!10%5E%7B-2%7D%24%20to%20%24%5Csim%5C!10%5E%7B3%7D%24.%20%60softplus%60%20is%20nearly%20linear%20for%20large%20arguments%2C%20so%20its%20Jacobian%20*saturates*%20%E2%80%94%20to%20move%20a%20parameter%20from%20%241%24%20to%20%241000%24%20the%20optimizer%20must%20take%20an%20enormous%20unconstrained%20step%2C%20and%20gradient%20descent%20crawls.%20%60exp%60%20makes%20the%20step%20*multiplicative*%3A%20a%20fixed%20step%20in%20unconstrained%20space%20is%20a%20fixed%20*ratio*%20in%20parameter%20space%2C%20so%20the%20optimizer%20traverses%20several%20decades%20as%20easily%20as%20one.%20When%20the%20data%20occupy%20only%20a%20narrow%20range%20(a%20shallow%20library)%2C%20%60softplus%60%20is%20perfectly%20adequate%20and%20you%20never%20notice%3B%20when%20the%20data%20reveal%20the%20full%20log-scale%20spread%2C%20%60exp%60%20becomes%20necessary.%0A%0A%20%20%20%20This%20is%20worth%20stating%20plainly%3A%20**there%20is%20no%20single%2C%20monolithic%20configuration%20of%20these%20models%20that%20fits%20every%20dataset.**%20The%20right%20parameterization%2C%20transform%2C%20optimizer%2C%20and%20depth%20filter%20all%20depend%20on%20what%20the%20data%20actually%20contain.%20That%20is%20not%20a%20defect%20of%20the%20approach%20%E2%80%94%20it%20is%20the%20reality%20of%20probabilistic%20modeling%20and%20optimization%20at%20this%20scale.%20The%20workflow%20is%20iterative%20by%20design%3A%20fit%2C%20read%20the%20diagnostics%2C%20adjust%20the%20piece%20they%20flag%2C%20refit.%20The%20diagnostics%20we%20run%20next%20are%20exactly%20how%20you%20catch%20a%20saturating%20transform%20or%20an%20under-converged%20fit%20*before*%20trusting%20any%20downstream%20biology.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20clear_caches%2C%20data_dir%2C%20gc%2C%20pickle%2C%20scribe%2C%20umi_thresh)%3A%0A%20%20%20%20%23%20Output%20directory%20shared%20with%20the%20previous%20tutorial's%20fits.%0A%20%20%20%20out_dir%20%3D%20data_dir%20%2F%20%22scribe_results%22%0A%20%20%20%20out_dir.mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20%23%20Clear%20GPU%20memory%20from%20previous%20fits%0A%20%20%20%20clear_caches()%0A%20%20%20%20gc.collect()%0A%0A%20%20%20%20%23%20Define%20parameterization%0A%20%20%20%20_parameterization%20%3D%20%22mean_odds%22%0A%0A%20%20%20%20%23%20Define%20output%20file%20path%0A%20%20%20%20_out_path%20%3D%20out_dir%20%2F%20f%22scribe_results_nbvcp-prob_%7B_parameterization%7D_%7Bumi_thresh%7Dumi.pkl%22%0A%0A%20%20%20%20if%20_out_path.exists()%3A%0A%20%20%20%20%20%20%20%20%23%20Load%20model%20from%20pkl%20file%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22rb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20results_nb%20%3D%20pickle.load(_f)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%23%20Fit%20basic%20model%20to%20data%20with%20variable%20capture%20probability%20fit%20per%20cell%0A%20%20%20%20%20%20%20%20results_nb%20%3D%20scribe.fit(%0A%20%20%20%20%20%20%20%20%20%20%20%20adata%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20variable_capture%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20parameterization%3D_parameterization%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20prob_prior%3D%22gaussian%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20unconstrained%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20positive_transform%3D%22exp%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_steps%3D100_000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20batch_size%3D2048%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_quad_nodes%3D256%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%23%20Save%20the%20fitted%20model%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22wb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pickle.dump(results_nb%2C%20_f)%0A%0A%20%20%20%20results_nb%0A%20%20%20%20return%20out_dir%2C%20results_nb%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20first%20diagnostic%2C%20as%20always%2C%20is%20the%20ELBO%20loss%20curve%3A%20a%20healthy%20run%20drops%20fast%20and%20then%20flattens.%20If%20the%20curve%20were%20still%20sliding%20downward%20at%20the%20end%20%E2%80%94%20a%20classic%20symptom%20of%20a%20saturating%20transform%20starving%20the%20optimizer%20%E2%80%94%20that%20would%20be%20the%20cue%20to%20revisit%20the%20configuration%20above%20before%20reading%20anything%20into%20the%20fit.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(results_nb%2C%20scribe)%3A%0A%20%20%20%20%23%20Plot%20ELBO%20loss%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_loss(results_nb%2C%20figsize%3D(6%2C%203))%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Next%2C%20posterior%20predictive%20checks%20across%20a%20grid%20of%20genes%20spanning%20the%20dynamic%20range.%20For%20the%20great%20majority%20of%20genes%20the%20negative%20binomial%20is%20an%20excellent%20model%20%E2%80%94%20it%20should%20track%20the%20observed%20histograms%20closely%20here.%20That%20is%20the%20point%20of%20the%20foil%3A%20it%20is%20*not*%20a%20strawman.%20Keep%20this%20in%20mind%20for%20later%2C%20because%20the%20genes%20where%20it%20fails%20are%20a%20specific%2C%20mechanistically%20meaningful%20minority%2C%20not%20a%20sign%20of%20general%20inadequacy.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20results_nb%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_nb%2C%0A%20%20%20%20%20%20%20%20adata%2C%0A%20%20%20%20%20%20%20%20n_genes%3D16%2C%0A%20%20%20%20%20%20%20%20n_rows%3D4%2C%0A%20%20%20%20%20%20%20%20figsize%3D(8%2C%208)%2C%0A%20%20%20%20%20%20%20%20n_samples%3D512%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20A%20global%20check%20complements%20the%20per-gene%20PPCs%3A%20the%20mean-calibration%20plot%20compares%20each%20gene's%20observed%20mean%20UMI%20count%20to%20the%20model's%20predicted%20mean%20from%20the%20MAP%20parameters.%20Points%20hugging%20the%20identity%20line%20mean%20the%20model%20reproduces%20overall%20expression%20levels%20across%20the%20whole%20transcriptome%2C%20not%20just%20the%20handful%20of%20genes%20that%20happen%20to%20land%20in%20the%20PPC%20grid.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20results_nb%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_mean_calibration(%0A%20%20%20%20%20%20%20%20results_nb%2C%20counts%3Dadata%2C%20figsize%3D(4%2C%204)%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Finally%2C%20the%20capture-scaling%20diagnostic%3A%20each%20cell's%20fitted%20capture%20probability%20%24%5Cnu_c%24%20against%20its%20library%20size.%20Because%20we%20deliberately%20kept%20only%20deep%20cells%2C%20this%20also%20confirms%20that%20the%20depth%20we%20selected%20for%20is%20being%20absorbed%20by%20the%20technical%20capture%20channel%20rather%20than%20leaking%20into%20the%20gene-level%20parameters.%20With%20the%20negative%20binomial%20foil%20validated%2C%20we%20can%20fit%20the%20two-state%20model%20and%20compare.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20results_nb%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_p_capture_scaling(%0A%20%20%20%20%20%20%20%20results_nb%2C%20counts%3Dadata%2C%20figsize%3D(4%2C%204)%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Fitting%20the%20two-state%20model%0A%0A%20%20%20%20Now%20we%20fit%20the%20two-state%20promoter%20model%20with%20the%20%60moment_delta%60%20configuration%20described%20earlier%20%E2%80%94%20and%2C%20for%20the%20same%20reason%20we%20just%20discussed%2C%20with%20%60positive_transform%3D%22exp%22%60%20on%20every%20positive%20parameter.%20We%20clear%20GPU%20memory%20before%20the%20fit%20and%20cache%20the%20result%20to%20a%20pickle%20so%20the%20(GPU-only)%20fit%20runs%20once.%20The%20two-state%20likelihood%20is%20heavier%20than%20the%20negative%20binomial%20(it%20integrates%20over%20the%20Beta%20latent%20at%20every%20evaluation)%2C%20so%20this%20is%20the%20most%20expensive%20fit%20in%20the%20notebook.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20clear_caches%2C%20gc%2C%20out_dir%2C%20pickle%2C%20scribe%2C%20umi_thresh)%3A%0A%20%20%20%20%23%20Clear%20GPU%20memory%20from%20previous%20fits%0A%20%20%20%20clear_caches()%0A%20%20%20%20gc.collect()%0A%0A%20%20%20%20%23%20Define%20the%20fit%20identity%20(used%20for%20the%20cache%20filename)%0A%20%20%20%20_model_type%20%3D%20%22twostatevcp%22%0A%20%20%20%20_parameterization%20%3D%20%22moment_delta%22%0A%20%20%20%20_inference_method%20%3D%20%22svi%22%0A%0A%20%20%20%20_out_path%20%3D%20out_dir%20%2F%20(%0A%20%20%20%20%20%20%20%20f%22scribe_results_%7B_model_type%7D-%22%0A%20%20%20%20%20%20%20%20f%22%7B_parameterization%7D-%22%0A%20%20%20%20%20%20%20%20f%22%7Bumi_thresh%7Dumi.pkl%22%0A%20%20%20%20)%0A%0A%20%20%20%20if%20_out_path.exists()%3A%0A%20%20%20%20%20%20%20%20%23%20Load%20model%20from%20pkl%20file%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22rb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20results_twostate%20%3D%20pickle.load(_f)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%23%20Fit%20the%20two-state%20promoter%20model%20with%20per-cell%20variable%20capture.%0A%20%20%20%20%20%20%20%20results_twostate%20%3D%20scribe.fit(%0A%20%20%20%20%20%20%20%20%20%20%20%20adata%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20model%3D_model_type%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20parameterization%3D_parameterization%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20inference_method%3D_inference_method%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20unconstrained%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20positive_transform%3D%22exp%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_steps%3D100_000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20batch_size%3D1024%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20priors%3D%7B%22inv_concentration%22%3A%20(0.0%2C%202.0)%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20optimizer_config%3D%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22name%22%3A%20%22clipped_adam%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22step_size%22%3A%201e-3%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22grad_clip_norm%22%3A%2010.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20early_stopping%3D%7B%22enabled%22%3A%20True%2C%20%22patience%22%3A%2010_000%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_quad_nodes%3D256%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%23%20Save%20the%20fitted%20model%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22wb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pickle.dump(results_twostate%2C%20_f)%0A%0A%20%20%20%20results_twostate%0A%20%20%20%20return%20(results_twostate%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20As%20always%2C%20the%20first%20plot%20to%20make%20is%20the%20ELBO%20loss%20curve%3A%20a%20clean%20run%20drops%20quickly%20and%20then%20flattens.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(results_twostate%2C%20scribe)%3A%0A%20%20%20%20%23%20Plot%20ELBO%20loss%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_loss(results_twostate%2C%20figsize%3D(6%2C%203))%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20results_twostate%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_ppc(%0A%0A%20%20%20%20%20%20%20%20results_twostate%2C%0A%20%20%20%20%20%20%20%20adata%2C%0A%20%20%20%20%20%20%20%20n_genes%3D16%2C%0A%20%20%20%20%20%20%20%20n_rows%3D4%2C%0A%20%20%20%20%20%20%20%20figsize%3D(8%2C%208)%2C%0A%20%20%20%20%20%20%20%20n_samples%3D512%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20fit%20to%20these%20example%20genes%20is%20great.%20That%20is%20a%20good%20sanity%20check%20for%20the%20two-state%20promoter%20model.%20Let's%20now%20look%20at%20the%20global%20diagnostic%20for%20the%20predicted%20and%20observed%20mean%20expression%20across%20all%20genes.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20results_twostate%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_mean_calibration(%0A%20%20%20%20%20%20%20%20results_twostate%2C%20counts%3Dadata%2C%20figsize%3D(4%2C%204)%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20mean%20calibration%20looks%20great.%20That%20means%20that%20that%20the%20combination%20of%20the%20two-state%20and%20the%20variable%20capture%20probability%20is%20able%20to%20reproduce%20the%20observed%20UMI%20counts.%0A%0A%20%20%20%20Because%20%60twostatevcp%60%20carries%20a%20per-cell%20capture%20parameter%20%24%5Cnu_c%24%2C%20the%20capture-scaling%20diagnostic%20from%20the%20previous%20tutorial%20applies%20unchanged%20%E2%80%94%20it%20confirms%20depth%20is%20being%20routed%20through%20the%20technical%20channel%20sensibly.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20results_twostate%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_p_capture_scaling(%0A%20%20%20%20%20%20%20%20results_twostate%2C%20counts%3Dadata%2C%20figsize%3D(4%2C%204)%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Reading%20the%20burstiness%20off%20the%20fit%0A%0A%20%20%20%20The%20two-state%20fit%20gives%20us%2C%20for%20every%20gene%2C%20a%20posterior%20over%20its%20biophysical%20parameters.%20The%20quantity%20that%20most%20directly%20measures%20%22how%20far%20from%20a%20negative%20binomial%22%20a%20gene%20is%20sits%20in%20the%20**OFF%20rate%20%24%5Chat%20k%5E-_g%24**%3A%0A%0A%20%20%20%20-%20%24%5Chat%20k%5E-_g%20%5Cgg%201%24%20%E2%80%94%20fast%20switching%2C%20the%20negative%20binomial%20limit.%20The%20gene%20is%20effectively%20constitutively%20bursty%3B%20its%20counts%20are%20unimodal.%0A%20%20%20%20-%20%24%5Chat%20k%5E-_g%20%5Clesssim%201%24%20%E2%80%94%20slow%20switching.%20The%20OFF%2FON%20dwell%20times%20are%20comparable%20to%20the%20mRNA%20lifetime%2C%20the%20lifetime%20cannot%20average%20over%20the%20switching%2C%20and%20the%20counts%20are%20genuinely%20**bimodal**.%0A%0A%20%20%20%20Even%20though%20we%20sampled%20in%20%60moment_delta%60%20coordinates%20%24(%5Cmu_g%2C%20F_g%2C%20%5Cdelta_g)%24%2C%20%60scribe%60%20materializes%20the%20derived%20biophysical%20parameters%20(%60k_off%60%2C%20%60k_on%60%2C%20%60burst_size%60%2C%20%60r_hat%60%2C%20%E2%80%A6)%20as%20deterministic%20sites%20in%20the%20posterior%2C%20so%20we%20can%20pull%20them%20directly%20with%20%60get_posterior_samples%60%20and%20take%20per-gene%20posterior%20medians.%0A%0A%20%20%20%20Here%20we%20fit%20both%20models%20on%20the%20full%20gene%20set%20(no%20%60gene_coverage%60%20pre-filter)%2C%20so%20the%20two-state%20and%20negative%20binomial%20panels%20share%20the%20same%20genes%20and%20line%20up%20directly%20%E2%80%94%20there%20is%20no%20pooled%20%60_other%60%20column%20to%20account%20for.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20np%2C%20results_twostate)%3A%0A%20%20%20%20%23%20Draw%20posterior%20samples%20of%20the%20(derived)%20biophysical%20parameters%20and%0A%20%20%20%20%23%20summarize%20each%20gene%20by%20its%20posterior%20median.%0A%20%20%20%20_post%20%3D%20results_twostate.get_posterior_samples(%0A%20%20%20%20%20%20%20%20n_samples%3D1_000%2C%20store_samples%3DFalse%2C%20convert_to_numpy%3DTrue%0A%20%20%20%20)%0A%0A%20%20%20%20k_off%20%3D%20np.median(_post%5B%22k_off%22%5D%2C%20axis%3D0)%0A%20%20%20%20excess_fano%20%3D%20np.median(_post%5B%22excess_fano%22%5D%2C%20axis%3D0)%0A%20%20%20%20mu_ts%20%3D%20np.median(_post%5B%22mu%22%5D%2C%20axis%3D0)%0A%20%20%20%20burst_size%20%3D%20np.median(_post%5B%22burst_size%22%5D%2C%20axis%3D0)%0A%0A%20%20%20%20gene_names_retained%20%3D%20np.asarray(adata.var_names)%0A%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22%7B(k_off%20%3C%201.0).sum()%7D%20have%20posterior-median%20k_off%20%3C%201%20%22%0A%20%20%20%20%20%20%20%20f%22(slow-switching%20%2F%20bursty).%22%0A%20%20%20%20)%0A%20%20%20%20return%20excess_fano%2C%20gene_names_retained%2C%20k_off%2C%20mu_ts%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Plotting%20the%20OFF%20rate%20%24k%5E-_g%24%20against%20the%20gene%20mean%20%24%5Cmu_g%24%2C%20coloured%20by%20the%20excess%20Fano%20factor%2C%20gives%20a%20map%20of%20the%20transcriptome's%20burstiness.%20Genes%20below%20the%20dashed%20line%20(%24k%5E-_g%20%3C%201%24)%20are%20the%20slow-switching%2C%20bimodal-capable%20genes%20%E2%80%94%20and%20they%20should%20also%20carry%20the%20largest%20excess%20overdispersion%2C%20since%20the%20two-state%20model%20attributes%20overdispersion%20to%20promoter%20switching%20(the%20excess%20variance%20is%20maximized%20at%20symmetric%20switching%20%24k%5E%2B_g%20%3D%20k%5E-_g%24%20and%20vanishes%20when%20either%20rate%20dominates).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(excess_fano%2C%20k_off%2C%20mu_ts%2C%20np%2C%20plt)%3A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(5.5%2C%204.5))%0A%20%20%20%20_sc%20%3D%20_ax.scatter(%0A%20%20%20%20%20%20%20%20mu_ts%2C%0A%20%20%20%20%20%20%20%20k_off%2C%0A%20%20%20%20%20%20%20%20c%3Dnp.log10(excess_fano%20%2B%201e-3)%2C%0A%20%20%20%20%20%20%20%20cmap%3D%22viridis%22%2C%0A%20%20%20%20%20%20%20%20s%3D8%2C%0A%20%20%20%20%20%20%20%20alpha%3D0.6%2C%0A%20%20%20%20%20%20%20%20edgecolor%3D%22none%22%2C%0A%20%20%20%20)%0A%20%20%20%20_ax.axhline(1.0%2C%20color%3D%22black%22%2C%20linestyle%3D%22--%22%2C%20lw%3D1)%0A%20%20%20%20_ax.set_xscale(%22log%22)%0A%20%20%20%20_ax.set_yscale(%22log%22)%0A%20%20%20%20_ax.set_xlabel(r%22gene%20mean%20expression%20%24%5Cmu_g%24%22)%0A%20%20%20%20_ax.set_ylabel(r%22OFF%20rate%20%24k%5E-_g%24%20(mRNA-lifetime%20units)%22)%0A%20%20%20%20_ax.text(%0A%20%20%20%20%20%20%20%20_ax.get_xlim()%5B0%5D%20*%201.5%2C%0A%20%20%20%20%20%20%20%200.6%2C%0A%20%20%20%20%20%20%20%20%22bursty%20%2F%20bimodal%20%20(%24k%5E-_g%20%3C%201%24)%22%2C%0A%20%20%20%20%20%20%20%20fontsize%3D9%2C%0A%20%20%20%20%20%20%20%20color%3D%22black%22%2C%0A%20%20%20%20)%0A%20%20%20%20_cb%20%3D%20plt.colorbar(_sc%2C%20ax%3D_ax)%0A%20%20%20%20_cb.set_label(r%22%24%5Clog_%7B10%7D%24%20excess%20Fano%20factor%20%24F_g%24%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Which%20genes%20leave%20the%20NB%20limit%20%E2%80%94%20and%20are%20they%20cell-cycle%20genes%3F%0A%0A%20%20%20%20On%20a%20monoculture%20we%20do%20not%20expect%20a%20forest%20of%20bimodal%20genes%20(the%20allele-averaging%20argument%20above)%2C%20so%20the%20question%20is%20sharper%3A%20of%20the%20genes%20the%20two-state%20model%20flags%20as%20*most*%20over-dispersed%20and%20slowest-switching%2C%20do%20any%20make%20biological%20sense%3F%20In%20a%20proliferating%20line%20like%20K562%20the%20most%20reliable%20source%20of%20genuine%20within-population%20spread%20is%20the%20**cell%20cycle**%20%E2%80%94%20histones%2C%20mitotic%20cyclins%20and%20kinetochore%20genes%20(%60HIST1H4C%60%2C%20%60CCNB1%60%2C%20%60CCNB2%60%2C%20%60MKI67%60%2C%20%60TOP2A%60%2C%20%60UBE2C%60%2C%20%60CENPF%60%2C%20%E2%80%A6)%20swing%20with%20cycle%20phase.%20If%20the%20model's%20%22bursty%22%20flags%20pick%20those%20up%2C%20that%20is%20a%20reassuring%20sanity%20check%20%E2%80%94%20with%20the%20honest%20caveat%20that%20the%20mechanism%20is%20cycle%20phase%2C%20not%20promoter%20telegraph.%0A%0A%20%20%20%20We%20pull%20the%20most-bursty%20genes%20(smallest%20%24%5Chat%20k%5E-_g%24)%2C%20note%20how%20many%20are%20cell-cycle%20markers%2C%20then%20run%20the%20head-to-head%20PPC%20on%20them%3A%20the%20gene-specific%20%24p_g%24%20negative%20binomial%20vs.%20the%20two-state%20model%2C%20the%20same%20genes%20by%20name%20in%20both%20panels%20(via%20%60plot_ppc%60's%20%60genes%3D%60%20argument).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gene_names_retained%2C%20k_off%2C%20mu_ts%2C%20np)%3A%0A%20%20%20%20%23%20Most%20bursty%20genes%20(smallest%20k_off)%20among%20adequately%20expressed%20genes.%0A%20%20%20%20_expressed%20%3D%20(mu_ts%20%3E%3D%201.0)%20%26%20(mu_ts%20%3C%3D%20500)%0A%20%20%20%20_cand%20%3D%20np.where(_expressed)%5B0%5D%0A%20%20%20%20_top%20%3D%20_cand%5Bnp.argsort(k_off%5B_cand%5D)%5D%5B%3A8%5D%0A%0A%20%20%20%20bursty_names%20%3D%20list(gene_names_retained%5B_top%5D)%0A%0A%20%20%20%20%23%20Cell-cycle%20markers%20we'd%20expect%20to%20be%20the%20most%20over-dispersed%20in%20a%0A%20%20%20%20%23%20proliferating%20monoculture%20(intersected%20with%20genes%20actually%20present).%0A%20%20%20%20_cc%20%3D%20%5B%22HIST1H4C%22%2C%20%22HIST1H1C%22%2C%20%22HIST1H2AC%22%2C%20%22CCNB1%22%2C%20%22CCNB2%22%2C%20%22MKI67%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%22TOP2A%22%2C%20%22UBE2C%22%2C%20%22CENPF%22%2C%20%22PCNA%22%2C%20%22TYMS%22%2C%20%22CKS1B%22%2C%20%22CKS2%22%5D%0A%20%20%20%20_name_to_idx%20%3D%20%7Bg%3A%20i%20for%20i%2C%20g%20in%20enumerate(gene_names_retained)%7D%0A%20%20%20%20_cc_here%20%3D%20%5Bg%20for%20g%20in%20_cc%20if%20g%20in%20_name_to_idx%5D%0A%20%20%20%20_cc_set%20%3D%20set(_cc_here)%0A%0A%20%20%20%20print(%22Most%20bursty%20genes%20(smallest%20k_off)%3A%22)%0A%20%20%20%20for%20_name%2C%20_ko%20in%20zip(bursty_names%2C%20k_off%5B_top%5D)%3A%0A%20%20%20%20%20%20%20%20_flag%20%3D%20%22%20%20%20%3C-%20cell-cycle%22%20if%20_name%20in%20_cc_set%20else%20%22%22%0A%20%20%20%20%20%20%20%20print(f%22%20%20%7B_name%3A%3E12s%7D%20%20%20k_off%20%3D%20%7B_ko%3A.3f%7D%7B_flag%7D%22)%0A%0A%20%20%20%20if%20_cc_here%3A%0A%20%20%20%20%20%20%20%20print(%22%5CnCell-cycle%20markers%20present%2C%20ranked%20by%20k_off%20(smaller%20%3D%20burstier)%3A%22)%0A%20%20%20%20%20%20%20%20for%20_g%20in%20sorted(_cc_here%2C%20key%3Dlambda%20g%3A%20k_off%5B_name_to_idx%5Bg%5D%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20print(f%22%20%20%7B_g%3A%3E12s%7D%20%20%20k_off%20%3D%20%7Bk_off%5B_name_to_idx%5B_g%5D%5D%3A.3f%7D%22)%0A%20%20%20%20return%20(bursty_names%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Negative%20binomial%20(gene-specific%20%24p_g%24).**%20Posterior%20predictive%20checks%20for%20the%20eight%20most%20bursty%20genes.%20Watch%20for%20the%20model's%20predictive%20band%20sitting%20*between*%20the%20data's%20two%20modes%20%E2%80%94%20a%20single%20hump%20trying%20to%20straddle%20a%20bimodal%20histogram.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20bursty_names%2C%20results_nb%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_nb%2C%0A%20%20%20%20%20%20%20%20adata%2C%0A%20%20%20%20%20%20%20%20genes%3Dbursty_names%2C%0A%20%20%20%20%20%20%20%20n_rows%3D2%2C%0A%20%20%20%20%20%20%20%20figsize%3D(9%2C%204)%2C%0A%20%20%20%20%20%20%20%20n_samples%3D512%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Two-state%20promoter%20model.**%20The%20same%20eight%20genes%2C%20in%20the%20same%20order.%20The%20Poisson%E2%80%93Beta%20likelihood%20can%20place%20mass%20at%20zero%20*and*%20at%20a%20second%20mode%20simultaneously%2C%20so%20its%20predictive%20band%20should%20track%20the%20bimodal%20shape%20the%20negative%20binomial%20could%20only%20average%20over.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20bursty_names%2C%20results_twostate%2C%20scribe)%3A%0A%20%20%20%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_twostate%2C%0A%20%20%20%20%20%20%20%20adata%2C%0A%20%20%20%20%20%20%20%20genes%3Dbursty_names%2C%0A%20%20%20%20%20%20%20%20n_rows%3D2%2C%0A%20%20%20%20%20%20%20%20figsize%3D(9%2C%204)%2C%0A%20%20%20%20%20%20%20%20n_samples%3D512%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20On%20K562%20the%20contrast%20is%20**muted%20compared%20with%20the%20synthetic%20showcase**%2C%20and%20that%20is%20the%20honest%2C%20expected%20result.%20Most%20genes%20sit%20so%20close%20to%20the%20negative-binomial%20limit%20that%20the%20two%20predictive%20bands%20are%20nearly%20indistinguishable%20%E2%80%94%20the%20two-state%20model%20is%20*correctly*%20reporting%20that%20bursting%20is%20not%20needed.%20The%20genes%20it%20does%20flag%20as%20most%20over-dispersed%20tend%20to%20be%20the%20cell-cycle%20markers%20above%3B%20for%20those%20the%20two-state%20band%20is%20modestly%20heavier-tailed%20and%20better-centered%20than%20the%20negative%20binomial%2C%20but%20the%20deep%2C%20clean%20valley%20we%20manufactured%20in%20Phase%201%20is%20absent.%20Allele%20averaging%20and%20capture%20noise%20have%20done%20exactly%20what%20we%20warned%20they%20would%20%E2%80%94%20blunted%20the%20telegraph%20signal%20in%20total%20counts.%0A%0A%20%20%20%20The%20takeaway%20is%20not%20%22the%20model%20failed.%22%20It%20is%20%22the%20data%20do%20not%20contain%20strong%20telegraph%20bimodality%2C%20and%20the%20model%20says%20so.%22%20That%20is%20the%20most%20useful%20thing%20a%20generative%20model%20can%20do%3A%20fit%20the%20bimodal%20genes%20when%20they%20exist%20(Phase%201)%2C%20reduce%20to%20the%20negative%20binomial%20when%20they%20do%20not%20(most%20of%20K562)%2C%20and%20hand%20you%20a%20per-gene%20burst-frequency%20%2F%20regime%20readout%20either%20way.%20Reading%20absolute%20switching%20rates%20off%20a%20monoculture's%20total%20counts%20is%20asking%20more%20than%20droplet%20data%20can%20give%20%E2%80%94%20that%20is%20a%20job%20for%20allele-resolved%20or%20imaging%20assays%2C%20not%20a%20different%20%60scribe%60%20knob.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20A%20closing%20thread%3A%20the%20same%20algebraic%20fault%20line%2C%20twice%0A%0A%20%20%20%20It%20is%20worth%20connecting%20this%20back%20to%20the%20previous%20tutorial's%20*other*%20extension%20%E2%80%94%20the%20hierarchical%20gene-specific%20%24p_g%24%20model%20%E2%80%94%20because%20the%20two%20are%20deeper%20cousins%20than%20they%20look.%0A%0A%20%20%20%20Recall%20what%20made%20the%20negative%20binomial%20Dirichlet%E2%80%93Multinomial%20(NBDM)%20factorization%20so%20clean%3A%20when%20every%20gene%20shares%20the%20same%20success%20probability%20%24p%24%2C%20the%20latent%20rates%20are%20Gammas%20with%20a%20*common*%20rate%20parameter%2C%20and%20independent%20Gammas%20with%20a%20shared%20rate%20normalize%20to%20a%20**Dirichlet**%20on%20the%20simplex.%20That%20single%20algebraic%20fact%20%E2%80%94%20and%20*only*%20that%20fact%20%E2%80%94%20is%20what%20lets%20%60scribe%60's%20differential-expression%20machinery%20operate%20on%20the%20dispersion%20vector%20%24%5Cunderline%7Br%7D%24%20alone%2C%20with%20%24p%24%20canceling%20out%20of%20the%20composition.%0A%0A%20%20%20%20Both%20extensions%20break%20that%20fact%2C%20for%20the%20same%20underlying%20reason%3A%0A%0A%20%20%20%20-%20**Hierarchical%20%24p_g%24**%20(%5B%60_hierarchical_p.qmd%60%5D(https%3A%2F%2Fgithub.com))%3A%20once%20each%20gene%20has%20its%20own%20%24p_g%24%2C%20the%20Gamma%20rates%20%24%5Ctheta_g%20%3D%20p_g%2F(1-p_g)%24%20are%20no%20longer%20shared%2C%20the%20total%20and%20the%20composition%20couple%2C%20and%20the%20closed-form%20Dirichlet%20is%20gone.%20%60scribe%60%20replaces%20it%20with%20a%20Gamma-based%20sampler%20that%20reduces%20*exactly*%20to%20the%20Dirichlet%20when%20all%20%24p_g%24%20coincide.%0A%20%20%20%20-%20**Two-state**%20(%5B%60_two_state_promoter.qmd%60%5D(https%3A%2F%2Fgithub.com))%3A%20the%20latent%20rate%20is%20now%20a%20**Beta**%2C%20and%20the%20Beta%20simply%20does%20not%20have%20the%20Gamma's%20normalize-to-a-Dirichlet%20property.%20The%20composition%20%24%5Crho_g%20%3D%20%5Chat%20r_g%20p_g%20%2F%20%5Csum_j%20%5Chat%20r_j%20p_j%24%20is%20still%20a%20perfectly%20valid%20simplex%20vector%2C%20and%20a%20sample-based%20CLR%20differential-expression%20pipeline%20still%20works%2C%20but%20the%20closed-form%20Dirichlet%20step%20disappears.%0A%0A%20%20%20%20So%20%22gene-specific%20%24p_g%24%22%20and%20%22bursty%20promoters%22%20are%20two%20faces%20of%20the%20same%20move%3A%20*leaving%20the%20special%20Gamma-with-shared-rate%20geometry%20that%20makes%20NBDM%20analytically%20clean.*%20The%20hierarchical%20model%20leaves%20it%20to%20let%20the%20mean%E2%80%93variance%20relationship%20vary%20across%20genes%3B%20the%20two-state%20model%20leaves%20it%20to%20let%20the%20count%20distribution%20change%20*shape*.%20Both%20pay%20the%20same%20price%20(no%20closed-form%20Dirichlet)%20and%20both%20recover%20NBDM%20in%20a%20limit%20(%24%5Csigma_p%20%5Cto%200%24%20for%20the%20hierarchy%2C%20%24k%5E-_g%20%5Cto%20%5Cinfty%24%20for%20the%20two-state%20model).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Conclusion%0A%0A%20%20%20%20This%20tutorial%20extended%20the%20negative%20binomial%20story%20of%20the%20previous%20one%20in%20the%20single%20direction%20the%20negative%20binomial%20cannot%20go%20on%20its%20own%20%E2%80%94%20**distribution%20shape**%20%E2%80%94%20and%20then%20tested%20that%20extension%20honestly%20against%20reality.%0A%0A%20%20%20%20%23%23%23%20What%20we%20did%0A%0A%20%20%20%201.%20**Identified%20the%20pivot.**%20Every%20%60scribe%60%20count%20model%20is%20%60Poisson(latent%20rate)%60%3B%20the%20prior%20on%20that%20rate%20(Gamma%20vs.%20Beta)%20sets%20both%20what%20shapes%20the%20counts%20can%20take%20and%20whether%20the%20Dirichlet%20composition%20survives.%0A%20%20%20%202.%20**Validated%20on%20synthetic%20data%20(Phase%201).**%20On%20data%20generated%20from%20a%20known%20two-state%20process%2C%20the%20negative%20binomial%20could%20not%20fit%20the%20bursty%20genes%20and%20the%20two-state%20model%20recovered%20them%20%E2%80%94%20proving%20the%20capability%20where%20we%20control%20the%20truth.%20We%20then%20mapped%20out%20*identifiability*%3A%20burst%20frequency%20%24k%5E%2B%24%20and%20the%20regime%20%24%5Cdelta%24%20are%20recovered%3B%20the%20OFF%20rate%20%24k%5E-%24%20is%20non-identifiable%20in%20the%20NB%20limit%20(and%20the%20fit%20says%20so%2C%20with%20a%20wide%20posterior)%3B%20the%20absolute%20rate%20%24%5Chat%20r%24%20is%20pinned%20only%20up%20to%20the%20capture%20scale.%0A%20%20%20%203.%20**Applied%20to%20a%20real%20monoculture%20(Phase%202).**%20On%20deeply-sequenced%20K562%20the%20two-state%20model%20reduced%20to%20the%20negative%20binomial%20for%20the%20vast%20majority%20of%20genes%20%E2%80%94%20the%20correct%2C%20honest%20verdict%2C%20since%20allele%20averaging%20and%20capture%20noise%20suppress%20telegraph%20bimodality%20in%20total%20counts.%20The%20genes%20it%20flagged%20as%20most%20over-dispersed%20lined%20up%20with%20the%20cell%20cycle.%0A%0A%20%20%20%20%23%23%23%20When%20to%20reach%20for%20it%0A%0A%20%20%20%20The%20two-state%20model%20earns%20its%20keep%20when%20you%20suspect%20genuine%20**bimodal%20expression%20from%20slow%20promoter%20switching**%20*within%20a%20homogeneous%20population*.%20It%20is%20not%20free%3A%20the%20likelihood%20requires%20a%20one-dimensional%20quadrature%20over%20the%20Beta%20latent%20(Gauss%E2%80%93Legendre%2C%20in%20%60scribe%60%2C%20with%20a%20configurable%20%60n_quad_nodes%60)%2C%20and%20you%20lose%20the%20closed-form%20Dirichlet%20composition.%20For%20genes%20in%20the%20fast-switching%20limit%20it%20simply%20reproduces%20the%20negative%20binomial%20%E2%80%94%20so%20let%20posterior%20predictive%20checks%2C%20not%20the%20loss%20curve%2C%20decide%20whether%20the%20extra%20structure%20paid%20off.%0A%0A%20%20%20%20%23%23%23%20The%20broader%20message%0A%0A%20%20%20%20Three%20threads%20run%20through%20this%20tutorial.%20First%2C%20**validate%20on%20data%20you%20control%20before%20trusting%20real%20data**%20%E2%80%94%20Phase%201%20is%20what%20lets%20us%20read%20Phase%202%20honestly%20instead%20of%20hopefully.%20Second%2C%20**a%20good%20generative%20model%20is%20as%20valuable%20when%20it%20says%20%22no%22%20as%20when%20it%20says%20%22yes%22**%3A%20reducing%20to%20the%20negative%20binomial%20across%20most%20of%20K562%20is%20information%2C%20not%20failure.%20Third%2C%20**the%20limits%20are%20often%20in%20the%20data%2C%20not%20the%20model**%20%E2%80%94%20total%20counts%20from%20a%20diploid%20monoculture%20are%20a%20blunt%20instrument%20for%20telegraph%20kinetics%20(the%20field%20uses%20allele-resolved%20or%20imaging%20assays%20for%20that)%2C%20and%20priors%20on%20weakly-identified%20coordinates%20have%20to%20be%20chosen%20deliberately%2C%20as%20our%20neutral%20regime%20prior%20was.%0A%0A%20%20%20%20%23%23%23%20What%20we%20did%20not%20cover%0A%0A%20%20%20%20-%20**Cell-type%20bimodality**%20%E2%80%94%20when%20the%20two%20modes%20are%20*different%20cell%20types*%20(as%20in%20PBMC)%2C%20the%20right%20tool%20is%20a%20**mixture%20model**%2C%20not%20a%20single%20two-state%20promoter%3B%20a%20single%20Beta%20cannot%20match%20a%20deep%20valley%20*and*%20a%20broad%2C%20over-dispersed%20ON%20population%20at%20once.%0A%20%20%20%20-%20**Differential%20expression%20with%20two-state%20compositions**%20%E2%80%94%20the%20sample-based%20CLR%20pipeline%20that%20replaces%20the%20closed-form%20Dirichlet.%0A%20%20%20%20-%20**Correlated%20two-state%20models%20(TSLN)**%20%E2%80%94%20a%20per-cell%20latent%20layer%20so%20bursty%20genes%20can%20covary%20(and%20the%20ON%20population%20can%20be%20over-dispersed)%2C%20fit%20via%20an%20SVI%E2%86%92Laplace%20cascade%20anchored%20on%20the%20fit%20we%20built%20here.%0A%20%20%20%20%22%22%22)%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
db41218960b1d595551267b3b0b7c11e