naginterfaces.library.univar.estim_​normal

naginterfaces.library.univar.estim_normal(method, x, xc, ic, xmu, xsig, tol, maxit)

estim_normal computes maximum likelihood estimates and their standard errors for parameters of the Normal distribution from grouped and/or censored data.

For full information please refer to the NAG Library document for g07bb

https://support.nag.com/numeric/nl/nagdoc_30/flhtml/g07/g07bbf.html

Parameters

methodstr, length 1

Indicates whether the Newton–Raphson or $EM$ algorithm should be used.

If $method=\text{'N'}$, the Newton–Raphson algorithm is used.

If $method=\text{'E'}$, the $EM$ algorithm is used.

xfloat, array-like, shape $\left(n\right)$

The observations $x_{\text{i}}$, $L_{\text{i}}$ or $U_{\text{i}}$, for $\text{i}=1,2,\dots ,n$.

If the observation is exactly specified – the exact value, $x_i$.

If the observation is right-censored – the lower value, $L_i$.

If the observation is left-censored – the upper value, $U_i$.

If the observation is interval-censored – the lower or upper value, $L_i$ or $U_i$, (see $xc$).

xcfloat, array-like, shape $\left(n\right)$

If the $\text{j}$th observation, for $\text{j}=1,2,\dots ,n$ is an interval-censored observation then $xc[j-1]$ should contain the complementary value to $x[j-1]$, that is, if $x[j-1]<xc[j-1]$, then $xc[j-1]$ contains upper value, $U_i$, and if $x[j-1]>xc[j-1]$, then $xc[j-1]$ contains lower value, $L_i$. Otherwise if the $j$th observation is exact or right - or left-censored $xc[j-1]$ need not be set.

Note: if $x[j-1]=xc[j-1]$ then the observation is ignored.

icint, array-like, shape $\left(n\right)$

$ic[\text{i}-1]$ contains the censoring codes for the $\text{i}$th observation, for $\text{i}=1,2,\dots ,n$.

If $ic[i-1]=0$, the observation is exactly specified.

If $ic[i-1]=1$, the observation is right-censored.

If $ic[i-1]=2$, the observation is left-censored.

If $ic[i-1]=3$, the observation is interval-censored.

xmufloat

If $xsig>0.0$ the initial estimate of the mean, $\mu$; otherwise $xmu$ need not be set.

xsigfloat

Specifies whether an initial estimate of $\mu$ and $\sigma$ are to be supplied.

$xsig>0.0$

$xsig$ is the initial estimate of $\sigma$ and $xmu$ must contain an initial estimate of $\mu$.

$xsig\le 0.0$

Initial estimates of $xmu$ and $xsig$ are calculated internally from:

1. the exact observations, if the number of exactly specified observations is $\ge 2$; or

2. the interval-censored observations; if the number of interval-censored observations is $\ge 1$; or

3. they are set to $0.0$ and $1.0$ respectively.

tolfloat

The relative precision required for the final estimates of $\mu$ and $\sigma$. Convergence is assumed when the absolute relative changes in the estimates of both $\mu$ and $\sigma$ are less than $tol$.

If $tol=0.0$, a relative precision of $0.000005$ is used.

maxitint

The maximum number of iterations.

If $maxit\le 0$, a value of $25$ is used.

Returns

xmufloat

The maximum likelihood estimate, $\hat{\mu }$, of $\mu$.

xsigfloat

The maximum likelihood estimate, $\hat{\sigma }$, of $\sigma$.

sexmufloat

The estimate of the standard error of $\hat{\mu }$.

sexsigfloat

The estimate of the standard error of $\hat{\sigma }$.

corrfloat

The estimate of the correlation between $\hat{\mu }$ and $\hat{\sigma }$.

devfloat

The maximized log-likelihood, $L(\hat{\mu },\hat{\sigma })$.

nobsint, ndarray, shape $\left(4\right)$

The number of the different types of each observation;

$nobs\left[0\right]$ contains number of right-censored observations.

$nobs\left[1\right]$ contains number of left-censored observations.

$nobs\left[2\right]$ contains number of interval-censored observations.

$nobs\left[3\right]$ contains number of exactly specified observations.

nitint

The number of iterations performed.

Raises

NagValueError

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 2$.

(errno $1$)

On entry, $method=⟨value⟩$.

Constraint: $method=\text{'N'}$ or $\text{'E'}$.

(errno $1$)

On entry, effective number of observations $<2$.

(errno $1$)

On entry, $i=⟨value⟩$ and $ic[i-1]=⟨value⟩$.

Constraint: $ic[i-1]=0$, $1$, $2$ or $3$.

(errno $1$)

On entry, $tol=⟨value⟩$.

Constraint: $\text{machine precision}<tol\le 1.0$ or $tol=0.0$.

(errno $2$)

The chosen method has not converged in $⟨value⟩$ iterations.

(errno $3$)

The process has diverged.

(errno $3$)

The EM process has failed.

(errno $4$)

Standard errors cannot be computed.

Notes

A sample of size $n$ is taken from a Normal distribution with mean $\mu$ and variance ${\sigma }^2$ and consists of grouped and/or censored data. Each of the $n$ observations is known by a pair of values $(L_i,U_i)$ such that:

$$L_i\le x_i\le U_i\text{.}$$

The data is represented as particular cases of this form:

exactly specified observations occur when $L_i=U_i=x_i$,

right-censored observations, known only by a lower bound, occur when $U_i\to \infty$,

left-censored observations, known only by a upper bound, occur when $L_i\to -\infty$,

and interval-censored observations when $L_i<x_i<U_i$.

Let the set $A$ identify the exactly specified observations, sets $B$ and $C$ identify the observations censored on the right and left respectively, and set $D$ identify the observations confined between two finite limits. Also let there be $r$ exactly specified observations, i.e., the number in $A$. The probability density function for the standard Normal distribution is

$$Z\left(x\right)=\frac{1}{\sqrt{2\pi }}exp(-\frac{1}{2}{x}^2)\text{,}-\infty <x<\infty$$

and the cumulative distribution function is

$$P\left(X\right)=1-Q\left(X\right)={\int }_{-\infty }^{X}Z\left(x\right)dx\text{.}$$

The log-likelihood of the sample can be written as:

$$L(\mu ,\sigma )=-r\log \left(\sigma \right)-\frac{1}{2}\sum _A{\left\{(x_i-\mu )/\sigma \right\}}^2+\sum _B\log \left(Q\left(l_i\right)\right)+\sum _C\log \left(P\left(u_i\right)\right)+\sum _D\log \left(p_i\right)$$

where $p_i=P\left(u_i\right)-P\left(l_i\right)$ and $u_i=(U_i-\mu )/\sigma \text{,}{l}_i=(L_i-\mu )/\sigma$.

Let

$$S\left(x_i\right)=\frac{Z\left(x_i\right)}{Q\left(x_i\right)}\text{,}{S}_1(l_i,u_i)=\frac{Z\left(l_i\right)-Z\left(u_i\right)}{p_i}$$

and

$$S_2(l_i,u_i)=\frac{u_{i}Z\left(u_i\right)-l_{i}Z\left(l_i\right)}{p_i}\text{,}$$

then the first derivatives of the log-likelihood can be written as:

$$\frac{\partial L(\mu ,\sigma )}{\partial \mu }=L_1(\mu ,\sigma )={\sigma }^{-2}\sum _A(x_i-\mu )+{\sigma }^{-1}\sum _{B}S\left(l_i\right)-{\sigma }^{-1}\sum _{C}S(-u_i)+{\sigma }^{-1}\sum _D{S}_1(l_i,u_i)$$

and

$$\frac{\partial L(\mu ,\sigma )}{\partial \sigma }=L_2(\mu ,\sigma )=-r{\sigma }^{-1}+{\sigma }^{-3}\sum _A{(x_i-\mu )}^2+{\sigma }^{-1}\sum _B{l}_{i}S\left(l_i\right)-{\sigma }^{-1}\sum _C{u}_{i}S(-u_i)$$

$$-{\sigma }^{-1}\sum _D{S}_2(l_i,u_i)$$

The maximum likelihood estimates, $\hat{\mu }$ and $\hat{\sigma }$, are the solution to the equations:

$$L_1(\hat{\mu },\hat{\sigma })=0$$

and

$$L_2(\hat{\mu },\hat{\sigma })=0$$

and if the second derivatives $\frac{{\partial }^{2}L}{{\partial }^2\mu }$, $\frac{{\partial }^{2}L}{\partial \mu \partial \sigma }$ and $\frac{{\partial }^{2}L}{{\partial }^2\sigma }$ are denoted by $L_{11}$, $L_{12}$ and $L_{22}$ respectively, then estimates of the standard errors of $\hat{\mu }$ and $\hat{\sigma }$ are given by:

$$se\left(\hat{\mu }\right)=\sqrt{\frac{-L_{22}}{L_{11}{L}_{22}-L_{12}^2}}\text{,}se\left(\hat{\sigma }\right)=\sqrt{\frac{-L_{11}}{L_{11}{L}_{22}-L_{12}^2}}$$

and an estimate of the correlation of $\hat{\mu }$ and $\hat{\sigma }$ is given by:

$$\frac{L_{12}}{\sqrt{L_{12}{L}_{22}}}\text{.}$$

To obtain the maximum likelihood estimates the equations (1) and (2) can be solved using either the Newton–Raphson method or the Expectation-maximization $\left(EM\right)$ algorithm of Dempster et al. (1977).

Newton–Raphson Method

This consists of using approximate estimates $\tilde{\mu }$ and $\tilde{\sigma }$ to obtain improved estimates $\tilde{\mu }+\delta \tilde{\mu }$ and $\tilde{\sigma }+\delta \tilde{\sigma }$ by solving

$$\begin{array}{c}\begin{array}{c}\delta \tilde{\mu }{L}_{11}+\delta \tilde{\sigma }{L}_{12}+L_1=0\text{,}\\ \\ \delta \tilde{\mu }{L}_{12}+\delta \tilde{\sigma }{L}_{22}+L_2=0\text{,}\end{array}\end{array}$$

for the corrections $\delta \tilde{\mu }$ and $\delta \tilde{\sigma }$.

EM Algorithm

The expectation step consists of constructing the variable $w_i$ as follows:

$$\text{if}i\in A\text{,}{w}_i=x_i$$

$$\text{if}i\in B\text{,}{w}_i=E(x_i|x_i>L_i)=\mu +\sigma S\left(l_i\right)$$

$$\text{if}i\in C\text{,}{w}_i=E(x_i|x_i<U_i)=\mu -\sigma S(-u_i)$$

$$\text{if}i\in D\text{,}{w}_i=E(x_i|L_i<x_i<U_i)=\mu +\sigma S_1(l_i,u_i)$$

the maximization step consists of substituting (3), (4), (5) and (6) into (1) and (2) giving:

$$\hat{\mu }=\sum _{i=1}^n{\hat{w}}_i/n$$

and

$${\hat{\sigma }}^2=\sum _{i=1}^n{({\hat{w}}_i-\hat{\mu })}^2/\{r+\sum _{B}T\left({\hat{l}}_i\right)+\sum _{C}T(-{\hat{u}}_i)+\sum _D{T}_1({\hat{l}}_i,{\hat{u}}_i)\}$$

where

$$T\left(x\right)=S\left(x\right)\{S\left(x\right)-x\}\text{,}{T}_1(l,u)=S_1^2(l,u)+S_2(l,u)$$

and where ${\hat{w}}_i$, ${\hat{l}}_i$ and ${\hat{u}}_i$ are $w_i$, $l_i$ and $u_i$ evaluated at $\hat{\mu }$ and $\hat{\sigma }$. Equations (3) and (8) are the basis of the $EM$ iterative procedure for finding $\hat{\mu }$ and ${\hat{\sigma }}^2$. The procedure consists of alternately estimating $\hat{\mu }$ and ${\hat{\sigma }}^2$ using (7) and (8) and estimating $\left\{{\hat{w}}_i\right\}$ using (3) and (6).

In choosing between the two methods a general rule is that the Newton–Raphson method converges more quickly but requires good initial estimates whereas the $EM$ algorithm converges slowly but is robust to the initial values. In the case of the censored Normal distribution, if only a small proportion of the observations are censored then estimates based on the exact observations should give good enough initial estimates for the Newton–Raphson method to be used. If there are a high proportion of censored observations then the $EM$ algorithm should be used and if high accuracy is required the subsequent use of the Newton–Raphson method to refine the estimates obtained from the $EM$ algorithm should be considered.

References

Dempster, A P, Laird, N M and Rubin, D B, 1977, Maximum likelihood from incomplete data via the $EM$ algorithm (with discussion), J. Roy. Statist. Soc. Ser. B (39), 1–38

Swan, A V, 1969, Algorithm AS 16. Maximum likelihood estimation from grouped and censored normal data, Appl. Statist. (18), 110–114

Wolynetz, M S, 1979, Maximum likelihood estimation from confined and censored normal data, Appl. Statist. (28), 185–195