Your function is giving NaN's during the optimization.
The R-forge version of optimx() has functionality specifically intended to deal
with this.
NOTE: the CRAN version does not, and the R-forge version still has some
glitches!
However, I easily ran the code you supplied by changing optim to optimx in the
penultimate
line. Here's the final output.
KKT condition testing
Number of parameters = 2
max abs gradient element = 0.004032794 test tol = 0.1018794
KKT1 result = TRUE
Hessian eigenvalues:
[1] 7.138974e+02 9.931285e-04
KKT2 result = TRUE
KKT results: gmax= 0.004032794 evratio= 1.391136e-06 KKT1 & 2: TRUE
TRUE
[1] 7.138974e+02 9.931285e-04
Save results from method BFGS
> zz
method par fvalues fns grs hes rs conv KKT1 KKT2
2 BFGS 0.2832468, 52.6096161 100.8794 10 1 0 -2 3 TRUE TRUE
1 NM 0.2827563, 54.3551491 100.8779 53 0 0 1 0 TRUE TRUE
mtilt xtimes meths
2 NA 0.11 BFGS
1 0.0005130808 0.57 NM
There were 50 or more warnings (use warnings() to see the first
50)>
Note the Hessian eigenvalues are pretty bad. And the warnings are the NaN's
produced.
This is with just the 2 default methods. When I tried "all.methods",
one (uobyqa) seemed
to lock up. This is a fairly ill-conditioned problem.
Best, JN
On 04/14/2012 06:00 AM, r-help-request at r-project.org
wrote:> Message: 13
> Date: Fri, 13 Apr 2012 03:54:43 -0700 (PDT)
> From: Christopher Kelvin <chris_kelvin2001 at yahoo.com>
> To: "r-help at r-project.org" <r-help at r-project.org>
> Subject: [R] R-Help: Censoring data
> Message-ID:
> <1334314483.27693.YahooMailNeo at web65408.mail.ac4.yahoo.com>
> Content-Type: text/plain; charset=iso-8859-1
>
> Hello,
> ?I want to estimate weibull parameters with 30% censored data. I have below
the code for the censoring
> ?but how it must be put into the likelihood equation to obtain the desire
estimate is where i have a problem with,
> ?can some body help?
> ?My likelihood equation is for a random type-I censoring where time for the
censored units is different for each unit.
> ?
> n=50;r=35
> p=0.8;b=1.5
> t<-rweibull(50,shape=p,scale=b)
> meantrue<-gamma(1+(1/p))*b
> meantrue
> d<-meantrue/0.30
>
> cen<- runif(50,min=0,max=d)
> cen
> s<-ifelse(t<=cen,1,0)
> s
>
> z<-function(p){?
> shape<-p[1]
> scale<-p[2]
>
log1<-(r*log(p[1])-r*(p[1])*log(p[2])+(p[1]-1)*sum(log(t))-sum((t/(p[2]))^(p[1])
> )-((n-r)*(sum(cen)/(p[2]))^(p[1])))
> return(-log1)
> }
>
> start <- c(1,1)
> zz<-optim(start,fn=z,hessian=T)
> zz
>
> Thanks in anticipation
>
> Chris Guure
> Researcher
> Institute for Mathematical Research
> UPM
