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.