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

μ

and variance

σ^{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} \leq x_{i} \leq U_{i} .

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 (x) = \frac{1}{\sqrt{2 π}} \exp (- \frac{1}{2} x^{2}), - \infty < x < \infty

and the cumulative distribution function is

P (X) = 1 - Q (X) = \int_{- \infty}^{X} Z (x) d x .

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

L (μ, σ) = - r \log σ - \frac{1}{2} \sum_{A} {\{(x_{i} - μ) / σ\}}^{2} + \sum_{B} \log (Q (l_{i})) + \sum_{C} \log (P (u_{i})) + \sum_{D} \log (p_{i})

where

p_{i} = P (u_{i}) - P (l_{i})

and

u_{i} = (U_{i} - μ) / σ, l_{i} = (L_{i} - μ) / σ

Let

S (x_{i}) = \frac{Z (x_{i})}{Q (x_{i})}, S_{1} (l_{i}, u_{i}) = \frac{Z (l_{i}) - Z (u_{i})}{p_{i}}

and

S_{2} (l_{i}, u_{i}) = \frac{u_{i} Z (u_{i}) - l_{i} Z (l_{i})}{p_{i}},

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

\frac{\partial L (μ, σ)}{\partial μ} = L_{1} (μ, σ) = σ^{- 2} \sum_{A} (x_{i} - μ) + σ^{- 1} \sum_{B} S (l_{i}) - σ^{- 1} \sum_{C} S (- u_{i}) + σ^{- 1} \sum_{D} S_{1} (l_{i}, u_{i})

and

\frac{\partial L (μ, σ)}{\partial σ} = L_{2} (μ, σ) = - r σ^{- 1} + σ^{- 3} \sum_{A} {(x_{i} - μ)}^{2} + σ^{- 1} \sum_{B} l_{i} S (l_{i}) - σ^{- 1} \sum_{C} u_{i} S (- u_{i})

- σ^{- 1} \sum_{D} S_{2} (l_{i}, u_{i})

The maximum likelihood estimates,

\hat{μ}

and

\hat{σ}

, are the solution to the equations:

L_{1} (\hat{μ}, \hat{σ}) = 0

(1)

and

L_{2} (\hat{μ}, \hat{σ}) = 0

(2)

and if the second derivatives

\frac{\partial^{2} L}{\partial^{2} μ}

\frac{\partial^{2} L}{\partial μ \partial σ}

and

\frac{\partial^{2} L}{\partial^{2} σ}

are denoted by

L_{11}

L_{12}

and

L_{22}

respectively, then estimates of the standard errors of

\hat{μ}

and

\hat{σ}

are given by:

se (\hat{μ}) = \sqrt{\frac{- L_{22}}{L_{11} L_{22} - L_{12}^{2}}}, se (\hat{σ}) = \sqrt{\frac{- L_{11}}{L_{11} L_{22} - L_{12}^{2}}}

and an estimate of the correlation of

\hat{μ}

and

\hat{σ}

is given by:

\frac{L_{12}}{\sqrt{L_{12} L_{22}}} .

To obtain the maximum likelihood estimates the equations (1) and (2) can be solved using either the Newton–Raphson method or the Expectation-maximization

(E M)

algorithm of Dempster et al. (1977).

Newton–Raphson Method

This consists of using approximate estimates

\tilde{μ}

and

\tilde{σ}

to obtain improved estimates

\tilde{μ} + δ \tilde{μ}

and

\tilde{σ} + δ \tilde{σ}

by solving

\begin{array}{l} δ \tilde{μ} L_{11} + δ \tilde{σ} L_{12} + L_{1} = 0, \\ δ \tilde{μ} L_{12} + δ \tilde{σ} L_{22} + L_{2} = 0, \end{array}

for the corrections

δ \tilde{μ}

and

δ \tilde{σ}

EM Algorithm

The expectation step consists of constructing the variable

w_{i}

as follows:

if i \in A, w_{i} = x_{i}

(3)

if i \in B, w_{i} = E (x_{i} ∣ x_{i} > L_{i}) = μ + σ S (l_{i})

(4)

if i \in C, w_{i} = E (x_{i} ∣ x_{i} < U_{i}) = μ - σ S (- u_{i})

(5)

if i \in D, w_{i} = E (x_{i} ∣ L_{i} < x_{i} < U_{i}) = μ + σ S_{1} (l_{i}, u_{i})

(6)

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

\hat{μ} = \sum_{i = 1}^{n} {\hat{w}}_{i} / n

(7)

and

{\hat{σ}}^{2} = \sum_{i = 1}^{n} {({\hat{w}}_{i} - \hat{μ})}^{2} / \{r + \sum_{B} T ({\hat{l}}_{i}) + \sum_{C} T (- {\hat{u}}_{i}) + \sum_{D} T_{1} ({\hat{l}}_{i}, {\hat{u}}_{i})\}

(8)

where

T (x) = S (x) \{S (x) - x\}, 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{μ}

and

\hat{σ}

. Equations (3) to (8) are the basis of the

E M

iterative procedure for finding

\hat{μ}

and

{\hat{σ}}^{2}

. The procedure consists of alternately estimating

\hat{μ}

and

{\hat{σ}}^{2}

using (7) and (8) and estimating

\{{\hat{w}}_{i}\}

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

E M

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

E M

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

E M

algorithm should be considered.

References

Dempster A P, Laird N M and Rubin D B (1977) Maximum likelihood from incomplete data via the

E M

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

E M

algorithm should be used.

method ='N'

, then the Newton–Raphson algorithm is used.

method ='E'

, then the

E M

algorithm is used.

Constraint:

method ='N'

'E'

2: $x (n)$ – double array

The observations

x_{i}

L_{i}

U_{i}

, for

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}

U_{i}

, (see xc).

3: $xc (n)$ – double array

If the

j

th observation, for

j = 1, 2, \dots, n

is an interval-censored observation then

xc (j)

should contain the complementary value to

x (j)

, that is, if

x (j) < xc (j)

, then

xc (j)

contains upper value,

U_{i}

, and if

x (j) > xc (j)

, then

xc (j)

contains lower value,

L_{i}

. Otherwise if the

j

th observation is exact or right- or left-censored

xc (j)

need not be set.

Note: if

x (j) = xc (j)

then the observation is ignored.

4: $ic (n)$ – int64int32nag_int array

ic (i)

contains the censoring codes for the

i

th observation, for

i = 1, 2, \dots, n

ic (i) = 0

, the observation is exactly specified.

ic (i) = 1

, the observation is right-censored.

ic (i) = 2

, the observation is left-censored.

ic (i) = 3

, the observation is interval-censored.

Constraint:

ic (i) = 0

1

2

3

, for

i = 1, 2, \dots, n

5: $xmu$ – double scalar

xsig > 0.0

the initial estimate of the mean,

μ

; otherwise xmu need not be set.

6: $xsig$ – double scalar

Specifies whether an initial estimate of

μ

and

σ

are to be supplied.

$xsig > 0.0$

xsig is the initial estimate of

σ

and xmu must contain an initial estimate of

μ

$xsig \leq 0.0$

Initial estimates of xmu and xsig are calculated internally from:

(a)	the exact observations, if the number of exactly specified observations is $\geq 2$ ; or
(b)	the interval-censored observations; if the number of interval-censored observations is $\geq 1$ ; or
(c)	they are set to $0.0$ and $1.0$ respectively.

7: $tol$ – double scalar

The relative precision required for the final estimates of

μ

and

σ

. Convergence is assumed when the absolute relative changes in the estimates of both

μ

and

σ

are less than tol.

tol = 0.0

, then a relative precision of

0.000005

is used.

Constraint:

machine precision < tol \leq 1.0

tol = 0.0

8: $maxit$ – int64int32nag_int scalar

The maximum number of iterations.

maxit \leq 0

, then a value of

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$ , the number of observations.

Constraint: $n \geq 2$ .

Output Parameters

1: $xmu$ – double scalar: The maximum likelihood estimate, $\hat{μ}$ , of $μ$ .
2: $xsig$ – double scalar: The maximum likelihood estimate, $\hat{σ}$ , of $σ$ .
3: $sexmu$ – double scalar: The estimate of the standard error of $\hat{μ}$ .
4: $sexsig$ – double scalar: The estimate of the standard error of $\hat{σ}$ .
5: $corr$ – double scalar: The estimate of the correlation between $\hat{μ}$ and $\hat{σ}$ .
6: $dev$ – double scalar: The maximized log-likelihood, $L (\hat{μ}, \hat{σ})$ .
7: $nobs (4)$ – int64int32nag_int array: The number of the different types of each observation;
$nobs (1)$ contains number of right-censored observations.

$nobs (2)$ contains number of left-censored observations.

$nobs (3)$ contains number of interval-censored observations.

$nobs (4)$ contains number of exactly specified observations.
8: $nit$ – int64int32nag_int scalar: The number of iterations performed.
9: $ifail$ – int64int32nag_int scalar: $ifail = 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$

On entry,	$method \neq'N'$ or $'E'$ ,
or	$n < 2$ ,
or	$ic (i) \neq 0$ , $1$ , $2$ or $3$ , for some $i$ ,
or	$tol < 0.0$ ,
or	$0.0 < tol < machine precision$ ,
or	$tol > 1.0$ .

$ifail = 2$: The chosen method failed to converge in maxit iterations. You should either increase tol or maxit or, if using the $E M$ 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$: The chosen method is diverging. This will be due to poor initial values. You should try different initial values.

$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.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

The accuracy is controlled by the argument tol.

If high precision is requested with the

E M

algorithm then there is a possibility that, due to the slow convergence, before the correct solution has been reached the increments of

\hat{μ}

and

\hat{σ}

may be smaller than tol and the process will prematurely assume convergence.

Further Comments

The process is deemed divergent if three successive increments of

μ

σ

increase.

Example

A sample of

18

observations and their censoring codes are read in and the Newton–Raphson method used to compute the estimates.

Open in the MATLAB editor: g07bb_example

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

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox