NAG Library Function Document

nag_censored_normal (g07bbc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

+− 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1 Purpose

nag_censored_normal (g07bbc) computes maximum likelihood estimates and their standard errors for arguments of the Normal distribution from grouped and/or censored data.

2 Specification

#include <nag.h>

#include <nagg07.h>

void	nag_censored_normal (Nag_CEMethod method, Integer n, const double x[], const double xc[], const Integer ic[], double xmu, double xsig, double tol, Integer maxit, double sexmu, double sexsig, double corr, double dev, Integer nobs[], Integer nit, NagError fail)

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

4 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

5 Arguments

1: method – Nag_CEMethodInput

On entry: indicates whether the Newton–Raphson or

E M

algorithm should be used.

method = Nag_CE_NR

, then the Newton–Raphson algorithm is used.

method = Nag_CE_EM

, then the

E M

algorithm is used.

Constraint:

method = Nag_CE_NR

Nag_CE_EM

2: n – IntegerInput

On entry:

n

, the number of observations.

Constraint:

n \geq 2

3: x[n] – const doubleInput

On entry: 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).

4: xc[n] – const doubleInput

On entry: if the

j

th observation, for

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.

5: ic[n] – const IntegerInput

On entry:

ic [i - 1]

contains the censoring codes for the

i

th observation, for

i = 1, 2, \dots, n

ic [i - 1] = 0

, the observation is exactly specified.

ic [i - 1] = 1

, the observation is right-censored.

ic [i - 1] = 2

, the observation is left-censored.

ic [i - 1] = 3

, the observation is interval-censored.

Constraint:

ic [i - 1] = 0

1

2

3

, for

i = 1, 2, \dots, n

6: xmu – double *Input/Output

On entry: if

xsig > 0.0

the initial estimate of the mean,

μ

; otherwise xmu need not be set.

On exit: the maximum likelihood estimate,

\hat{μ}

, of

μ

7: xsig – double *Input/Output

On entry: 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.

On exit: the maximum likelihood estimate,

\hat{σ}

, of

σ

8: tol – doubleInput

On entry: 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

9: maxit – IntegerInput

On entry: the maximum number of iterations.

maxit \leq 0

, then a value of

25

is used.

10: sexmu – double *Output

On exit: the estimate of the standard error of

\hat{μ}

11: sexsig – double *Output

On exit: the estimate of the standard error of

\hat{σ}

12: corr – double *Output

On exit: the estimate of the correlation between

\hat{μ}

and

\hat{σ}

13: dev – double *Output

On exit: the maximized log-likelihood,

L (\hat{μ}, \hat{σ})

14: nobs[ $4$ ] – IntegerOutput

On exit: the number of the different types of each observation;

nobs [0]

contains number of right-censored observations.

nobs [1]

contains number of left-censored observations.

nobs [2]

contains number of interval-censored observations.

nobs [3]

contains number of exactly specified observations.

15: nit – Integer *Output

On exit: the number of iterations performed.

16: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $⟨value⟩$ had an illegal value.
NE_CONVERGENCE: Method has not converged in $⟨value⟩$ iterations.
NE_DIVERGENCE: Process has diverged.
NE_EM_PROCESS: The EM process has failed.
NE_INT: On entry, $n = ⟨value⟩$ .
Constraint: $n \geq 2$ .
NE_INT_ARRAY: On entry, $ic [⟨value⟩]$ is invalid, it contains $⟨value⟩$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_OBSERVATIONS: On entry, effective number of observations $< 2$ .
NE_REAL: On entry, tol is invalid: $tol = ⟨value⟩$ .
NE_STANDARD_ERRORS: Standard errors cannot be computed.

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

8 Parallelism and Performance

Not applicable.

9 Further Comments

The process is deemed divergent if three successive increments of

μ

σ

increase.

10 Example

A sample of

18

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

NAG Library Function Documentnag_censored_normal (g07bbc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_censored_normal (g07bbc)