Integer, Intent (In)	::	n, ic(*), maxit
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	nit
Real (Kind=nag_wp), Intent (In)	::	x(n), tol
Real (Kind=nag_wp), Intent (Inout)	::	gamma
Real (Kind=nag_wp), Intent (Out)	::	beta, sebeta, segam, corr, dev, wk(n)
Character (1), Intent (In)	::	cens

C Header Interface

#include <nag.h>

void	g07bef_ (const char cens, const Integer n, const double x[], const Integer ic[], double beta, double gamma, const double tol, const Integer maxit, double sebeta, double segam, double corr, double dev, Integer nit, double wk[], Integer ifail, const Charlen length_cens)

The routine may be called by the names g07bef or nagf_univar_estim_weibull.

3 Description

g07bef computes maximum likelihood estimates of the parameters of the Weibull distribution from exact or right-censored data.

For

n

realizations,

y_{i}

, from a Weibull distribution a value

x_{i}

is observed such that

x_{i} \leq y_{i} .

There are two situations:

(a)exactly specified observations, when $x_{i} = y_{i}$
(b)right-censored observations, known by a lower bound, when $x_{i} < y_{i}$ .

The probability density function of the Weibull distribution, and hence the contribution of an exactly specified observation to the likelihood, is given by:

f (x; λ, γ) = λ γ x^{γ - 1} \exp (- λ x^{γ}), x > 0,   for ​ λ, γ > 0;

while the survival function of the Weibull distribution, and hence the contribution of a right-censored observation to the likelihood, is given by:

S (x; λ, γ) = \exp (- λ x^{γ}), x > 0,   for ​ λ, γ > 0 .

d

of the

n

observations are exactly specified and indicated by

i \in D

and the remaining

(n - d)

are right-censored, then the likelihood function,

Like ​ (λ, γ)

is given by

Like (λ, γ) \propto {(λ γ)}^{d} (\prod_{i \in D} x_{i}^{γ - 1}) \exp (- λ \sum_{i = 1}^{n} x_{i}^{γ}) .

To avoid possible numerical instability a different parameterisation

β, γ

is used, with

β = \log (λ)

. The kernel log-likelihood function,

L (β, γ)

, is then:

L (β, γ) = d \log (γ) + d β + (γ - 1) \sum_{i \in D} \log (x_{i}) - e^{β} \sum_{i = 1}^{n} x_{i}^{γ} .

If the derivatives

\frac{\partial L}{\partial β}

\frac{\partial L}{\partial γ}

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

\frac{\partial^{2} L}{\partial β \partial γ}

and

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

are denoted by

L_{1}

L_{2}

L_{11}

L_{12}

and

L_{22}

, respectively, then 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)

Estimates of the asymptotic 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}}} .

An estimate of the correlation coefficient of

\hat{β}

and

\hat{γ}

is given by:

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

Note: if an estimate of the original parameter

λ

is required, then

\hat{λ} = \exp (\hat{β}) and se (\hat{λ}) = \hat{λ} se (\hat{β}) .

The equations (1) and (2) are solved by the Newton–Raphson iterative method with adjustments made to ensure that

\hat{γ} > 0.0

4 References

Gross A J and Clark V A (1975) Survival Distributions: Reliability Applications in the Biomedical Sciences Wiley

Kalbfleisch J D and Prentice R L (1980) The Statistical Analysis of Failure Time Data Wiley

5 Arguments

1: $cens$ – Character(1) Input

On entry: indicates whether the data is censored or non-censored.

$cens ='N'$: Each observation is assumed to be exactly specified. ic is not referenced.
$cens ='C'$: Each observation is censored according to the value contained in $ic (i)$ , for $i = 1, 2, \dots, n$ .

Constraint:

cens ='N'

'C'

2: $n$ – Integer Input

On entry:

n

, the number of observations.

Constraint:

n \geq 1

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

On entry:

x (i)

contains the

i

th observation,

x_{i}

, for

i = 1, 2, \dots, n

Constraint:

x (i) > 0.0

, for

i = 1, 2, \dots, n

4: $ic (*)$ – Integer array Input

Note: the dimension of the array ic must be at least

n

cens ='C'

, and at least

1

otherwise.

On entry: if

cens ='C'

ic (i)

contains the censoring codes for the

i

th observation, for

i = 1, 2, \dots, n

ic (i) = 0

, the

i

th observation is exactly specified.

ic (i) = 1

, the

i

th observation is right-censored.

cens ='N'

, ic is not referenced.

Constraint: if

cens ='C'

, then

ic (i) = 0

1

, for

i = 1, 2, \dots, n

5: $beta$ – Real (Kind=nag_wp) Output

On exit: the maximum likelihood estimate,

\hat{β}

, of

β

6: $gamma$ – Real (Kind=nag_wp) Input/Output

On entry: indicates whether an initial estimate of

γ

is provided.

gamma > 0.0

, it is taken as the initial estimate of

γ

and an initial estimate of

β

is calculated from this value of

γ

gamma \leq 0.0

, initial estimates of

γ

and

β

are calculated, internally, providing the data contains at least two distinct exact observations. (If there are only two distinct exact observations, the largest observation must not be exactly specified.) See Section 9 for further details.

On exit: contains the maximum likelihood estimate,

\hat{γ}

, of

γ

7: $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 \leq tol \leq 1.0

tol = 0.0

8: $maxit$ – Integer Input

On entry: the maximum number of iterations allowed.

maxit \leq 0

, a value of

25

is used.

9: $sebeta$ – Real (Kind=nag_wp) Output

On exit: an estimate of the standard error of

\hat{β}

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

On exit: an estimate of the standard error of

\hat{γ}

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

On exit: an estimate of the correlation between

\hat{β}

and

\hat{γ}

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

On exit: the maximized kernel log-likelihood,

L (\hat{β}, \hat{γ})

13: $nit$ – Integer Output

On exit: the number of iterations performed.

14: $wk (n)$ – Real (Kind=nag_wp) array Workspace

15: $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, $cens = ⟨ value ⟩$ .
Constraint: $cens ='N'$ or $'C'$ .

On entry, $n = ⟨ value ⟩$ .
Constraint: $n \geq 1$ .

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

$ifail = 2$: On entry, $i = ⟨ value ⟩$ and $ic (i) = ⟨ value ⟩$ .
Constraint: $ic (i) = 0$ or $1$ .

On entry, $i = ⟨ value ⟩$ and $x (i) = ⟨ value ⟩$ .
Constraint: $x (i) > 0.0$ .

$ifail = 3$: On entry, there are no exactly specified observations.

Unable to calculate initial values. This is due to there being either less than two distinct exactly specified observations or exactly two and the largest observation is one of the exact observations.

$ifail = 4$: The chosen method has not converged in $⟨ value ⟩$ iterations. You should either increase tol or maxit.

$ifail = 5$: Hessian matrix of the Newton–Raphson process is singular. Either different initial estimates should be provided or the data should be checked to see if the Weibull distribution is appropriate.

The process has diverged. The process is deemed divergent if three successive increments of $β$ or $γ$ increase. Either different initial estimates should be provided or the data should be checked to see if the Weibull distribution is appropriate.

$ifail = 6$: Potential overflow detected. This is an unlikely error exit usually caused by a large input estimate of $γ$ .

$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

Given that the Weibull distribution is a suitable model for the data and that the initial values are reasonable the convergence to the required accuracy, indicated by tol, should be achieved.

8 Parallelism and Performance

g07bef is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

The initial estimate of

γ

is found by calculating a Kaplan–Meier estimate of the survival function,

\hat{S} (x)

, and estimating the gradient of the plot of

\log (- \log (\hat{S} (x)))

against

x

. This requires the Kaplan–Meier estimate to have at least two distinct points.

The initial estimate of

\hat{β}

, given a value of

\hat{γ}

, is calculated as

\hat{β} = \log (\frac{d}{\sum_{i = 1}^{n} x_{i}^{\hat{γ}}}) .

10 Example

In a study,

20

patients receiving an analgesic to relieve headache pain had the following recorded relief times (in hours):

1.1 1.4 1.3 1.7 1.9 1.8 1.6 2.2 1.7 2.7 4.1 1.8 1.5 1.2 1.4 3.0 1.7 2.3 1.6 2.0

(See Gross and Clark (1975).) This data is read in and a Weibull distribution fitted assuming no censoring; the parameter estimates and their standard errors are printed.

g07be: FL CL CPP AD

NAG FL Interfaceg07bef (estim_​weibull)

▸▿ 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
g07bef (estim_weibull)