Prev: Help w/ Difference-in-Difference
Next: Datastep question
From: JamesZHENG on 14 Apr 2010 15:28
There is something wrong with the following code. But I don't know
where the error is.

The beta1 and beta2 value calculated from the code do not match with
the values got from other software packages and procedures.

Please help me with this. Thanks!




options charcode nocenter linesize=80;
data president;
input obs dead survival;
cards;

1 1 10.0
2 1 29.0
3 1 26.0
4 1 28.0
5 1 15.0
6 1 23.0
7 1 17.0
8 1 25.0
9 1 0.1
10 1 20.0
11 1 4.0
12 1 1.0
13 1 24.0
14 1 16.0
15 1 12.0
16 1 4.0
17 1 10.0
18 1 17.0
19 1 16.0
20 1 0.1
21 1 7.0
22 1 24.0
23 1 12.0
24 1 4.0
25 1 18.0
26 1 21.0
27 1 11.0
28 1 2.0
29 1 9.0
30 1 36.0
31 1 12.0
32 1 28.0
33 1 3.0
34 1 16.0
35 1 9.0
;


proc iml;
start mle;
use president;
read all into y var{dead};
read all into z var{survival};
print, "Newton-Raphson Estimation of Weibull Model";

/* Starting values for beta1 and beta2 */
beta = {0.5,0.5};
b1 = beta[1,];
b2 = beta [2,];
x = y||z;
crit = 1;
n=35;
result = j(15,7,0);


/* Use Newton's method to maximize likelihood function.*/

*** the following defines the log-likelihood function for the Weibull
distribution;
/* b1=theta (scale parameter); b2= beta (location parameter);
thet - location parameter, setting to 0 */
do iter = 1 to 15 while (crit > 1.0e-10);
p = 35; m = p - 0;
sum1 = 0.; sum2 = 0.;
do i = 1 to p;
temp = x[i];
if i <= m then sum1 = sum1 + log(temp);
sum2 = sum2 + (temp /b1)##b2;
end;
logL = m*log(b2) - m*b2*log(b1) + (b2-1)*sum1 - sum2;


***The following code specifies a gradient module,
i.e. the derivatives respect to theta and beta;

g = j(1,2,0.);
p = 35; m = p - 0;
sum1 = 0.; sum2 = 0.; sum3 = 0.;
do i = 1 to p;
temp = x[i];
if i <= m then sum1 = sum1 + log(temp);
sum2 = sum2 + (temp / b1)##b2;
sum3 = sum3 + ((temp / b1)##b2) * (log(temp / b1));
end;
g[1] = -m * b2 / b1 + sum2 * b2 / b1;
g[2] = m / b2 - m * log(b1) + sum1 - sum3;



/* ***The following code specifies elements of The Hessian matrix,
i.e. the second derivatives with respect to theta and beta */

h = j(2,2,0.);
p = 35; m = p - 0;
sum1 = 0.; sum2 = 0.; sum3 = 0.; sum4 = 0.;sum5 =0.;
do i = 1 to p;
temp = x[i];
if i <= m then
sum1 = sum1 + log(temp);
sum2 = sum2 + (temp / b1)##b2;
sum3 = sum3 + ((temp / b1)##b2) * (log(temp / b1));
sum4 = sum4 + b2* ((temp / b1)##b2) * (log(temp / b1));
sum5 = sum5 + ((temp / b1)##b2) * ((log(temp / b1))##2);
end;
h[1,1] = (b2 / b1##2) *(m -(b2-1)* sum2);
h[1,2] = (1 / b2) * ( m - sum2 - sum4);
h[2,1] = (1 / b2) * ( m - sum2 - sum4);
h[2,2] = (m / b2##2) - sum5;


/* The Newton two step */
gt = T(g);
db= -inv(h)*gt;
betanew = beta + db;
crit = sqrt (ssq(betanew - beta));
beta = betanew;

/* Place the result of this iteration in the result matrix */
result[iter,] = iter||g||T(beta)||crit||logL;

end;

cnames = {iter,g1,g2,beta1,beta2,crit,logl};
print "Results of Newton Iterations" result [colname=cnames];

cov = -inv(h); * cov from Hessian;

print "Beta Coefficients" beta;
print "Hessian Form of Covariance Matrix" cov;

finish;
run mle;
quit;
 | 
Pages: 1
Prev: Help w/ Difference-in-Difference
Next: Datastep question