Message-ID: <1348147023.81022.YahooMailNeo@web120603.mail.ne1.yahoo.com>
Date: 2012-09-20T13:17:03Z
From: Christopher Kelvin
Subject: Problem with Newton_Raphson
In-Reply-To: <6415E33E-5A3B-4EF9-BE5C-C30D7A24EA20@xs4all.nl>
By your suggestion if i understand you, i have changed (p[1]) and (p[2]) and have also corrected the error sum(t), but if i run it, the parameters estimates are very large.
Can you please run it and help me out? The code is given below.
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(shape)+ n*sum(s)*log(scale)+(shape-1)*sum(s)*log(1-((exp(-(scale)*sum(x)))))
-(scale)*sum(x) + (shape)*log((exp(-(scale)*sum(x))))-
(shape)*sum(s)*log((exp(-(scale)*sum(x))))
? return(-log1)
? }
? start <- c(1,1)
? zz<-optim(start,fn=z,data=q,hessian=T)
? zz
Thank you
Chris Kelvin
----- Original Message -----
From: Berend Hasselman <bhh at xs4all.nl>
To: Christopher Kelvin <chris_kelvin2001 at yahoo.com>
Cc: "r-help at r-project.org" <r-help at r-project.org>
Sent: Thursday, September 20, 2012 8:52 PM
Subject: Re: [R] Problem with Newton_Raphson
On 20-09-2012, at 13:46, Christopher Kelvin wrote:
> 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.
>
Newton-Raphson? is not a method for optim.
> 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]
>
Running your code given above gives an error message:
Error in sum(t) : invalid 'type' (closure) of argument
Where is object 't'?
Why are you defining function z within the rr loop? Only the last definition is given to optim.
Why use p[1] and p[2] explicitly in the calculation of log1 in the body of function z when you can use shape and scale defined in the lines before log1 <-.
Berend