poisson.lognormal.mcmc<-function(sumx,n,prior=list("a"=1,"b"=1,"tau"=1000),I=1000, do.tuning=FALSE, tuning=50, do.burnin=FALSE, burnin=100){ mu<-rep(NA,I) lambda<-rep(NA,I) mu[1]<-1 lambda[1]<-sumx/n tau=prior[["tau"]] a=prior[["a"]] b=prior[["b"]] #Berechnung der log Full Conditional von mu log.fc.mu <- function(mu,lambda,sumx,tau) { return(mu*log(lambda)-mu^2-mu^2/tau) } log.fc.lambda <- function(lambda,mu,sumx,n){ return((sumx+mu-1)*log(lambda)-n*lambda-0.5*log(lambda)^2) } rw.sd <- .1 akzeptanz <- 0 #Ziehen der Zufallszahlen mittel MH-Algorithmus i=1 while (i<=I) { i<-i+1 mustern <- rnorm(1,mu[i-1],rw.sd) sumx<-sum(data$n) logalpha<-log.fc.mu(mustern,lambda[i-1],sumx,tau)-log.fc.mu(mu[i-1],lambda[i-1],sumx,tau) alpha<-exp(logalpha) ifelse(runif(1).5){ rw.sd <-rw.sd * ak.rate * 2 akzeptanz = 0 i = 1 lambda[1] <- lambda[tuning] mu[1] <- mu[tuning] } if (ak.rate<.3){ rw.sd <-rw.sd * ak.rate akzeptanz = 0 i = 1 lambda[1] <- lambda[tuning] mu[1] <- mu[tuning] } cat(paste("Akzeptanzrate",ak.rate,"\n")) } if(i==burnin)if (do.burnin) { i=1 do.burnin=FALSE } } lambda<-lambda[-1] mu<-mu[-1] return(coda::mcmc(cbind(mu,lambda))) }