lineage.states.StateDistributionCR ================================== .. py:module:: lineage.states.StateDistributionCR .. autoapi-nested-parse:: Competing-risks state distributions. The Gamma/GaPhs emissions in this package treat a cell's fate (Bernoulli) and its phase duration (Gamma) as independent observations, and discard the duration of any cell that dies. That throws away every death time and, more subtly, mis-states the likelihood of a time-censored cell: "no event yet at time t" is written as P(division > t) when it should be P(division > t AND death > t). Here each phase instead carries two latent clocks, T_D ~ Gamma(a, s) division / transition T_X ~ Gamma(a_x, s_x) death and we observe min(T_D, T_X) together with an indicator of which fired. The three likelihood cases are the standard competing-risks ones: transition seen at t f_D(t) * S_X(t) death seen at t f_X(t) * S_D(t) censored at t S_D(t) * S_X(t) The division probability is then *derived*, P(divide) = int f_D(t) S_X(t) dt, rather than fit as a free Bernoulli parameter, so the death fraction and the death timing are forced to agree. With the G1 death clock pinned to a constant hazard (shape 1) this costs no degrees of freedom relative to the Bernoulli/Gamma model. Attributes ---------- .. autoapisummary:: lineage.states.StateDistributionCR.TIME_FLOOR Classes ------- .. autoapisummary:: lineage.states.StateDistributionCR.StateDistribution lineage.states.StateDistributionCR.StateDistributionPhase Functions --------- .. autoapisummary:: lineage.states.StateDistributionCR._gamma_logpdf lineage.states.StateDistributionCR._gamma_logsf lineage.states.StateDistributionCR.event_masks lineage.states.StateDistributionCR.exponential_estimator lineage.states.StateDistributionCR.fit_clocks lineage.states.StateDistributionCR.atonce_estimator Module Contents --------------- .. py:data:: TIME_FLOOR :value: 1e-10 .. py:function:: _gamma_logpdf(t, a, scale) log Gamma(a, scale) pdf, vectorized directly against scipy.special. .. py:function:: _gamma_logsf(t, a, scale) log Gamma(a, scale) survival function, vectorized directly against scipy.special. .. py:function:: event_masks(x) Split cells into the three competing-risks cases. ``x`` has the usual three per-phase columns ``[fate, duration, censoring]``, where fate is 1 for surviving the phase, 0 for dying in it, and NaN when unknown, and the censoring flag is 1 when the phase was seen through to its end. Cells whose duration is negative have been masked for cross validation, and cells with a NaN duration never entered the phase; both are excluded everywhere. :return: boolean masks for (division observed, death observed, censored) .. py:function:: exponential_estimator(obs, events, weights, param_idx, K) Weighted right-censored MLE for exponential scales, one per group. With a constant hazard the MLE is total time at risk over events observed, which needs no iteration. A pseudocount keeps a group with no observed deaths finite. .. py:function:: fit_clocks(distributions, x_list, gammas_list, state_j) Fit both clocks of one state, sharing each shape across conditions. ``distributions`` is either a single :class:`StateDistribution` or a list of them, one per condition. The shape parameters are shared across conditions and the scales are free, mirroring how ``atonce_estimator`` treats the Gamma model. .. py:function:: atonce_estimator(all_tHMMobj, x_list, gammas_list, phase) M step across several conditions at once, matching the Gamma model's interface. .. py:class:: StateDistribution(gamma_a = 7.0, gamma_scale = 4.5, death_a = 1.0, death_scale = 40.0, fixed_death_shape = True) One cell-cycle phase with competing division and death clocks. ``params`` is ``[bern_p, gamma_a, gamma_scale, death_a, death_scale]``. The first three entries keep the meaning they have in :class:`~lineage.states.StateDistributionGamma.StateDistribution`, so downstream figure code that indexes them positionally continues to work; ``bern_p`` is now derived from the two clocks rather than fit. .. py:attribute:: atonce_estimator .. py:attribute:: fixed_death_shape :value: True .. py:attribute:: params .. py:property:: div_clock .. py:property:: death_clock .. py:method:: division_probability() P(the division clock fires first) = int f_D(t) S_X(t) dt. ``quad`` evaluates the integrand at ~100-200 points per call, and this runs once per state per M step. Using the frozen ``rv_continuous`` machinery for each point (as ``div.pdf`` / ``death.sf`` do) costs orders of magnitude more than the arithmetic itself, so the integrand is written directly against ``scipy.special`` instead. .. py:method:: rvs(size, rng=None) Draw min(T_D, T_X) and record which clock fired. .. py:method:: dist(other) Wasserstein distance between the division clocks of two states. Kept on the division clock alone so that the number is comparable with the Bernoulli/Gamma model's. .. py:method:: dof() Two division-clock parameters plus the death clock's scale, and its shape when that is not pinned. The Bernoulli is derived, so it is not counted. .. py:method:: logpdf(x) Competing-risks log likelihood of each cell's phase observation. Called once per state on every E step, over the whole lineage array, so it is vectorized directly against ``scipy.special`` rather than going through the frozen ``rv_continuous`` machinery (``self.div_clock.logpdf`` etc.), which pays generic validation/broadcasting overhead on every call regardless of array size. .. py:method:: estimator(x, gammas) Weighted M step for a single condition. The two clocks separate in the complete-data likelihood, so each is just a weighted right-censored Gamma fit over every timed cell. .. py:method:: censor_lineage_array(censor_condition, tree, obs, states, desired_experiment_time=2000000000000.0) Applies censoring to array representation directly. .. py:class:: StateDistributionPhase(gamma_a1 = 7.0, gamma_scale1 = 3.0, gamma_a2 = 14.0, gamma_scale2 = 6.0, death_scale1 = 40.0, death_a2 = 3.0, death_scale2 = 20.0) G1 and G2 phases, each with its own pair of competing clocks. ``params`` is ``[bern_p1, bern_p2, a1, s1, a2, s2, death_a1, death_s1, death_a2, death_s2]``. The leading six entries match :class:`~lineage.states.StateDistributionGaPhs.StateDistribution` exactly. .. py:attribute:: atonce_estimator .. py:attribute:: G1 .. py:attribute:: G2 .. py:attribute:: params .. py:method:: _sync() Mirror the sub-distributions' parameters into the flat ``params`` array. .. py:method:: rvs(size, rng=None) .. py:method:: dist(other) .. py:method:: dof() .. py:method:: logpdf(x) .. py:method:: estimator(x, gammas) .. py:method:: censor_lineage_array(censor_condition, tree, obs, states, desired_experiment_time=2000000000000.0)