mu.true<-0.5
nu2.true<-0.5
N<-50

loglambda<-mu.true+rnorm(N,0,nu2.true)
lambda<-exp(loglambda)
x<-rpois(N,lambda)

print(mean(x))
print(var(x))

plot(x)
plot(density(x))

mu<-log(mean(x))
eps<-rep(0,N)
nu2<-1

rw.var<-1
sumx<-sum(x)
tau2<-1000
mu.save<-nu2.save<-eps.save<-NULL
mu.save2<-nu2.save2<-eps.save2<-NULL

for (i in 1:100)
{
  mu.star<-rnorm(1,mu,rw.var)  
  logalpha<-sumx*mu.star-sum(exp(mu.star+eps))-mu.star*mu.star/tau2/2
  logalpha<-logalpha-sumx*mu+sum(exp(mu+eps))+mu*mu/tau2/2
  alpha<-exp(logalpha)
  if (runif(1)<alpha)mu<-mu.star

  for (j in 1:N)
  {
    eps.star<-rnorm(1,eps[j],rw.var)  
    logalpha<-eps.star*x[j]-exp(mu+eps.star)-eps.star*eps.star/nu2/2
    logalpha<-logalpha-eps[j]*x[j]+exp(mu+eps[j])+eps[j]*eps[j]/nu2/2
    alpha<-exp(logalpha)
    if (runif(1)<alpha)eps[j]<-eps.star
  }

  a<-0.1+N/2
  b<-0.1+0.5*sum(eps*eps)
  nu2<-1/rgamma(1,a,b)
}

mu.start<-mu
eps.start<-eps
nu2.start<-nu2

for (i in 1:1000)
{
  mu.star<-rnorm(1,mu,rw.var)  
  logalpha<-sumx*mu.star-sum(exp(mu.star+eps))-mu.star*mu.star/tau2/2
  logalpha<-logalpha-sumx*mu+sum(exp(mu+eps))+mu*mu/tau2/2
  alpha<-exp(logalpha)
  if (runif(1)<alpha)mu<-mu.star
  mu.save<-c(mu.save,mu)
  
  for (j in 1:N)
  {
    eps.star<-rnorm(1,eps[j],rw.var)  
    logalpha<-eps.star*x[j]-exp(mu+eps.star)-eps.star*eps.star/nu2/2
    logalpha<-logalpha-eps[j]*x[j]+exp(mu+eps[j])+eps[j]*eps[j]/nu2/2
    alpha<-exp(logalpha)
    if (runif(1)<alpha)eps[j]<-eps.star
  }
  eps.save<-rbind(eps.save,eps)
  
  a<-1+N/2
  b<-1+0.5*sum(eps*eps)
  nu2<-1/rgamma(1,a,b)
  nu2.save<-c(nu2.save,nu2)
}
plot(mu.save,type="l")
plot(apply(eps.save,1,mean),type="l")
plot(eps.save[,23],type="l")
plot(nu2.save,type="l")

mu<-mu.start
eps<-eps.start
nu2<-nu2.start

for (i in 1:1000)
{
  mu.star<-rnorm(1,mu,rw.var)  
  logalpha<-sumx*mu.star-sum(exp(mu.star+eps))-mu.star*mu.star/tau2/2
  logalpha<-logalpha-sumx*mu+sum(exp(mu+eps))+mu*mu/tau2/2
  alpha<-exp(logalpha)
  if (runif(1)<alpha)mu<-mu.star
  mu.save2<-c(mu.save2,mu)
  
  for (j in 1:N)
  {
    eps.star<-rnorm(1,eps[j],rw.var)  
    logalpha<-eps.star*x[j]-exp(mu+eps.star)-eps.star*eps.star/nu2/2
    logalpha<-logalpha-eps[j]*x[j]+exp(mu+eps[j])+eps[j]*eps[j]/nu2/2
    alpha<-exp(logalpha)
    if (runif(1)<alpha)eps[j]<-eps.star
  }
  eps.save2<-rbind(eps.save2,eps)
  
  a<-0.1+N/2
  b<-0.1+0.5*sum(eps*eps)
  nu2<-1/rgamma(1,a,b)
  nu2.save2<-c(nu2.save2,nu2)
}

coda::gelman.diag(list(coda::mcmc(nu2.save),coda::mcmc(nu2.save2)))
coda::gelman.diag(list(coda::mcmc(eps.save),coda::mcmc(eps.save2)))
coda::gelman.diag(list(coda::mcmc(mu.save),coda::mcmc(mu.save2)))

plot(density(c(mu.save,mu.save2)))
plot(density(c(nu2.save,nu2.save)))
plot(density(c(eps.save[,1],eps.save2[,1])))
for (i in 2:50)lines(density(c(eps.save[,i],eps.save2[,i])))     

image(MASS::kde2d(eps.save[,1],eps.save[,2]))