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

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 arguments 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$ 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 $\left({L}_{i},{U}_{i}\right)$ such that:
 $Li≤xi≤Ui.$
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
 $Zx=12πexp-12x2 , -∞
and the cumulative distribution function is
 $PX= 1-QX=∫-∞XZxdx.$
The log-likelihood of the sample can be written as:
 $L μ,σ =-r log⁡σ - 1 2 ∑ A x i -μ / σ 2 +∑B log Q l i + ∑ C log P u i + ∑ D log p i$
where ${p}_{i}=P\left({u}_{i}\right)-P\left({l}_{i}\right)$ and ${u}_{i}=\left({U}_{i}-\mu \right)/\sigma \text{, }{l}_{i}=\left({L}_{i}-\mu \right)/\sigma$.
Let
 $Sxi=Zxi Qxi , S1li,ui=Zli-Zuipi$
and
 $S2li,ui=uiZui-liZlipi,$
then the first derivatives of the log-likelihood can be written as:
 $∂Lμ,σ ∂μ =L1μ,σ=σ-2∑Axi-μ+σ-1∑BSli-σ-1∑CS-ui+σ-1∑DS1li,ui$
and
 $∂Lμ,σ ∂σ =L2μ,σ=-rσ-1+σ-3∑A xi-μ 2+σ-1∑BliSli-σ-1∑CuiS-ui$
 $-σ-1∑DS2li,ui$
The maximum likelihood estimates, $\stackrel{^}{\mu }$ and $\stackrel{^}{\sigma }$, are the solution to the equations:
 $L1μ^,σ^=0$ (1)
and
 $L2μ^,σ^=0$ (2)
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 $\stackrel{^}{\mu }$ and $\stackrel{^}{\sigma }$ are given by:
 $seμ^=-L22 L11L22-L122 , seσ^=-L11 L11L22-L122$
and an estimate of the correlation of $\stackrel{^}{\mu }$ and $\stackrel{^}{\sigma }$ is given by:
 $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 $\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,$
for the corrections $\delta \stackrel{~}{\mu }$ and $\delta \stackrel{~}{\sigma }$.
EM Algorithm
The expectation step consists of constructing the variable ${w}_{i}$ as follows:
 $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$ (4)
 $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$ (6)
the maximization step consists of substituting (3), (4), (5) and (6) into (1) and (2) giving:
 $μ^=∑i=1nw^i/n$ (7)
and
 $σ^2=∑i=1nw^i-μ^2/ r+∑BTl^i+∑CT-u^i+∑DT1l^i,u^i$ (8)
where
 $Tx=SxSx-x , T1l,u=S12l,u+S2l,u$
and where ${\stackrel{^}{w}}_{i}$, ${\stackrel{^}{l}}_{i}$ and ${\stackrel{^}{u}}_{i}$ are ${w}_{i}$, ${l}_{i}$ and ${u}_{i}$ evaluated at $\stackrel{^}{\mu }$ and $\stackrel{^}{\sigma }$. Equations (3) to (8) are the basis of the $EM$ iterative procedure for finding $\stackrel{^}{\mu }$ and ${\stackrel{^}{\sigma }}^{2}$. The procedure consists of alternately estimating $\stackrel{^}{\mu }$ and ${\stackrel{^}{\sigma }}^{2}$ using (7) and (8) and estimating $\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$ 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

## Parameters

### Compulsory Input Parameters

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

### Optional Input Parameters

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

### Output Parameters

1:     $\mathrm{xmu}$ – double scalar
The maximum likelihood estimate, $\stackrel{^}{\mu }$, of $\mu$.
2:     $\mathrm{xsig}$ – double scalar
The maximum likelihood estimate, $\stackrel{^}{\sigma }$, of $\sigma$.
3:     $\mathrm{sexmu}$ – double scalar
The estimate of the standard error of $\stackrel{^}{\mu }$.
4:     $\mathrm{sexsig}$ – double scalar
The estimate of the standard error of $\stackrel{^}{\sigma }$.
5:     $\mathrm{corr}$ – double scalar
The estimate of the correlation between $\stackrel{^}{\mu }$ and $\stackrel{^}{\sigma }$.
6:     $\mathrm{dev}$ – double scalar
The maximized log-likelihood, $L\left(\stackrel{^}{\mu },\stackrel{^}{\sigma }\right)$.
7:     $\mathrm{nobs}\left(4\right)$int64int32nag_int array
The number of the different types of each observation;
${\mathbf{nobs}}\left(1\right)$ contains number of right-censored observations.
${\mathbf{nobs}}\left(2\right)$ contains number of left-censored observations.
${\mathbf{nobs}}\left(3\right)$ contains number of interval-censored observations.
${\mathbf{nobs}}\left(4\right)$ contains number of exactly specified observations.
8:     $\mathrm{nit}$int64int32nag_int scalar
The number of iterations performed.
9:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{method}}\ne \text{'N'}$ or $\text{'E'}$, or ${\mathbf{n}}<2$, or ${\mathbf{ic}}\left(i\right)\ne 0$, $1$, $2$ or $3$, for some $i$, or ${\mathbf{tol}}<0.0$, or , or ${\mathbf{tol}}>1.0$.
${\mathbf{ifail}}=2$
The chosen method failed to converge in maxit iterations. You should either increase tol or maxit or, if using the $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.
${\mathbf{ifail}}=3$
The chosen method is diverging. This will be due to poor initial values. You should try different initial values.
${\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.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The accuracy is controlled by the argument tol.
If high precision is requested with the $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.

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

## Example

A sample of $18$ observations and their censoring codes are read in and the Newton–Raphson method used to compute the estimates.
```function g07bb_example

fprintf('g07bb example results\n\n');

% Data
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];

n  = numel(x);
xc = zeros(n,1);
ic = zeros(n,1,'int64');
xc(n,1) = 2.5;
ic(n-5:n,1) = [1; 1; 1; 2; 2; 3];

% Parameters
method = 'N';
xmu  = 4;
xsig = 1;
tol  = 5e-05;
maxit = int64(50);

% Calculate estimates
[xmu, xsig, sexmu, sexsig, corr, dev, nobs, nit, ifail] = ...
g07bb( ...
method, x, xc, ic, xmu, xsig, tol, maxit);

% Display results
fprintf(' Mean                                     = %8.4f\n', xmu);
fprintf(' Standard deviation                       = %8.4f\n', xsig);
fprintf(' Standard error of mean                   = %8.4f\n', sexmu);
fprintf(' Standard error of sigma                  = %8.4f\n', sexsig);
fprintf(' Correlation coefficient                  = %8.4f\n', corr);
fprintf(' Number of right censored observations    = %3d\n', nobs(1));
fprintf(' Number of left censored observations     = %3d\n', nobs(2));
fprintf(' Number of interval censored observations = %3d\n', nobs(3));
fprintf(' Number of exactly specified observations = %3d\n', nobs(4));
fprintf(' Number of iterations                     = %3d\n', nit);
fprintf(' Log-likelihood                           = %8.4f\n', dev);

```
```g07bb example results

Mean                                     =   4.4924
Standard deviation                       =   1.0196
Standard error of mean                   =   0.2606
Standard error of sigma                  =   0.1940
Correlation coefficient                  =   0.0160
Number of right censored observations    =   3
Number of left censored observations     =   2
Number of interval censored observations =   1
Number of exactly specified observations =  12
Number of iterations                     =   5
Log-likelihood                           = -22.2817
```