[R] Problem with Newton_Raphson

Christopher Kelvin chris_kelvin2001 at yahoo.com
Thu Sep 20 13:46:32 CEST 2012


Hello,
I have being trying to estimate the parameters of the generalized exponential distribution. The random number generation for the GE distribution is x<-(-log(1-U^(1/p1))/b), where U stands for uniform dist. The data i have generated to estimate the parameters is right censored and the code is given below; The problem is that, the newton-Raphson approach isnt working and i do not know what is wrong. Can somebody suggest something or help in identifying what the problem might be. 

p1<-0.6;b<-2
n=20;rr=5000
U<-runif(n,0,1)
for (i in 1:rr){

x<-(-log(1-U^(1/p1))/b)

 meantrue<-gamma(1+(1/p1))*b
  meantrue
  d<-meantrue/0.30
  cen<- runif(n,min=0,max=d)
  s<-ifelse(x<=cen,1,0)
  q<-c(x,cen)

    z<-function(data, p){ 
    shape<-p[1]
    scale<-p[2]
    log1<-n*sum(s)*log(p[1])+ n*sum(s)*log(p[2])+(p[1]-1)*sum(s)*log(1-((exp(-(p[2])*sum(x)))))
-(p[2])*sum(t) + (p[1])*log((exp(-(p[2])*sum(x))))-
(p[1])*sum(s)*log((exp(-(p[2])*sum(x))))
  return(-log1)
  }
}
  start <- c(1,1)
  zz<-optim(start,fn=z,data=q,hessian=T)
  zz
  m1<-zz$par[2]
  p<-zz$par[1] 


Thank you
Chris Kelvin
INSPEM. UPM




More information about the R-help mailing list