牛顿法

Aim: Find oˆ such that

Problem: Analytic solution of likelihood equations not always available.

 Example: Censored exponentially distributed observations

 Suppose that  and that the censored times

 

are observed. Let m be the number of uncensored observations. Then

with first and second derivative

Thus we obtain for the observed and expected information

Thus the MLE can be obtained be the Newton-Raphson iteration

Numerical example: Choose starting value in (0, 1)

Implementation in R:

[html]  view plain copy

1. #Statistics 24600 - Spring 2004  
2. #Instructor: Michael Eichler  
3. #  
4. #Method : Newton-Raphson method  
5. #Example: Exponential distribution  
6. #----------------------------------  
7. #Log-likelihood with first and second derivative  
8. ln<-function(p,Y,R) {  
9.   m<-sum(R==1)  
10.   ln<-m*log(p)-p*sum(Y)  
11.   attr(ln,"gradient")<-m/p-sum(Y)  
12.   attr(ln,"hessian")<--m/p^2  
13.   ln  
14. }  
15. #Newton-Raphson method  
16. newmle<-function(p,ln) {  
17.   l<-ln(p)  
18.   pnew<-p-attr(l,"gradient")/attr(l,"hessian")  
19.   pnew  
20. }  
21. #Simulate censored exponentially distributed data  
22. Y<-rexp(10,1/5)  
23. R<-ifelse(Y>10,0,1)  
24. Y[R==0]=10  
25. #Plot first derivative of the log-likelihood  
26. x<-seq(0.05,0.6,0.01)  
27. plot(x,attr(ln(x,Y,R),"gradient"),type="l",  
28.   xlab=expression(theta),ylab="Score function")  
29. abline(0,0)  
30. #Apply Newton-Raphson iteration 3 times  
31. p<-newmle(p,ln,Y=Y,R=R)  
32. p  
33. p<-newmle(p,ln,Y=Y,R=R)  
34. p  
35. p<-newmle(p,ln,Y=Y,R=R)  
36. p