NAG Library Routine Document
G07BEF computes maximum likelihood estimates for parameters of the Weibull distribution from data which may be right-censored.
|SUBROUTINE G07BEF (
||CENS, N, X, IC, BETA, GAMMA, TOL, MAXIT, SEBETA, SEGAM, CORR, DEV, NIT, WK, IFAIL)
||N, IC(*), MAXIT, NIT, IFAIL
||X(N), BETA, GAMMA, TOL, SEBETA, SEGAM, CORR, DEV, WK(N)
G07BEF computes maximum likelihood estimates of the parameters of the Weibull distribution from exact or right-censored data.
, from a Weibull distribution a value
is observed such that
There are two situations:
||exactly specified observations, when
||right-censored observations, known by a lower bound, when .
The probability density function of the Weibull distribution, and hence the contribution of an exactly specified observation to the likelihood, is given by:
while the survival function of the Weibull distribution, and hence the contribution of a right-censored observation to the likelihood, is given by:
observations are exactly specified and indicated by
and the remaining
are right-censored, then the likelihood function,
is given by
To avoid possible numerical instability a different parameterisation
is used, with
. The kernel log-likelihood function,
, is then:
If the derivatives
are denoted by
, respectively, then the maximum likelihood estimates,
, are the solution to the equations:
Estimates of the asymptotic standard errors of
are given by:
An estimate of the correlation coefficient of
is given by:
if an estimate of the original parameter
is required, then
The equations (1)
are solved by the Newton–Raphson iterative method with adjustments made to ensure that
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
- 1: CENS – CHARACTER(1)Input
: indicates whether the data is censored or non-censored.
- Each observation is assumed to be exactly specified. IC is not referenced.
- Each observation is censored according to the value contained in
, for .
- 2: N – INTEGERInput
On entry: , the number of observations.
- 3: X(N) – REAL (KIND=nag_wp) arrayInput
On entry: contains the th observation, , for .
, for .
- 4: IC() – INTEGER arrayInput
the dimension of the array IC
must be at least
, and at least
contains the censoring codes for the
th observation, for
If , the th observation is exactly specified.
If , the th observation is right-censored.
, then IC
is not referenced.
if , then or , for .
- 5: BETA – REAL (KIND=nag_wp)Output
On exit: the maximum likelihood estimate, , of .
- 6: GAMMA – REAL (KIND=nag_wp)Input/Output
: indicates whether an initial estimate of
If , it is taken as the initial estimate of and an initial estimate of is calculated from this value of .
, then initial estimates of
are calculated, internally, providing the data contains at least two distinct exact observations. (If there are only two distinct exact observations, then the largest observation must not be exactly specified.) See Section 8
for further details.
On exit: contains the maximum likelihood estimate, , of .
- 7: TOL – REAL (KIND=nag_wp)Input
: the relative precision required for the final estimates of
. Convergence is assumed when the absolute relative changes in the estimates of both
are less than TOL
If , then a relative precision of is used.
- 8: MAXIT – INTEGERInput
: the maximum number of iterations allowed.
If , then a value of is used.
- 9: SEBETA – REAL (KIND=nag_wp)Output
On exit: an estimate of the standard error of .
- 10: SEGAM – REAL (KIND=nag_wp)Output
On exit: an estimate of the standard error of .
- 11: CORR – REAL (KIND=nag_wp)Output
On exit: an estimate of the correlation between and .
- 12: DEV – REAL (KIND=nag_wp)Output
On exit: the maximized kernel log-likelihood, .
- 13: NIT – INTEGEROutput
On exit: the number of iterations performed.
- 14: WK(N) – REAL (KIND=nag_wp) arrayWorkspace
- 15: IFAIL – INTEGERInput/Output
must be set to
. If you are unfamiliar with this parameter you should refer to Section 3.3
in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
. When the value is used it is essential to test the value of IFAIL on exit.
unless the routine detects an error or a warning has been flagged (see Section 6
6 Error Indicators and Warnings
If on entry
, explanatory error messages are output on the current error message unit (as defined by X04AAF
Errors or warnings detected by the routine:
|On entry,|| or ,|
|On entry,||the th observation, , for some ,|
|or||the th censoring code, or , for some and .|
On entry, there are no exactly specified observations, or the routine was requested to calculate initial values and there are either less than two distinct exactly specified observations or there are exactly two and the largest observation is one of the exact observations.
The method has failed to converge in MAXIT
iterations. You should increase TOL
Process has diverged. The process is deemed divergent if three successive increments of or increase or if the 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.
A potential overflow has been detected. This is an unlikely exit usually caused by a large input estimate of .
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.
The initial estimate of is found by calculating a Kaplan–Meier estimate of the survival function, , and estimating the gradient of the plot of against . This requires the Kaplan–Meier estimate to have at least two distinct points.
The initial estimate of
, given a value of
, is calculated as
In a study,
patients receiving an analgesic to relieve headache pain had the following recorded relief times (in hours):
(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.
9.1 Program Text
Program Text (g07befe.f90)
9.2 Program Data
Program Data (g07befe.d)
9.3 Program Results
Program Results (g07befe.r)