Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_univar_estim_normal (g07bb)

## Purpose

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

## Syntax

[xmu, xsig, sexmu, sexsig, corr, dev, nobs, nit, ifail] = g07bb(method, x, xc, ic, xmu, xsig, tol, maxit, 'n', n)
[xmu, xsig, sexmu, sexsig, corr, dev, nobs, nit, ifail] = nag_univar_estim_normal(method, x, xc, ic, xmu, xsig, tol, maxit, 'n', n)

## Description

A sample of size n$n$ is taken from a Normal distribution with mean μ$\mu$ and variance σ2${\sigma }^{2}$ and consists of grouped and/or censored data. Each of the n$n$ observations is known by a pair of values (Li,Ui)$\left({L}_{i},{U}_{i}\right)$ such that:
 Li ≤ xi ≤ Ui. $Li≤xi≤Ui.$
The data is represented as particular cases of this form:
• exactly specified observations occur when Li = Ui = xi${L}_{i}={U}_{i}={x}_{i}$,
• right-censored observations, known only by a lower bound, occur when Ui${U}_{i}\to \infty$,
• left-censored observations, known only by a upper bound, occur when Li${L}_{i}\to -\infty$,
• and interval-censored observations when Li < xi < Ui${L}_{i}<{x}_{i}<{U}_{i}$.
Let the set A$A$ identify the exactly specified observations, sets B$B$ and C$C$ identify the observations censored on the right and left respectively, and set D$D$ identify the observations confined between two finite limits. Also let there be r$r$ exactly specified observations, i.e., the number in A$A$. The probability density function for the standard Normal distribution is
 Z(x) = 1/(sqrt(2π))exp( − (1/2)x2) ,   − ∞ < x < ∞ $Z(x)=12πexp(-12x2) , -∞
and the cumulative distribution function is
 X P(X) = 1 − Q(X) = ∫ Z(x)dx. − ∞
$P(X)= 1-Q(X)=∫-∞XZ(x)dx.$
The log-likelihood of the sample can be written as:
 L (μ,σ) = − r logσ − (1/2) ∑A {(xi − μ) / σ}2 + ∑B log(Q(li)) + ∑C log(P(ui)) + ∑D log(pi) $L (μ,σ) =-r log⁡σ - 1 2 ∑ A { ( x i -μ ) / σ } 2 +∑B log( Q ( l i ) ) + ∑ C log( P ( u i ) ) + ∑ D log( p i )$
where pi = P(ui)P(li)${p}_{i}=P\left({u}_{i}\right)-P\left({l}_{i}\right)$ and ui = (Uiμ) / σ,  li = (Liμ) / σ${u}_{i}=\left({U}_{i}-\mu \right)/\sigma \text{, }{l}_{i}=\left({L}_{i}-\mu \right)/\sigma$.
Let
 S(xi) = (Z(xi))/(Q(xi)),  S1(li,ui) = (Z(li) − Z(ui))/(pi) $S(xi)=Z(xi) Q(xi) , S1(li,ui)=Z(li)-Z(ui)pi$
and
 S2(li,ui) = (uiZ(ui) − liZ(li))/(pi), $S2(li,ui)=uiZ(ui)-liZ(li)pi,$
then the first derivatives of the log-likelihood can be written as:
 ( ∂ L(μ,σ))/( ∂ μ) = L1(μ,σ) = σ − 2∑A(xi − μ) + σ − 1∑BS(li) − σ − 1∑CS( − ui) + σ − 1∑DS1(li,ui) $∂L(μ,σ) ∂μ =L1(μ,σ)=σ-2∑A(xi-μ)+σ-1∑BS(li)-σ-1∑CS(-ui)+σ-1∑DS1(li,ui)$
and
 ( ∂ L(μ,σ))/( ∂ σ) = L2(μ,σ) = − rσ − 1 + σ − 3∑A(xi − μ)2 + σ − 1∑BliS(li) − σ − 1∑CuiS( − ui) $∂L(μ,σ) ∂σ =L2(μ,σ)=-rσ-1+σ-3∑A (xi-μ) 2+σ-1∑BliS(li)-σ-1∑CuiS(-ui)$
 − σ − 1∑DS2(li,ui) $-σ-1∑DS2(li,ui)$
The maximum likelihood estimates, μ̂$\stackrel{^}{\mu }$ and σ̂$\stackrel{^}{\sigma }$, are the solution to the equations:
 L1(μ̂,σ̂) = 0 $L1(μ^,σ^)=0$ (1)
and
 L2(μ̂,σ̂) = 0 $L2(μ^,σ^)=0$ (2)
and if the second derivatives (2 L)/(2 μ) $\frac{{\partial }^{2}L}{{\partial }^{2}\mu }$, (2 L)/(μ σ) $\frac{{\partial }^{2}L}{\partial \mu \partial \sigma }$ and (2L)/(2σ) $\frac{{\partial }^{2}L}{{\partial }^{2}\sigma }$ are denoted by L11${L}_{11}$, L12${L}_{12}$ and L22${L}_{22}$ respectively, then estimates of the standard errors of μ̂$\stackrel{^}{\mu }$ and σ̂$\stackrel{^}{\sigma }$ are given by:
 se(μ̂) = sqrt(( − L22)/(L11L22 − L122)),  se(σ̂) = sqrt(( − L11)/(L11L22 − L122)) $se(μ^)=-L22 L11L22-L122 , se(σ^)=-L11 L11L22-L122$
and an estimate of the correlation of μ̂$\stackrel{^}{\mu }$ and σ̂$\stackrel{^}{\sigma }$ is given by:
 (L12)/(sqrt(L12L22)). $L12L12L22 .$
To obtain the maximum likelihood estimates the equations (1) and (2) can be solved using either the Newton–Raphson method or the Expectation-maximization (EM)$\left(EM\right)$ algorithm of Dempster et al. (1977).
Newton–Raphson Method
This consists of using approximate estimates μ̃$\stackrel{~}{\mu }$ and σ̃$\stackrel{~}{\sigma }$ to obtain improved estimates μ̃ + δμ̃$\stackrel{~}{\mu }+\delta \stackrel{~}{\mu }$ and σ̃ + δσ̃$\stackrel{~}{\sigma }+\delta \stackrel{~}{\sigma }$ by solving
 δμ̃L11 + δσ̃L12 + L1 = 0, δμ̃L12 + δσ̃L22 + L2 = 0,
$δμ~L11+δσ~L12+L1=0, δμ~L12+δσ~L22+L2=0,$
for the corrections δμ̃$\delta \stackrel{~}{\mu }$ and δσ̃$\delta \stackrel{~}{\sigma }$.
EM Algorithm
The expectation step consists of constructing the variable wi${w}_{i}$ as follows:
 if   i ∈ A,   wi = xi E ( ∣ Li < xi < Ui) = μ + σ S1 (li,ui) $if i∈A, wi= xi E (∣ Li< xi< Ui )= μ+σ S1 (li,ui)$ (3)
 if   i ∈ B,   wi = E (xi ∣ xi > Li) = μ + σS (li) S1 (li,ui) $if i∈B, wi= E ( xi∣ xi> Li )=μ+σS (li) S1 (li,ui)$ (4)
 if   i ∈ C,   wi = E (xi ∣ xi < Ui) = μ − σS ( − ui) (li,ui) $if i∈C, wi= E ( xi∣ xi< Ui )=μ-σS (-ui) (li,ui)$ (5)
 if   i ∈ D,   wi = E (xi ∣ Li < xi < Ui) = μ + σ S1 (li,ui) $if i∈D, wi= E ( xi∣ Li< xi< Ui )=μ+σ S1 (li,ui)$ (6)
the maximization step consists of substituting (3), (4), (5) and (6) into (1) and (2) giving:
 n μ̂ = ∑ ŵi / n i = 1
$μ^=∑i=1nw^i/n$
(7)
and
 n σ̂2 = ∑ (ŵi − μ̂)2 / {r + ∑BT(l̂i) + ∑CT( − ûi) + ∑DT1(l̂i,ûi)} i = 1
$σ^2=∑i=1n(w^i-μ^)2/ {r+∑BT(l^i)+∑CT(-u^i)+∑DT1(l^i,u^i)}$
(8)
where
 T(x) = S(x){S(x) − x} ,   T1(l,u) = S12(l,u) + S2(l,u) $T(x)=S(x){S(x)-x} , T1(l,u)=S12(l,u)+S2(l,u)$
and where i${\stackrel{^}{w}}_{i}$, i${\stackrel{^}{l}}_{i}$ and i${\stackrel{^}{u}}_{i}$ are wi${w}_{i}$, li${l}_{i}$ and ui${u}_{i}$ evaluated at μ̂$\stackrel{^}{\mu }$ and σ̂$\stackrel{^}{\sigma }$. Equations (3) to (8) are the basis of the EM$EM$ iterative procedure for finding μ̂$\stackrel{^}{\mu }$ and σ̂2${\stackrel{^}{\sigma }}^{2}$. The procedure consists of alternately estimating μ̂$\stackrel{^}{\mu }$ and σ̂2${\stackrel{^}{\sigma }}^{2}$ using (7) and (8) and estimating {i}$\left\{{\stackrel{^}{w}}_{i}\right\}$ using (3) to (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$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$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$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$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

## Parameters

### Compulsory Input Parameters

1:     method – string (length ≥ 1)
Indicates whether the Newton–Raphson or EM$EM$ algorithm should be used.
If method = 'N'${\mathbf{method}}=\text{'N'}$, then the Newton–Raphson algorithm is used.
If method = 'E'${\mathbf{method}}=\text{'E'}$, then the EM$EM$ algorithm is used.
Constraint: method = 'N'${\mathbf{method}}=\text{'N'}$ or 'E'$\text{'E'}$.
2:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n2${\mathbf{n}}\ge 2$.
The observations xi${x}_{\mathit{i}}$, Li${L}_{\mathit{i}}$ or Ui${U}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
If the observation is exactly specified – the exact value, xi${x}_{i}$.
If the observation is right-censored – the lower value, Li${L}_{i}$.
If the observation is left-censored – the upper value, Ui${U}_{i}$.
If the observation is interval-censored – the lower or upper value, Li${L}_{i}$ or Ui${U}_{i}$, (see xc).
3:     xc(n) – double array
n, the dimension of the array, must satisfy the constraint n2${\mathbf{n}}\ge 2$.
If the j$\mathit{j}$th observation, for j = 1,2,,n$\mathit{j}=1,2,\dots ,n$ is an interval-censored observation then xc(j)${\mathbf{xc}}\left(j\right)$ should contain the complementary value to x(j)${\mathbf{x}}\left(j\right)$, that is, if x(j) < xc(j)${\mathbf{x}}\left(j\right)<{\mathbf{xc}}\left(j\right)$, then xc(j)${\mathbf{xc}}\left(j\right)$ contains upper value, Ui${U}_{i}$, and if x(j) > xc(j)${\mathbf{x}}\left(j\right)>{\mathbf{xc}}\left(j\right)$, then xc(j)${\mathbf{xc}}\left(j\right)$ contains lower value, Li${L}_{i}$. Otherwise if the j$j$th observation is exact or right- or left-censored xc(j)${\mathbf{xc}}\left(j\right)$ need not be set.
Note: if x(j) = xc(j)${\mathbf{x}}\left(j\right)={\mathbf{xc}}\left(j\right)$ then the observation is ignored.
4:     ic(n) – int64int32nag_int array
n, the dimension of the array, must satisfy the constraint n2${\mathbf{n}}\ge 2$.
ic(i)${\mathbf{ic}}\left(\mathit{i}\right)$ contains the censoring codes for the i$\mathit{i}$th observation, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
If ic(i) = 0${\mathbf{ic}}\left(i\right)=0$, the observation is exactly specified.
If ic(i) = 1${\mathbf{ic}}\left(i\right)=1$, the observation is right-censored.
If ic(i) = 2${\mathbf{ic}}\left(i\right)=2$, the observation is left-censored.
If ic(i) = 3${\mathbf{ic}}\left(i\right)=3$, the observation is interval-censored.
Constraint: ic(i) = 0${\mathbf{ic}}\left(\mathit{i}\right)=0$, 1$1$, 2$2$ or 3$3$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
5:     xmu – double scalar
If xsig > 0.0${\mathbf{xsig}}>0.0$ the initial estimate of the mean, μ$\mu$; otherwise xmu need not be set.
6:     xsig – double scalar
Specifies whether an initial estimate of μ$\mu$ and σ$\sigma$ are to be supplied.
xsig > 0.0${\mathbf{xsig}}>0.0$
xsig is the initial estimate of σ$\sigma$ and xmu must contain an initial estimate of μ$\mu$.
xsig0.0${\mathbf{xsig}}\le 0.0$
Initial estimates of xmu and xsig are calculated internally from:
 (a) the exact observations, if the number of exactly specified observations is ≥ 2$\text{}\ge 2$; or (b) the interval-censored observations; if the number of interval-censored observations is ≥ 1$\text{}\ge 1$; or (c) they are set to 0.0$0.0$ and 1.0$1.0$ respectively.
7:     tol – double scalar
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${\mathbf{tol}}=0.0$, then a relative precision of 0.000005$0.000005$ is used.
Constraint: machine precision < tol1.0 or tol = 0.0${\mathbf{tol}}=0.0$.
8:     maxit – int64int32nag_int scalar
The maximum number of iterations.
If maxit0${\mathbf{maxit}}\le 0$, then a value of 25$25$ is used.

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays x, xc, ic. (An error is raised if these dimensions are not equal.)
n$n$, the number of observations.
Constraint: n2${\mathbf{n}}\ge 2$.

wk

### Output Parameters

1:     xmu – double scalar
The maximum likelihood estimate, μ̂$\stackrel{^}{\mu }$, of μ$\mu$.
2:     xsig – double scalar
The maximum likelihood estimate, σ̂$\stackrel{^}{\sigma }$, of σ$\sigma$.
3:     sexmu – double scalar
The estimate of the standard error of μ̂$\stackrel{^}{\mu }$.
4:     sexsig – double scalar
The estimate of the standard error of σ̂$\stackrel{^}{\sigma }$.
5:     corr – double scalar
The estimate of the correlation between μ̂$\stackrel{^}{\mu }$ and σ̂$\stackrel{^}{\sigma }$.
6:     dev – double scalar
The maximized log-likelihood, L(μ̂,σ̂)$L\left(\stackrel{^}{\mu },\stackrel{^}{\sigma }\right)$.
7:     nobs(4$4$) – int64int32nag_int array
The number of the different types of each observation;
nobs(1)${\mathbf{nobs}}\left(1\right)$ contains number of right-censored observations.
nobs(2)${\mathbf{nobs}}\left(2\right)$ contains number of left-censored observations.
nobs(3)${\mathbf{nobs}}\left(3\right)$ contains number of interval-censored observations.
nobs(4)${\mathbf{nobs}}\left(4\right)$ contains number of exactly specified observations.
8:     nit – int64int32nag_int scalar
The number of iterations performed.
9:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, method ≠ 'N'${\mathbf{method}}\ne \text{'N'}$ or 'E'$\text{'E'}$, or n < 2${\mathbf{n}}<2$, or ic(i) ≠ 0${\mathbf{ic}}\left(i\right)\ne 0$, 1$1$, 2$2$ or 3$3$, for some i$i$, or tol < 0.0${\mathbf{tol}}<0.0$, or 0.0 < tol < machine precision, or tol > 1.0${\mathbf{tol}}>1.0$.
ifail = 2${\mathbf{ifail}}=2$
The chosen method failed to converge in maxit iterations. You should either increase tol or maxit or, if using the EM$EM$ algorithm try using the Newton–Raphson method with initial values those returned by the current call to nag_univar_estim_normal (g07bb). All returned values will be reasonable approximations to the correct results if maxit is not very small.
ifail = 3${\mathbf{ifail}}=3$
The chosen method is diverging. This will be due to poor initial values. You should try different initial values.
ifail = 4${\mathbf{ifail}}=4$
nag_univar_estim_normal (g07bb) was unable to calculate the standard errors. This can be caused by the method starting to diverge when the maximum number of iterations was reached.

## Accuracy

The accuracy is controlled by the parameter tol.
If high precision is requested with the EM$EM$ algorithm then there is a possibility that, due to the slow convergence, before the correct solution has been reached the increments of μ̂$\stackrel{^}{\mu }$ and σ̂$\stackrel{^}{\sigma }$ may be smaller than tol and the process will prematurely assume convergence.

## Further Comments

The process is deemed divergent if three successive increments of μ$\mu$ or σ$\sigma$ increase.

## Example

```function nag_univar_estim_normal_example
method = 'N';
x = [4.5;
5.4;
3.9;
5.1;
4.6;
4.8;
2.9;
6.3;
5.5;
4.6;
4.1;
5.2;
3.2;
4;
3.1;
5.1;
3.8;
2.2];
xc = [0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
2.5];
ic = [int64(0);0;0;0;0;0;0;0;0;0;0;0;1;1;1;2;2;3];
xmu = 4;
xsig = 1;
tol = 5e-05;
maxit = int64(50);
[xmuOut, xsigOut, sexmu, sexsig, corr, dev, nobs, nit, ifail] = ...
nag_univar_estim_normal(method, x, xc, ic, xmu, xsig, tol, maxit)
```
```

xmuOut =

4.4924

xsigOut =

1.0196

sexmu =

0.2606

sexsig =

0.1940

corr =

0.0160

dev =

-22.2817

nobs =

3
2
1
12

nit =

5

ifail =

0

```
```function g07bb_example
method = 'N';
x = [4.5;
5.4;
3.9;
5.1;
4.6;
4.8;
2.9;
6.3;
5.5;
4.6;
4.1;
5.2;
3.2;
4;
3.1;
5.1;
3.8;
2.2];
xc = [0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
0;
2.5];
ic = [int64(0);0;0;0;0;0;0;0;0;0;0;0;1;1;1;2;2;3];
xmu = 4;
xsig = 1;
tol = 5e-05;
maxit = int64(50);
[xmuOut, xsigOut, sexmu, sexsig, corr, dev, nobs, nit, ifail] = ...
g07bb(method, x, xc, ic, xmu, xsig, tol, maxit)
```
```

xmuOut =

4.4924

xsigOut =

1.0196

sexmu =

0.2606

sexsig =

0.1940

corr =

0.0160

dev =

-22.2817

nobs =

3
2
1
12

nit =

5

ifail =

0

```

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013