Worked Nonlinear Gaussian Likelihood Estimation

Why Work Through a Nonlinear Gaussian Fit?

A nonlinear regression can describe a response whose mean changes in a curved or saturating way, while allowing its variability to change across observations. Maximum likelihood estimation handles both features in one model, but the likelihood, score, and uncertainty calculation must all reflect the variance specification. This tutorial develops those calculations for a concrete Gaussian model and shows how to interpret a contrast between fitted means.

Suppose observations \((y_i,x_i,z_i)\), \(i=1,\ldots,n\), are conditionally independent, with mean \(m_i(\beta)\) and variance \(v_i(\gamma)>0\). The parameter vector \(\theta\) contains the mean parameters \(\beta\) and variance parameters \(\gamma\). For illustration, use the nonlinear mean \(m_i=a\exp(bx_i)\) and the variance \(v_i=\exp(\eta+\lambda z_i)\). The logarithmic variance specification guarantees positive variances for every real \(\eta,\lambda\).

Model and assumptions: Conditional on the covariates, the responses are independent and \(Y_i\sim N(m_i(\beta),v_i(\gamma))\). The mean and variance functions are correctly specified, the variances are positive, and the true parameter is an interior point of the parameter space. For the local covariance approximation, the likelihood is assumed sufficiently regular and identifiable near the solution.

The Gaussian assumption describes the full conditional distribution, not only the first two moments. In particular, each observation contributes a normal density with its own mean and variance. If the variance were constant, the variance function would simplify; here it is part of the parameter estimation problem.

Likelihood, Log Likelihood, and Score

For one observation, the conditional density is \(f(y_i\mid x_i,z_i;\theta)=(2\pi v_i)^{-1/2}\exp\{-(y_i-m_i)^2/(2v_i)\}\). Independence makes the sample likelihood the product of these densities. Taking logarithms gives the expression used for estimation:

\ell(\theta) =-\frac{1}{2}\sum_{i=1}^n \left\{\log(2\pi)+\log v_i(\gamma) +\frac{[y_i-m_i(\beta)]^2}{v_i(\gamma)}\right\}.

The maximum-likelihood estimator (MLE) is a parameter value that maximizes \(\ell(\theta)\). Define the residual \(r_i=y_i-m_i(\beta)\), the mean derivative \(d_i=\partial m_i/\partial\beta\), and the variance derivative \(q_i=\partial\log v_i/\partial\gamma\), each understood as a vector when there are multiple parameters. Differentiating the log likelihood gives separate score contributions:

U_\beta(\theta)=\frac{\partial\ell}{\partial\beta} =\sum_{i=1}^n\frac{r_i}{v_i}d_i, \qquad U_\gamma(\theta)=\frac{\partial\ell}{\partial\gamma} =\frac{1}{2}\sum_{i=1}^n \left(\frac{r_i^2}{v_i}-1\right)q_i.

The mean score weights each mean derivative by the residual divided by its variance. The variance score compares each squared standardized residual \(r_i^2/v_i\) with its model expectation of one. At an interior maximum, the full score is zero. In a nonlinear model these equations generally have no closed-form solution, so an optimizer searches for a solution numerically.

Worked Example: Evaluating a Gaussian Log Likelihood

Consider three hypothetical observations at \(x=(0,1,2)\), with responses \(y=(1,2,4)\). Evaluate the log likelihood at the candidate mean parameters \(a=2,b=0\) and variance parameters \(\eta=0,\lambda=0\), using \(z_i=x_i\). This is a likelihood evaluation at a candidate parameter value, not a claim that the value is the MLE.

The mean is \(m_i=2\exp(0x_i)=2\) for each observation, and the variance is \(v_i=\exp(0+0x_i)=1\). Thus the residuals are \((-1,0,2)\), and their squared values sum to \(1+0+4=5\). Substitution into the log likelihood gives

\ell =-\frac{3}{2}\log(2\pi)-\frac{3}{2}\log(1)-\frac{1}{2}(5) =-\frac{3}{2}\log(2\pi)-2.5 \approx -5.2568.

The \(\log(1)\) term is zero. The numerical value follows from \(\log(2\pi)\approx1.83788\), so the first term is approximately \(-2.75682\), and adding \(-2.5\) gives approximately \(-5.25682\). This calculation illustrates why the normalizing constant and the residual term both belong in the likelihood.

Fisher Scoring as a Numerical Update

One way to find an MLE is to use Fisher scoring. At an iterate \(\theta^{(t)}\), it updates the parameter by multiplying the score by the inverse expected information:

\theta^{(t+1)} =\theta^{(t)}+ \mathcal{I}(\theta^{(t)})^{-1}U(\theta^{(t)}), \qquad \mathcal{I}(\theta) =E_\theta\!\left[ U(\theta)U(\theta)^\mathsf{T} \right].

For the independent Gaussian model, the expected information blocks for mean and log-variance parameters are \(\sum_i d_i d_i^\mathsf{T}/v_i\) and \(\frac12\sum_i q_i q_i^\mathsf{T}\), respectively; the expected cross-block is zero. This follows because a centered normal residual has zero third moment, so the mean and variance score contributions have zero expected product. The blocks are evaluated at the current iterate. This is a local update, not a guarantee that one step reaches the maximum.

Worked Example: One Fisher-Scoring Step

Use the model \(m_i=a\exp(bx_i)\), \(v_i=\exp(\eta+\lambda x_i)\) with hypothetical data \(x=(0,1,2)\) and \(y=(2,3,1)\). Start at \(a=2,b=0,\eta=0,\lambda=0\). At this point all fitted means and variances are 2 and 1, respectively, and the residuals are \((0,1,-1)\). The mean derivatives with respect to \(a,b\) are \((1,1,1)\) and \((0,2,4)\).

The mean score is therefore \(U_{(a,b)}=(0+1-1,\;0+2-4)^\mathsf{T}=(0,-2)^\mathsf{T}\). The expected information for these parameters is

\mathcal{I}_{(a,b)} = \begin{pmatrix}1&1&1\\0&2&4\end{pmatrix} \begin{pmatrix}1&0\\1&2\\1&4\end{pmatrix} = \begin{pmatrix}3&6\\6&20\end{pmatrix}.

Its inverse is \(\frac{1}{24}\begin{pmatrix}20&-6\\-6&3\end{pmatrix}\), so the mean-parameter update is \((0.5,-0.25)^\mathsf{T}\). The variance score uses \(q_i=(1,x_i)^\mathsf{T}\). Since \(r_i^2/v_i=(0,1,1)\), it is \(U_{(\eta,\lambda)}=\frac12(-1+0+0,\;0+0+0)^\mathsf{T}=(-0.5,0)^\mathsf{T}\).

The variance expected information is \(\frac12\sum_i q_iq_i^\mathsf{T} =\begin{pmatrix}1.5&1.5\\1.5&2.5\end{pmatrix}\). Its inverse is \(\frac{1}{1.5}\begin{pmatrix}2.5&-1.5\\-1.5&1.5\end{pmatrix}\), giving the update \((-0.8333,0.5)^\mathsf{T}\). Thus this single Fisher-scoring step gives \((a,b,\eta,\lambda)=(2.5,-0.25,-0.8333,0.5)\). It is only an intermediate iterate; further iterations and a convergence check are needed before treating a value as a fitted MLE.

Covariance and Uncertainty for Fitted Means

After obtaining an interior MLE \(\hat\theta\), a standard local covariance estimate under the correctly specified likelihood is the inverse observed information, where the observed information is the negative Hessian of the log likelihood:

\widehat{\operatorname{Cov}}(\hat\theta) = \left[ -\frac{\partial^2\ell(\theta)} {\partial\theta\,\partial\theta^\mathsf{T}} \bigg|_{\theta=\hat\theta} \right]^{-1}.

Under the usual regularity conditions, this matrix estimates the covariance of the MLE in a large sample. The inverse expected information evaluated at the MLE is another common large-sample approximation. These are local approximations: they describe uncertainty near the fitted parameter and do not remove the need to check model specification. A fitted mean is a function of estimated parameters, so its uncertainty must account for their covariance.

For a scalar contrast \(c(\theta)\), the delta method uses its gradient \(g=\partial c(\theta)/\partial\theta\) at \(\hat\theta\). If \(\widehat V\) is the estimated parameter covariance matrix, then the estimated variance of the contrast is \(g^\mathsf{T}\widehat Vg\), and its standard error is the square root of that quantity. This propagates uncertainty in the parameters to uncertainty in the contrast.

Worked Example: Contrasting Two Fitted Values

Suppose a hypothetical fitted model has mean \(m(x;a,b)=a\exp(-bx)\), with estimates \(\hat a=66/7\) response units and \(\hat b=\log 2\) per unit of \(x\). Suppose the mean-parameter block of the estimated covariance matrix from the fit is diagonal, with \(\widehat{\operatorname{Var}}(\hat a)=1/25\) and \(\widehat{\operatorname{Var}}(\hat b)=1451/108900\); its covariance between \(\hat a\) and \(\hat b\) is zero. These values are hypothetical inputs for illustrating the contrast calculation.

Consider the fitted-value contrast \(c=m(0)-m(1)\). Its estimate is \(\hat c=\hat a(1-\exp(-\hat b))=(66/7)(1-1/2)=33/7\approx4.714\) response units. The gradient with respect to \((a,b)\) is \((1-\exp(-b),a\exp(-b))\); at the estimates it is \((1/2,33/7)\). Therefore,

\widehat{\operatorname{Var}}(\hat c) =\frac{1}{4}\frac{1}{25} +\left(\frac{33}{7}\right)^2\frac{1451}{108900} =\frac{49}{4900}+\frac{1451}{4900} =\frac{15}{49}.

The standard error is \(\sqrt{15}/7\approx0.5533\) response units. A normal-approximation 95% interval is \(33/7\pm1.96\sqrt{15}/7\), or approximately \([3.630,5.799]\) response units. Under the large-sample approximation, this procedure has 95% coverage over repeated samples. It does not mean there is a 95% probability that this particular fixed contrast lies in the computed interval.

Practical Interpretation and a Common Pitfall

A nonlinear Gaussian likelihood fit ties together three parts of the analysis: the assumed shape of the conditional mean, the pattern of conditional variances, and the resulting uncertainty. Changing the variance function can change both the score weights and the fitted mean parameters. It is therefore not generally valid to fit the mean as if all observations had equal variance and then treat that fit as the joint MLE under a parameter-dependent variance model.

Pitfall: dropping the log-variance term. The term \(-\frac12\sum_i\log v_i\) depends on variance parameters. Omitting it can make larger variances appear preferable simply because they reduce the standardized residual penalty.
Pitfall: treating one update as the estimate. A scoring step moves toward a score root locally. The updated value is not automatically an MLE; it must be iterated to convergence and evaluated as a candidate maximum.
Pitfall: using a parameter standard error for a fitted contrast. A contrast generally depends on several parameters. Its variance uses the full gradient and covariance matrix, including covariance terms when they are nonzero.
Pitfall: interpreting the normal interval as a probability statement about a fixed parameter. The interval procedure has repeated-sampling coverage under its assumptions; the fitted interval does not assign probability to a fixed parameter value.

Key takeaway: In nonlinear Gaussian maximum likelihood, derive the likelihood from the complete conditional normal density, include every variance-dependent term, and use the fitted covariance matrix to propagate uncertainty to quantities such as contrasts between predicted means.

Check Your Understanding

Answer each question using the model and methods developed above.

Why must the Gaussian log likelihood include \(\log v_i\) when the variance depends on parameters?
For a Gaussian model, what does the term \(r_i^2/v_i\) measure, and how does it enter the variance score?
In a Fisher-scoring update, why is the parameter increment not necessarily the final MLE?
What matrix is inverted to form the usual observed-information covariance estimate at an interior fitted solution?
For a fitted-value contrast \(c(\theta)\), what gradient and covariance calculation give its delta-method standard error?