$ontext

  Maximum likelihood estimation of weibull distribution

  Erwin Kalvelagen, 2008

  References:

  SAS/IML User's Guide, Example 11.4: MLEs for Two-Parameter Weibull Distribution

  Lawless, J. F. (1982), Statistical Models and Methods For Lifetime Data, New York: Wiley

  Pike, M. C. [1966]. A suggested method of analysis of a certain class of experiments
  in carcinogenesis. Biometrics 22, 142-61


  Data:
  Number of days until the appearance of a carcinoma in 19
  rats painted with carcinogen DMBA.

$offtext


set i /i1*i19/;
table data(i,*)
             days   censored
    i1        143
    i2        164
    i3        188
    i4        188
    i5        190
    i6        192
    i7        206
    i8        209
    i9        213
    i10       216
    i11       220
    i12       227
    i13       230
    i14       234
    i15       246
    i16       265
    i17       304
    i18       216     1
    i19       244     1
;

set k(i) 'not censored';
k(i)$(data(i,'censored')=0) = yes;

parameter x(i);
x(i) = data(i,'days');

scalars
  p 'number of observations'
  m 'number of uncensored observations'
;

p = card(i);
m = card(k);

display p,m;


*-------------------------------------------------------------------------------
* estimation
*-------------------------------------------------------------------------------

scalar theta 'location parameter' /0/;

variables
     sigma 'scale parameter'
     c     'shape parameter'
     loglik 'log likelihood'
;

equation eloglike;

c.lo = 0.001;
sigma.lo = 0.001;


eloglike.. loglik =e= m*log(c) - m*c*log(sigma)
                 + (c-1)*sum(k,log(x(k)-theta))
                 - sum(i,((x(i)-theta)/sigma)**c);

model mle /eloglike/;
solve mle maximizing loglik using nlp;

*-------------------------------------------------------------------------------
* get hessian
*-------------------------------------------------------------------------------
option nlp=convert;
$onecho > convert.opt
hessian
$offecho
mle.optfile=1;
solve mle minimizing loglik using nlp;


*
* gams cannot add elements at runtime so we declare the necessary elements here
*
set dummy /e1,x1,x2/;

parameter h(*,*,*) '-hessian';
execute_load "hessian.gdx",h;
display h;

set j /sigma,c/;
parameter h0(j,j);
h0('sigma','sigma') = h('e1','x1','x1');
h0('c','c') = h('e1','x2','x2');
h0('sigma','c') = h('e1','x1','x2');
h0('c','sigma') = h('e1','x1','x2');
display h0;


*-------------------------------------------------------------------------------
* invert hessian
*-------------------------------------------------------------------------------

execute_unload "h.gdx",j,h0;
execute "=invert.exe h.gdx j h0 invh.gdx invh";
parameter invh(j,j);
execute_load "invh.gdx",invh;
display invh;


*-------------------------------------------------------------------------------
* quantile of normal distribution
*-------------------------------------------------------------------------------

* find
*   p = 0.05
*   q = probit(1-p/2)

scalar prob /.05/;

* we don't have the inverse error function so we calculate it
* using a small cns model
equations e;
variables probit;
e.. Errorf(probit) =e= 1-prob/2;
model inverterrorf /e/;
solve inverterrorf using cns;

display probit.l;

* verification:
*> qnorm(0.975);
*[1] 1.959964
*>

*-------------------------------------------------------------------------------
* calculate standard errors and confidence intervals
*-------------------------------------------------------------------------------

parameter result(j,*);
result('c','estimate') = c.l;
result('sigma','estimate') = sigma.l;

result(j,'stderr') = sqrt(abs(invh(j,j)));

result(j,'conf lo') = result(j,'estimate') - probit.l*result(j,'stderr');
result(j,'conf up') = result(j,'estimate') + probit.l*result(j,'stderr');

display result;