Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_univar_estim_weibull (g07be)

Purpose

nag_univar_estim_weibull (g07be) computes maximum likelihood estimates for parameters of the Weibull distribution from data which may be right-censored.

Syntax

[beta, gamma, sebeta, segam, corr, dev, nit, ifail] = g07be(cens, x, ic, gamma, tol, maxit, 'n', n)
[beta, gamma, sebeta, segam, corr, dev, nit, ifail] = nag_univar_estim_weibull(cens, x, ic, gamma, tol, maxit, 'n', n)

Description

nag_univar_estim_weibull (g07be) computes maximum likelihood estimates of the parameters of the Weibull distribution from exact or right-censored data.
For n$n$ realisations, yi${y}_{i}$, from a Weibull distribution a value xi${x}_{i}$ is observed such that
 xi ≤ yi. $xi≤yi.$
There are two situations:
 (a) exactly specified observations, when xi = yi${x}_{i}={y}_{i}$ (b) right-censored observations, known by a lower bound, when xi < yi${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γ − 1exp( − λxγ),  x > 0,   for ​λ,γ > 0; $f(x;λ,γ)=λγxγ-1exp(-λ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. $S(x;λ,γ)=exp(-λ xγ), x> 0, for ​ λ ,γ> 0.$
If d$d$ of the n$n$ observations are exactly specified and indicated by iD$i\in D$ and the remaining (nd)$\left(n-d\right)$ are right-censored, then the likelihood function, Like ​(λ,γ)$\text{Like ​}\left(\lambda ,\gamma \right)$ is given by
Like(λ,γ) ∝ (λγ)d (∏i ∈ Dxiγ − 1)
 exp ( n ) − λ ∑ xiγ i = 1
.
$Like(λ,γ)∝(λγ)d (∏i∈Dxiγ-1) exp(-λ∑i=1nxiγ) .$
To avoid possible numerical instability a different parameterisation β,γ$\beta ,\gamma$ is used, with β = log(λ)$\beta =\mathrm{log}\left(\lambda \right)$. The kernel log-likelihood function, L(β,γ)$L\left(\beta ,\gamma \right)$, is then:
 n L(β,γ) = dlog(γ) + dβ + (γ − 1)∑i ∈ Dlog(xi) − eβ ∑ xiγ. i = 1
$L(β,γ)=dlog(γ)+dβ+(γ-1)∑i∈Dlog(xi)-eβ∑i=1nxiγ.$
If the derivatives (L)/(β) $\frac{\partial L}{\partial \beta }$, (L)/(γ) $\frac{\partial L}{\partial \gamma }$, (2L)/(β2) $\frac{{\partial }^{2}L}{{\partial \beta }^{2}}$, (2L)/(β γ) $\frac{{\partial }^{2}L}{\partial \beta \partial \gamma }$ and (2L)/(γ2) $\frac{{\partial }^{2}L}{{\partial \gamma }^{2}}$ are denoted by L1${L}_{1}$, L2${L}_{2}$, L11${L}_{11}$, L12${L}_{12}$ and L22${L}_{22}$, respectively, then the maximum likelihood estimates, β̂$\stackrel{^}{\beta }$ and γ̂$\stackrel{^}{\gamma }$, are the solution to the equations:
 L1(β̂,γ̂) = 0 $L1(β^,γ^)=0$ (1)
and
 L2(β̂,γ̂) = 0 $L2(β^,γ^)=0$ (2)
Estimates of the asymptotic standard errors of β̂$\stackrel{^}{\beta }$ and γ̂$\stackrel{^}{\gamma }$ are given by:
 se(β̂) = sqrt(( − L22)/(L11L22 − L122)),  se(γ̂) = sqrt(( − L11)/(L11L22 − L122)). $se(β^)=-L22 L11L22-L122 , se(γ^)=-L11 L11L22-L122 .$
An estimate of the correlation coefficient of β̂$\stackrel{^}{\beta }$ and γ̂$\stackrel{^}{\gamma }$ is given by:
 (L12)/(sqrt(L12L22)). $L12L12L22 .$
Note:  if an estimate of the original parameter λ$\lambda$ is required, then
 λ̂ = exp(β̂)  and  se(λ̂) = λ̂se(β̂). $λ^=exp(β^) and se(λ^)=λ^se(β^).$
The equations (1) and (2) are solved by the Newton–Raphson iterative method with adjustments made to ensure that γ̂ > 0.0$\stackrel{^}{\gamma }>0.0$.

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

Parameters

Compulsory Input Parameters

1:     cens – string (length ≥ 1)
Indicates whether the data is censored or non-censored.
cens = 'N'${\mathbf{cens}}=\text{'N'}$
Each observation is assumed to be exactly specified. ic is not referenced.
cens = 'C'${\mathbf{cens}}=\text{'C'}$
Each observation is censored according to the value contained in ic(i)${\mathbf{ic}}\left(\mathit{i}\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
Constraint: cens = 'N'${\mathbf{cens}}=\text{'N'}$ or 'C'$\text{'C'}$.
2:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
x(i)${\mathbf{x}}\left(\mathit{i}\right)$ contains the i$\mathit{i}$th observation, xi${x}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
Constraint: x(i) > 0.0${\mathbf{x}}\left(\mathit{i}\right)>0.0$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
3:     ic( : $:$) – int64int32nag_int array
Note: the dimension of the array ic must be at least n${\mathbf{n}}$ if cens = 'C'${\mathbf{cens}}=\text{'C'}$, and at least 1$1$ otherwise.
If cens = 'C'${\mathbf{cens}}=\text{'C'}$, then ic(i)${\mathbf{ic}}\left(\mathit{i}\right)$ contains the censoring codes for the i$\mathit{i}$th observation, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
If ic(i) = 0${\mathbf{ic}}\left(i\right)=0$, the i$i$th observation is exactly specified.
If ic(i) = 1${\mathbf{ic}}\left(i\right)=1$, the i$i$th observation is right-censored.
If cens = 'N'${\mathbf{cens}}=\text{'N'}$, then ic is not referenced.
Constraint: if cens = 'C'${\mathbf{cens}}=\text{'C'}$, then ic(i) = 0${\mathbf{ic}}\left(\mathit{i}\right)=0$ or 1$1$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
4:     gamma – double scalar
Indicates whether an initial estimate of γ$\gamma$ is provided.
If gamma > 0.0${\mathbf{gamma}}>0.0$, it is taken as the initial estimate of γ$\gamma$ and an initial estimate of β$\beta$ is calculated from this value of γ$\gamma$.
If gamma0.0${\mathbf{gamma}}\le 0.0$, then initial estimates of γ$\gamma$ and β$\beta$ 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 [Further Comments] for further details.
5:     tol – double scalar
The relative precision required for the final estimates of β$\beta$ and γ$\gamma$. Convergence is assumed when the absolute relative changes in the estimates of both β$\beta$ and γ$\gamma$ are less than tol.
If tol = 0.0${\mathbf{tol}}=0.0$, then a relative precision of 0.000005$0.000005$ is used.
Constraint: machine precisiontol1.0 or tol = 0.0${\mathbf{tol}}=0.0$.
6:     maxit – int64int32nag_int scalar
The maximum number of iterations allowed.
If maxit0${\mathbf{maxit}}\le 0$, then a value of 25$25$ is used.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the array x.
n$n$, the number of observations.
Constraint: n1${\mathbf{n}}\ge 1$.

wk

Output Parameters

1:     beta – double scalar
The maximum likelihood estimate, β̂$\stackrel{^}{\beta }$, of β$\beta$.
2:     gamma – double scalar
Contains the maximum likelihood estimate, γ̂$\stackrel{^}{\gamma }$, of γ$\gamma$.
3:     sebeta – double scalar
An estimate of the standard error of β̂$\stackrel{^}{\beta }$.
4:     segam – double scalar
An estimate of the standard error of γ̂$\stackrel{^}{\gamma }$.
5:     corr – double scalar
An estimate of the correlation between β̂$\stackrel{^}{\beta }$ and γ̂$\stackrel{^}{\gamma }$.
6:     dev – double scalar
The maximized kernel log-likelihood, L(β̂,γ̂)$L\left(\stackrel{^}{\beta },\stackrel{^}{\gamma }\right)$.
7:     nit – int64int32nag_int scalar
The number of iterations performed.
8:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{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${\mathbf{ifail}}=1$
 On entry, cens ≠ 'N'${\mathbf{cens}}\ne \text{'N'}$ or 'C'$\text{'C'}$, or n < 1${\mathbf{n}}<1$, or tol < 0.0${\mathbf{tol}}<0.0$, or 0.0 < tol < machine precision, or tol > 1.0${\mathbf{tol}}>1.0$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, the i$i$th observation, x(i) ≤ 0.0${\mathbf{x}}\left(i\right)\le 0.0$, for some i = 1,2, … ,n$i=1,2,\dots ,n$, or the i$i$th censoring code, ic(i) ≠ 0${\mathbf{ic}}\left(i\right)\ne 0$ or 1$1$, for some i = 1,2, … ,n$i=1,2,\dots ,n$ and cens = 'C'${\mathbf{cens}}=\text{'C'}$.
ifail = 3${\mathbf{ifail}}=3$
On entry, there are no exactly specified observations, or the function 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.
ifail = 4${\mathbf{ifail}}=4$
The method has failed to converge in maxit iterations. You should increase tol or maxit.
ifail = 5${\mathbf{ifail}}=5$
Process has diverged. The process is deemed divergent if three successive increments of β$\beta$ or γ$\gamma$ 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.
ifail = 6${\mathbf{ifail}}=6$
A potential overflow has been detected. This is an unlikely exit usually caused by a large input estimate of γ$\gamma$.

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.

The initial estimate of γ$\gamma$ is found by calculating a Kaplan–Meier estimate of the survival function, (x)$\stackrel{^}{S}\left(x\right)$, and estimating the gradient of the plot of log(log((x)))$\mathrm{log}\left(-\mathrm{log}\left(\stackrel{^}{S}\left(x\right)\right)\right)$ against x$x$. This requires the Kaplan–Meier estimate to have at least two distinct points.
The initial estimate of β̂$\stackrel{^}{\beta }$, given a value of γ̂$\stackrel{^}{\gamma }$, is calculated as
 β̂ = log(d/( ∑ i = 1nxiγ̂)) . $β^=log(d∑i=1nxiγ^ ) .$

Example

function nag_univar_estim_weibull_example
cens = 'No censor';
x = [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;
1.7;
2.3;
1.6;
2];
ic = [int64(0)];
gamma = 0;
tol = 0;
maxit = int64(0);
[beta, gammaOut, sebeta, segam, corr, dev, nit, ifail] = ...
nag_univar_estim_weibull(cens, x, ic, gamma, tol, maxit)

beta =

-2.1073

gammaOut =

2.7870

sebeta =

0.4627

segam =

0.4273

corr =

-0.8755

dev =

-20.5864

nit =

5

ifail =

0

function g07be_example
cens = 'No censor';
x = [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;
1.7;
2.3;
1.6;
2];
ic = [int64(0)];
gamma = 0;
tol = 0;
maxit = int64(0);
[beta, gammaOut, sebeta, segam, corr, dev, nit, ifail] = g07be(cens, x, ic, gamma, tol, maxit)

beta =

-2.1073

gammaOut =

2.7870

sebeta =

0.4627

segam =

0.4273

corr =

-0.8755

dev =

-20.5864

nit =

5

ifail =

0