#Esimerkki Metropolis-algoritmin käytöstä

metro<-function(f,start,simlkm=1000, Sigma=0.1*diag(length(start)), ...)
# Funktio Metropolis-algoritmin ajamiseksi
# Funktio antaa tuloksena matriisin, jonka riveillä ovat generoidut satunnaisvektorit
# Parametrit: 
# f: simuloitavan jakauman tiheysfunktion logaritmi
# start: alkuarvot
# simlkm: simulointien lkm
# Sigma: hyppäysjakauman kovarianssimatriisi
# ...: Funktion f parametrit 
{
  k<-length(start)
  C<-t(chol(Sigma))
  zmat<-matrix(0,nr=simlkm,nc=k)
  zmat[1,]<-start  

  for(i in 2:simlkm)
  {
      zold <- zmat[i-1,]
      znew <-C%*%rnorm(k)+zold 
      prob<-min(1,exp(f(znew, ...)-f(zold, ...)))
      zmat[i,] <- if(rbinom(1,1,prob)==1) znew else zold
  }
  zmat
}


#Esimerkki: generoidaan 100 havaintoa Weibull(0.3,10)-jakaumasta ja pyritään 
#estimoimaan parametrit

data<-rweibull(100,0.3,10)
plot(density(data))

#Logaritmoitu posteriorijakauman tiheysfunktio:

lposterior <- function(param,y)
{
  delta<-param[1]
  beta <-param[2]
  n <- length(y)
 if(delta <=0 || beta<=0)
   return(-100000)
 (n-1)*log(delta)-(n*delta-1)*log(beta)+(delta-1)*sum(log(y))-sum(y^delta)/beta^delta
}

# Varsinainen simulointi ja suppenemisen tarkastelu 

sim <-metro(lposterior,simlkm=10000,start=c(1,1),y=data,Sigma=diag(c(0.01,1)))
colnames(sim)<-c("delta","beta")
plot(ts(sim))
Sigma<-2.4^2/2*var(sim[1000:10000,])
sim[,] <-metro(lposterior,simlkm=10000,start=sim[10000,],y=data,Sigma)
plot(ts(sim))
library(MCMCpack)
cumuplot(mcmc(sim))