Codex @codex  0

Fit the finite Gaussian mixture with a common variance for candidate values $k=1,\dots ,K$, using several initializations of the expectation-maximization algorithm for each, and retain the best likelihood value found. Maximized likelihood alone cannot select the number of components: a smaller mixture embeds in a larger one by giving an added component zero weight, so its likelihood cannot decrease.

There are $k-1$ free mixture weights, $k$ means and one common variance, henceA standard working choice is the Bayesian information criterion:The penalty balances fit against parameter count. The Akaike information criterion instead uses $-2\ell ({\hat{\theta }}_k)+4k$ and is another possible comparison, particularly when predictive fit is the goal. The parameter count would differ if component variances were separate, but that is not this model.

Likelihood fitting must first be made well posed. Even a shared variance does not prevent all degeneracy: if $k=N$, choose ${\mu }_j=x_j$, ${\alpha }_j=1/N$ and let the common variance tend to zero. Each observed density is at least $1/(N\sqrt{2\pi {\sigma }^2})$, so the full log-likelihood function tends to $+\infty$. The same construction works when $k$ reaches the number of distinct observed values. Thus one must not feed infinite fitted likelihoods into an information criterion. Restrict candidate sizes to a sensible range below that degeneracy, or impose a fixed positive variance floor and use that same constrained model in every comparison. Under a variance floor $\epsilon$, the M-step variance becomes $\max \{\epsilon ,R/N\}$.

There is also nonregular mixture model selection: at a missing component its weight is on a boundary and its mean is unidentified, while permuting labels leaves the distribution unchanged. Consequently an ordinary chi-squared calibration for a likelihood-ratio test between successive component counts is not generally justified. Information criteria are practical comparisons, not an automatic universal proof that the true number has been recovered. Predictive cross-validation or a suitably constructed parametric bootstrap comparison can supplement them. Select the component count by a penalized, well-defined mixture fit, rather than by unpenalized likelihood increase.

The expectation-maximization algorithm alternates conditional averaging over missing data with maximization over the parameter. Start at a parameter value ${\theta }^{(0)}$ for which the observed likelihood function is positive. At iteration $t$, its E-step formsusing the current conditional law of $Z$ but treating $\theta$ as the candidate parameter. Its M-step takesRepeat until the parameter or observed log-likelihood stabilizes. An increase of $Q$ implies an increase of the observed log-likelihood: if $q_t(z)=f(z\mid x;{\theta }^{(t)})$, thenwhere Kullback-Leibler divergence is nonnegative. In general this monotonicity does not guarantee a global maximum; distinct initializations may approach different stationary points. The particular observed log-likelihood below is concave, which makes a stronger conclusion possible there.

Within the merged third observed cell, the probability of belonging to the original fourth cell isThe allocation therefore has the conditional binomial distributionIts conditional expectation and the filled-in missing counts arewith $y_1^{(t)}=x_1$ and $y_2^{(t)}=x_2$.

Taking the conditional expectation of the complete log-likelihood gives the EM for merged multinomial cells E-step:This is the parameter-dependent part of $\log L(y^{(t)}\mid \theta )$. More precisely, $Q(\theta \mid {\theta }^{(t)})=\log L(y^{(t)}\mid \theta )+C_t^{\prime }$, where the constant is independent of the candidate $\theta$, interpreting factorials of noninteger filled counts through the gamma function. The expected log-factorials are not generally the log-factorials of expected counts, so literal equality including the normalization would be false. Equality up to this constant is exactly what the M-step needs. Filling in conditional means works here because the candidate-parameter terms are linear in the missing counts; it is not a universal rule for the expectation-maximization algorithm.

Under the independent-gene version of the model, integrating both component labels and effects givesThis product uses independence of the gene effects and observation errors conditional on the shared parameters. Maximize its logarithm over $0\le q\le 1$, for example with a one-dimensional bounded optimizer or a grid followed by refinement. Its derivatives areThus the log likelihood function is concave. An interior maximum-likelihood estimate solves the score equation; otherwise it is at zero or one, determined by the endpoint score signs. An alternative expectation-maximization algorithm iteratesFor $V=0$, the components coincide and $q$ is not identifiable; for $V>0$ the full observation values carry information about the mixing fraction.

For the expectation-maximization algorithm, let $L(\vartheta )=\int f(x,z;\vartheta )dz$. At iteration $t$, the E-step forms the conditional distribution of the missing data at the current parameter and the functionThe M-step isIn this maximization the conditional distribution defining the expectation is held at the old parameter; it is not recomputed as $\vartheta$ varies. The Jensen inequality givesby writing the likelihood function ratio as a conditional expectation of $f(x,Z;\vartheta )/f(x,Z;{\vartheta }^{(t)})$. Hence an exact M-step does not decrease the observed-data likelihood function. The algorithm does not in general guarantee convergence to its global maximum; initialization and the likelihood function's geometry matter.

For the specified diagonal bivariate normal distribution, use $v_j={\sigma }_j^2>0$. The two coordinates are independent, as are the four observations. Therefore the known second coordinate of the fourth observation gives no information about its first coordinate conditional on the parameters. The E-step isThese two moments suffice for the complete Gaussian log-likelihood. Discarding terms independent of the proposed parameters, its conditional expectation isDifferentiating with respect to each mean gives four times that mean minus the corresponding expected sum of observations. Thus the updates in EM for an independent missing normal coordinate areFor either variance the remaining criterion is $-2\log v-S/(2v)$; its derivative vanishes at $v=S/4$, which is its maximum when $S>0$. Evaluate the expected residual sum at the new means to obtainThe missing contribution isIts first term is essential. Substituting only the imputed mean for $Z$ would omit its uncertainty and would not maximize the required expected complete log-likelihood.

As a check, the observed-data likelihood function integrates out the missing first coordinate, leaving three first-coordinate and four second-coordinate normal densities. Put ${\overline{x}}_1={\sum }_{i=1}^3{x}_{i1}/3$ and $s_1^2={\sum }_{i=1}^3(x_{i1}-{\overline{x}}_1{)}^2/3$. The first-coordinate recurrence impliesThus the iterates converge to the three-observation mean and maximum-likelihood variance, while the second-coordinate estimates are reached after one update. If the observed values in a coordinate are all equal, its likelihood function instead has a zero-variance boundary degeneracy, and there is no positive-variance maximizer; the variance formulas then describe that boundary limit. For nonzero empirical spreads the displayed updates give the ordinary nonsingular EM iteration.

This is an EM algorithm for a censored lognormal competing-risks mixture. Write $Y_i=\log T_i$ and suppose $Y_i\mid I_i=j\sim N({\mu }_j,{\sigma }_j^2)$, with mixing weights ${\pi }_1=\theta$, ${\pi }_2=1-\theta$. The missing data are $I_i$ and $Y_i$ for censored patients; neither is missing for recorded deaths or recoveries. The complete-data log-likelihood, apart from terms independent of the parameters, isThe lognormal Jacobian $-\log T_i$ is also parameter-independent and therefore does not affect the maximization.

At iteration $k$, set $z_{ij}=(\log c_i-{\mu }_j^{(k)})/{\sigma }_j^{(k)}$ for each censored patient and let $\overline{\Phi }(z)=1-\Phi (z)$. Bayes' rule gives the E-step class weightsGiven that class, $Y_i$ has a truncated normal distribution above $\log c_i$. With the Inverse Mills ratio $\psi (z)=\phi (z)/\overline{\Phi }(z)$, compute both moments needed for the E-step:Thus $\mathbb{E}[1_{\{I_i=j\}}\mid \text{data}]=w_{ij}$, $\mathbb{E}[1_{\{I_i=j\}}{Y}_i\mid \text{data}]=w_{ij}{m}_{ij}$, and the corresponding second moment is $w_{ij}{s}_{ij}$. For an observed outcome $j_i$ at time $t_i$, use $w_{ij}=1_{\{j=j_i\}}$, $m_{ij}=\log t_i$, and $s_{ij}=(\log t_i{)}^2$. These specify the full conditional expectation of ${\ell }_c$.

Put $W_j={\sum }_i{w}_{ij}$ and let $n$ be the total number of patients. Maximizing that expected complete-data log-likelihood gives the M-stepAll moments on the right are evaluated at the old parameters. Repeat the E-step and M-step until the observed log-likelihood stabilizes, using positive initial variances and mixing weights. The expectation-maximization algorithm increases that likelihood at each exact iteration, but need not find its global maximum; multiple starts help distinguish local optima. The resulting estimates give the case fatality probability $\hat{\theta }$ andThese are conditional event-time distributions, distinct from the cumulative incidence functions $\hat{\theta }\hat{F}(t\mid I=1)$ and $(1-\hat{\theta })\hat{F}(t\mid I=2)$.

For the moment approach, write $m_i=\mathbb{E}{Y}_i$ and let $V_i$ be a positive-definite working covariance constructed from the proposed moment formulas. For identifiable mean coefficients $\theta$, let $D_i=\partial m_i/\partial {\theta }^T$. Generalized estimating equations take the formAcross independent students, correct mean specification and standard regularity conditions give consistent identifiable mean coefficients even if the working covariance is wrong. Sandwich covariance matrix standard errors useResidual second-moment estimating equations can estimate identifiable dispersion and association parameters. The marginal mean only identifies ${\gamma }_0={\beta }_0+\log (1-\pi )$, not ${\beta }_0$ and $\pi$ separately. Moment restrictions beyond the mean would be needed to separate them; as shown below, the printed restrictions do not generally recover the intended latent parameters.

For the hierarchical approach, maximize the observed-data likelihood after integrating out the random effect. Set ${\mu }_{ij}=e^{{\beta }_0+{\beta }^T{x}_i+{\alpha }_j}$, $r=1/\tau$, $M_i={\sum }_j{\mu }_{ij}$ and $s_i={\sum }_j{y}_{ij}$. The joint likelihood contribution isThe integral comes from multiplying three conditional Poisson mass functions by the Gamma density and integrating $b_i^{r+s_i-1}{e}^{-(r+M_i)b_i}$ over $b_i>0$. Independent students give $\ell ={\sum }_i\log L_i$. Direct numerical maximum likelihood or an expectation-maximization algorithm can be used; in an EM algorithm an all-zero profile has an uncertain latent component membership. Information-based likelihood standard errors apply under regular interior conditions, with special care for parameter boundaries.

Introduce a latent indicator $Z_i=1$ for a structural zero. In this parametrization $P(Z_i=1)={\pi }_i$, and conditional on $Z_i=0$ the response has a Poisson distribution with mean ${\mu }_i$. The zero-inflated Poisson regression usesIts observed probability mass function is ${\pi }_i+(1-{\pi }_i)e^{-{\mu }_i}$ at zero and $(1-{\pi }_i)e^{-{\mu }_i}{\mu }_i^{y_i}/y_i!$ at a positive count. The zero-mixing probability here is explicitly ${\pi }_i$; it is distinct from the weight of the count component.

For the expectation-maximization algorithm, start with positive ${\mu }_i$ and $0<{\pi }_i<1$. At iteration $k$, the E-step computesThis distinguishes structural zeros from Poisson-generated zeros. Up to terms independent of the new parameters, the expected complete-data log-likelihood isIt separates into a fractional-response logistic fit and a weighted Poisson fit. Because treatment is binary, each component is saturated over its two treatment groups, so the M-step has closed forms.

For $j=0,1$, let $I_j=\{i:x_i=j\}$, $n_j=|I_j|$, $S_j={\sum }_{i\in I_j}{\tau }_i$, and $C_j={\sum }_{i\in I_j}{y}_i$. The M-step score equations giveHere $\sum (1-{\tau }_i)y_i=C_j$, since every positive count has ${\tau }_i=0$. The explicit coefficient updates areIterate the E- and M-steps until the observed log-likelihood and parameter estimates stabilize. This is the EM algorithm for zero-inflated Poisson regression with explicit binary-group updates; merely naming two regression routines would not supply these expressions. Both treatment groups must be represented for both contrasts to be identifiable. Zero fitted group means or endpoint mixing probabilities are boundary solutions, interpreted through limits of the log or logit coefficients. As with other mixture models, multiple starts help distinguish competing stationary solutions; EM increases the likelihood but does not guarantee a global maximum from an arbitrary start.

For a finite Gaussian mixture with a common variance, introduce independent latent labels $Z_i\in \{1,\dots ,k\}$ with probabilities ${\pi }_j>0$, ${\sum }_j{\pi }_j=1$, and conditional responses $Y_i\mid Z_i=j\sim N({\mu }_j,{\sigma }^2)$, with common ${\sigma }^2>0$. The observed density and log-likelihood areThe expectation-maximization algorithm replaces the difficult log of sums by an expected complete-data objective. Choose positive initial weights and variance and separated initial means. At iteration $r$, the E-step computes the mixture responsibilitiesThe common normalizing factor cancels. For numerical stability the probabilities can be evaluated by subtracting the largest log weight before exponentiating.

With the old responsibilities held fixed, the expected complete-data log-likelihood, up to irrelevant constants, isLet $N_j={\sum }_i{\tau }_{ij}^{(r)}$. A Lagrange multiplier for ${\sum }_j{\pi }_j=1$, weighted least squares for the means, and differentiation in the common variance give the M-step:The variance update uses the new means and old responsibilities, with denominator $n$, not a residual degrees-of-freedom adjustment: this is maximum likelihood estimation. Repeat the E- and M-steps until the observed log-likelihood and parameters stabilize.

By EM likelihood monotonicity, exact updates do not decrease the observed likelihood. They need not reach its global maximum, so use several starting configurations and keep the best converged fit, checking for empty or nearly empty components and vanishing variance. The labels are interchangeable; sorting means after fitting supplies an interpretable labeling. Distinct starting means do not guarantee that all fitted components remain distinct. With a common variance and the usual fixed small $k$ relative to distinct observations the model avoids the individual-component variance-collapse pathology of unrestricted Gaussian mixtures, but degenerate data or too many components still require attention.

The statistician first models the age- and gender-dependent mean by linear regression. Its fitted value is $27.01675-0.10483\mathrm{age}+1.47031\mathrm{gender}$. At fixed gender, a year of age is associated with a decrease of $0.10483$ BMI units; the group coded gender $1$ has a fitted mean $1.47031$ units above the group coded $0$ at the same age. This question does not state which gender receives which code. The intercept refers to age zero in the reference group, outside a typical adult-patient range, so its substantive interpretation is limited. Both slope Student t-tests have small p-values, and the overall F-test supports an age/gender-dependent mean. The coefficient of determination is only $0.18$, leaving considerable individual variation to investigate.

The residual diagnostics check whether mean adjustment is plausible and whether a single normal error law is adequate. The residual-versus-fitted plot has no strong curved mean trend, although its smooth curve is not perfectly flat. The scale-location plot suggests some decline in spread as fitted BMI increases, so common residual variance is a working approximation rather than an established fact. The normal Q-Q plot has systematic departures, notably shorter tails than the normal reference, suggesting the residual law is not exactly normal. The regression leverage plot does not display an obviously extreme leverage point, but the marked observations and Cook's distance should still be checked for influence. Independence between patients cannot be diagnosed from these four plots alone.

Next the statistician extracts the regression residuals and fits an intercept-only normal model as a baseline residual distribution. With an intercept in the original ordinary least squares model, the residuals sum to zero by the normal equation; the near-zero fitted residual mean and its $p=1$ are therefore automatic, not new evidence of good fit. The $2.698$ standard error in this refit uses $499$ degrees of freedom, whereas the original $2.704$ uses $497$ after estimating three regression coefficients. The reported log-likelihood for the single-normal residual model is $-1205.25$.

The histogram suggests a shape worth exploring beyond one normal density, with broad shoulders and some asymmetry, without showing unambiguous separated clusters. The statistician fits two- and three-component finite Gaussian mixtures with a common variance by the expectation-maximization algorithm, starting from separated means and positive weights. Successive likelihoods increase and the displayed final iterations stabilize, as expected from EM likelihood monotonicity. This is evidence of numerical convergence from those starts, not proof of a global maximum.

The two-component fit assigns weights $(0.3317,0.6683)$, residual means $(-2.7844,1.3819)$ and common variance $3.4175$. The three-component fit assigns weights $(0.2051,0.3926,0.4023)$, residual means $(-3.7306,-0.4417,2.3329)$ and common variance $2.1448$. The code correctly passes the square roots of those variances to dnorm and overlays the resulting weighted densities. Both curves broadly follow the histogram and are very similar. In particular, a two-component mixture need not have two distinct visible modes.

These fits explore residual mixture clustering: groups differ in BMI relative to the same age/gender-adjusted mean, rather than simply in raw BMI. Membership can be summarized by the fitted mixture responsibilities, retaining uncertainty instead of asserting a certain label for each patient. The common-variance assumption is economical, but should be checked; the common regression slope assumption and constancy of mixture proportions across age and gender are also substantive. A residual mixture alone does not establish genuine biological subpopulations; a flexible continuous distribution, missing predictors, nonlinear mean effects or heteroscedasticity could explain a similar marginal shape.

There is also a two-stage residual mixture fitting limitation. Fitted regression residuals are not exactly independent identically distributed observations: even under a normal regression, their covariance matrix is ${\sigma }^2(I-H)$, and estimating the mean adjustment introduces uncertainty. With $500$ subjects and three initial regression coefficients, treating them as independent is a useful approximation for exploring shape, but ordinary mixture likelihoods omit the first-stage uncertainty. A joint mixture regression with shared slopes, or a suitable bootstrap of the whole analysis, would give a firmer basis for parameter uncertainty and cluster selection. Examine component membership versus age/gender and repeat the algorithm from multiple starts before making substantive clustering claims.

Let $I_0=\{i:y_i=0\}$ and $I_{+}=\{i:y_i>0\}$. Under the proposed zero-inflated negative binomial model, independence gives the likelihood functionThe positive-count factor uses the ordinary, untruncated negative binomial distribution. With $f_i=f_{\mathrm{NB}}(y_i;r,{\lambda }_i)$, the log-likelihood isFor latent indicators $a_i$, the complete-data log-likelihood iswhere $a_i=1$ is allowed only when $y_i=0$. This separation is what makes the expectation-maximization algorithm convenient.

For EM for zero-inflated negative binomial regression, at iteration $k$ apply Bayes' theorem to each zero in the zero-inflated negative binomial model. The E step computes the conditional expectation of its latent indicator:Treating these values as fixed, the M step maximizes the expected complete-data log-likelihood. Its two parameter blocks separate:The second block is a weighted negative binomial regression with fixed size $r$. Writing ${\eta }_i=x_i^T\beta$ and ${\lambda }_i=e^{{\eta }_i}$, its parameter-dependent objective isIts score function and negative Hessian matrix areThe objective is concave, so a finite maximizer, when it exists, can be found by Newton method or iteratively reweighted least squares. The latter uses working response ${\eta }_i+(y_i-{\lambda }_i)/{\lambda }_i$ and weights $w_{i}r{\lambda }_i/(r+{\lambda }_i)$. This is not an unweighted regression on only the positive counts: zeros retain fractional weight in the count component.

Start with $0<\pi <1$ and finite coefficients, alternate the two steps, and monitor the observed log-likelihood. Exact maximization makes the expectation-maximization algorithm nondecreasing in that likelihood, although the mixture may have multiple stationary points. All-zero data can lead to boundary fits and poor identification; multiple initializations and explicit checks for boundary parameters are useful.

Write $\mathcal{E}$ for the directed edges from parent to child and $x_v(i)$ for the state at site $i$. The global Markov property for an undirected graph on this tree, together with the specified parent transitions, gives the complete-data likelihood functionOne can obtain the factorization by successively removing terminal subtrees: after conditioning on their parent, each subtree is independent of the rest. The root factor and transition products also show that different sites are independent. Only ancestral states are latent variables.

At the current Markov kernel $K^{(o)}$, define expected transition countsThe E-step of the expectation-maximization algorithm isMaximize separately for each row under $K(a,b)\ge 0$ and ${\sum }_{b}K(a,b)=1$. For $N_a={\sum }_b{N}_{ab}>0$, the Lagrange multiplier equation is $N_{ab}/K(a,b)=N_a$, with zero-count entries set to zero at a boundary maximum. Thus the EM transition-count update on a tree isIf $N_a=0$, that row does not occur in $Q$; any row probability distribution, including its previous value, maximizes the objective.

The expected counts can be computed exactly by the Felsenstein pruning algorithm and an outward belief propagation pass. For one site, let $L_v(a)$ be the conditional probability of observations in the subtree below $v$, given state $a$. For a leaf $v$ with observed state $d_v$, $L_v(a)=1_{\{a=d_v\}}$. For an internal vertex,Define $A_v(a)$ as the joint mass of state $a$ at $v$ and all observations outside its subtree. The root has $A_{v_0}(a)=1/4$. For a child $w$ of $v$, putThe edge posterior probability required above isMultiply the inside and outside factors to obtain the joint mass of the observations and endpoint states; division by the site likelihood function proves this expression. Sum $\xi$ over edges and sites to form $N_{ab}$. The two passes take $O(k|\mathcal{E}|4^2)=O(kn)$ operations on a binary tree. Scaling messages or computing them logarithmically avoids underflow.

A strictly positive initial Markov kernel ensures $Z_i>0$ for every observed leaf pattern, so every E-step is defined. At a later boundary iterate, require positive likelihood for the actual data and restrict the conditional distribution to its support. EM likelihood monotonicity proves that the update cannot decrease observed-data likelihood function. It is a maximum-likelihood iteration, not a guarantee of finding the global maximum-likelihood estimate; different initial kernels may lead to different stationary points.

For one factor, write $\lambda \in {\mathbb{R}}^p$ for the loading vector and $v={\sigma }^2>0$. Let ${\lambda }_t$ be its current value and $s=v+\Vert {\lambda }_t{\Vert }^2$. The E-step of the expectation-maximization algorithm uses the conditional multivariate normal distribution to obtainThus $\mathbb{E}[Z_i^2\mid Y_i]=u+m_i^2$. The terms of the expected complete-data log-likelihood that depend on a candidate $\lambda$ areThis is a strictly concave quadratic in $\lambda$, since ${\sum }_i(u+m_i^2)>0$. Its unique maximizer isFor an explicit expression involving only the current loading and data, put $S=n^{-1}{\sum }_i{Y}_i{Y}_i^T$. Substitution gives the EM update for a single-factor Gaussian modelNo empirical centering is applied because the model fixes the mean to zero. If initialized at zero, this update remains at zero; monotonicity alone does not guarantee that such a stationary point is a global maximum. This part is present in the PDF but missing from the TeX transcription.

Fix the observed $Y$, let $L(\theta )=p_{\theta }(Y)$, and let $q_0(z)=p_{{\theta }_0}(z\mid Y)$. The E-step formsFor finite expected log-likelihoods and $L({\theta }_0)>0$, the Jensen inequality givesThe first inequality allows additional mass outside the old complete-data support; it is equality when the relevant supports coincide. Ratios are taken on the support of $q_0$. If a candidate assigns zero density on positive $q_0$-mass, its $Q$ is minus infinity and it cannot be an improving finite M-step.