I am facing a similar problem
The code which i am pasting below has been taken from the book "Bayesian
Computation with R" by Jim Albert.
laplace=function (logpost, mode, par)
{
options(warn=-1)
fit=optim(mode, logpost, gr = NULL, par, hessian=TRUE,
control=list(fnscale=-1))
options(warn=0)
mode=fit$par
h=-solve(fit$hessian)
p=length(mode)
int = p/2 * log(2 * pi) + 0.5 * log(det(h)) +
logpost(mode, par)
stuff = list(mode = mode, var = h, int = int,
converge=fit$convergence==0)
return(stuff)
}
A<-read.table(file="chemotherapy.txt",head=TRUE)
A
library(survival)
Model_01=survreg(Surv(A$time,A$status)~factor(A$treat)+A$age,dist="weibull")
weibullregpost=function (theta, A)
{
logf=function(t,c,x,sigma,mu,beta)
{
z=(log(t)-mu-x%*%beta)/sigma
f=1/sigma*exp(z-exp(z))
S=exp(-exp(z))
c*log(f)+(1-c)*log(S)
}
k = dim(A)[2]
p = k - 2
t = A[, 1]
c = A[, 2]
X = A[, 3:k]
sigma = exp(theta[1])
mu = theta[2]
beta = array(theta[3:k], c(p,1))
return(sum(logf(t,c,X,sigma,mu,beta)))
}
start=array(c(-Model_01$scale, Model_01$coefficients),c(1,4))
d=cbind(A$time,A$status,A$treat-1,A$age)
fit=laplace(weibullregpost, start, d)
fit
proposal=list(var=fit$var,scale=1.5)
rwmetrop=function (logpost, proposal, start, m, par)
{
pb = length(start)
Mpar = array(0, c(m, pb))
b = matrix(t(start))
lb = logpost(start, par)
a = chol(proposal$var)
scale = proposal$scale
accept = 0
for (i in 1:m) {
bc = b + scale * t(a) %*% array(rnorm(pb), c(pb, 1))
lbc = logpost(t(bc), par)
prob = exp(lbc - lb)
if (is.na(prob) == FALSE) {
if (runif(1) < prob) {
lb = lbc
b = bc
accept = accept + 1
}
}
Mpar[i, ] = b
}
accept = accept/m
stuff = list(par = Mpar, accept = accept)
return(stuff)
}
bayesfit=rwmetrop(weibullregpost,proposal,fit$mode,10000,d)
bayesfit$accept
par(mfrow=c(2,2))
sigma=exp(bayesfit$par[,1])
mu=bayesfit$par[,2]
beta1=bayesfit$par[,3]
beta2=bayesfit$par[,4]
hist(beta1,xlab="treatment")
hist(beta2,xlab="age",main="")
hist(sigma,xlab="sigma",main="")
The code is running fine for the data given in the book, but when i am using
it on my data which is quite large compared to book's data, i am gettin same
error.
Then i tried changing the parameters and tested a lot of them but the error
remains.
The other error which i got was
Error in chol.default(proposal$var) :
the leading minor of order 1 is not positive definite
In addition: Warning message:
In log(det(h)) : NaNs produced
Can anybody suggest some way to get the right parameters?
It is really very urgent and i desperately need some help!!!
Thanks a lot
sandsky wrote:
>
> I have an error for a simple optimization problem. Is there anyone knowing
> about this error?
>
> lambda1=-9
> lambda2=-6
> L<-function(a){
> s2i2f<-(exp(-lambda1*(250^a)-lambda2*(275^a-250^a))
> -exp(-lambda1*(250^a)-lambda2*(300^a-250^a)))
> logl<-log(s2i2f)
> return(-logl)}
> optim(1,L)
>
> Error in optim(1, L) : function cannot be evaluated at initial parameters
>
> Thank you in advance
>
--
View this message in context:
http://www.nabble.com/Error%3A--function-cannot-be-evaluated-at-initial-parameters-tp19565286p24846455.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
