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

2 Specification

Fortran Interface

Subroutine g07bbf (

method, n, x, xc, ic, xmu, xsig, tol, maxit, sexmu, sexsig, corr, dev, nobs, nit, wk, ifail)

Integer, Intent (In)	::	n, ic(n), maxit
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	nobs(4), nit
Real (Kind=nag_wp), Intent (In)	::	x(n), xc(n), tol
Real (Kind=nag_wp), Intent (Inout)	::	xmu, xsig
Real (Kind=nag_wp), Intent (Out)	::	sexmu, sexsig, corr, dev, wk(2*n)
Character (1), Intent (In)	::	method

C Header Interface

#include <nag.h>

void

g07bbf_ (const char *method, const Integer *n, const double x[], const double xc[], const Integer ic[], double *xmu, double *xsig, const double *tol, const Integer *maxit, double *sexmu, double *sexsig, double *corr, double *dev, Integer nobs[], Integer *nit, double wk[], Integer *ifail, const Charlen length_method)

The routine may be called by the names g07bbf or nagf_univar_estim_normal.

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$ – Character(1) Input

On entry: indicates whether the Newton–Raphson or

E M

algorithm should be used.

method ='N'

, the Newton–Raphson algorithm is used.

method ='E'

, the

E M

algorithm is used.

Constraint:

method ='N'

'E'

2: $n$ – Integer Input

On entry:

n

, the number of observations.

Constraint:

n \geq 2

3: $x (n)$ – Real (Kind=nag_wp) array Input

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)$ – Real (Kind=nag_wp) array Input

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

5: $ic (n)$ – Integer array Input

On entry:

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

6: $xmu$ – Real (Kind=nag_wp) 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$ – Real (Kind=nag_wp) 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$ – Real (Kind=nag_wp) Input

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

, a relative precision of

0.000005

is used.

Constraint:

machine precision < tol \leq 1.0

tol = 0.0

9: $maxit$ – Integer Input

On entry: the maximum number of iterations.

maxit \leq 0

, a value of

25

is used.

10: $sexmu$ – Real (Kind=nag_wp) Output

On exit: the estimate of the standard error of

\hat{μ}

11: $sexsig$ – Real (Kind=nag_wp) Output

On exit: the estimate of the standard error of

\hat{σ}

12: $corr$ – Real (Kind=nag_wp) Output

On exit: the estimate of the correlation between

\hat{μ}

and

\hat{σ}

13: $dev$ – Real (Kind=nag_wp) Output

On exit: the maximized log-likelihood,

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

14: $nobs (4)$ – Integer array Output

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

15: $nit$ – Integer Output

On exit: the number of iterations performed.

16: $wk (2 \times n)$ – Real (Kind=nag_wp) array Workspace

17: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

- 1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

- 1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

- 1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

0

is recommended. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: On entry, effective number of observations $< 2$ .

On entry, $i = 〈value〉$ and $ic (i) = 〈value〉$ .
Constraint: $ic (i) = 0$ , $1$ , $2$ or $3$ .

On entry, $method = 〈value〉$ .
Constraint: $method ='N'$ or $'E'$ .

On entry, $n = 〈value〉$ .
Constraint: $n \geq 2$ .

On entry, $tol = 〈value〉$ .
Constraint: $machine precision < tol \leq 1.0$ or $tol = 0.0$ .

$ifail = 2$: The chosen method has not converged in $〈value〉$ 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 g07bbf. All returned values will be reasonable approximations to the correct results if maxit is not very small.

$ifail = 3$: The EM process has failed. Different initial values should be tried.

The process has diverged. Different initial values should be tried.

$ifail = 4$: Standard errors cannot be computed. 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.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

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

g07bbf is not threaded in any implementation.

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.

g07bb: FL CL CPP AD

NAG FL Interfaceg07bbf (estim_​normal)

▸▿ 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 FL Interface
g07bbf (estim_normal)