# Bayes'sche Statistik: Gauss'sches Modell. A priori Verteilung fuer
# theta (Erwartungswert) ist N(xi,ka2), A-priori Verteilung fuer 1/sigma^2
# (sogenannte Praezision) ist Gamma(gam,lam), und die beiden Parameter sind
# a priori unabhängig. 
x <- rnorm(4) #Erzeugung der Daten
xb <- mean(x)
n <- length(x)
xs <- var(x)*(n-1)/n

#Wahl der Hyperparameter
#1. Nicht-informativ
xi <- 0
ka2 <- 100*xs
gam <- 1
lam <- 1
#2. Informativ fuer theta, mit den Daten uebereinstimmend
xi <- xb
ka2 <- xs/n
gam <- 1
lam <- 1
#3. Informativ fuer theta, im Konflikt mit den Daten
xi <- xb+3*sqrt(xs/n)
ka2 <- xs/n
gam <- 1
lam <- 1 


f.post.contour <- function(xi,ka2,gam,lam,n,xb,xs)
{
  ## Purpose: Plotting contours of the posterior. The main task is to define
  ##          a suitable domain where the posterior is concentrated 
  ## -------------------------------------------------------------------------
  ## Arguments: xi,ka2,al,lam=hyperparameters, xb,xs=mle
  ##            n=sample size (n=0 gives the prior !)
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  1 Nov 1999, 18:26
  s2 <- lam/gam #gam/lam ist das a priori Mittel von 1/sigma^2
  a <- s2/(s2+n*ka2)
  mu <- a*xi+(1-a)*xb #a posteriori Mittel von theta falls sigma^2=s2
  th <- seq(mu-3*sqrt(a*ka2),mu+3*sqrt(a*ka2),length=32)
                                        #Bereich fuer theta festlegen
  lam <- lam+0.5*n*xs #Skalenparameter der a posteriori Verteilung von
                      #sigma^2 fuer theta=xb
  gam <- gam+0.5*n #Formparameter dito
  si <- sqrt(lam)*seq(1/sqrt(qgamma(0.99,gam)),1/sqrt(qgamma(0.01,gam)),
                      length=32)#Bereich fuer sigma festlegen
  z <- 0.5*outer((th-xi)^2,rep(1,32))/ka2+(2*gam+1)*outer(rep(1,32),log(si))+
    outer(lam+0.5*n*(th-xb)^2,si^2,FUN="/")
                                        #log posterior auf einem Gitter
  contour(th,si,-z,nlevels=20,xlab="theta",ylab="sigma")
  }


#Simulation von der a posteriori Verteilung mit Metropolis
#(log Skala für die Standardabweichung)

#Allgemeine Funktion fuer Metropolis random walk
f.metrop <- function(x0,scale,nrep,r.target, ...)
{
  ## Purpose: Metropolis algorithm with Gaussian random walk proposals
  ## -------------------------------------------------------------------------
  ## Arguments: x0 = (vector of) starting values, scale = vector of
  ##            standard deviations for the random walk, nrep=number of
  ##            replicates, r.target = ratio of target densities
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 22 Jan 2002, 19:21
  d <- length(x0)
  x <- matrix(rnorm(d*nrep),nrow=nrep,ncol=d)
  y <- matrix(0,nrow=nrep-1,ncol=d)
  x <- x%*%diag(scale,nrow=d)
  x[1,] <- x0
  nrej <- 0
  for (i in (2:nrep)) {
    x[i,] <- x[i-1,]+x[i,]
    y[i-1,] <- x[i,]
    if (runif(1) > r.target(x[i,],x[i-1,],...)) {
      x[i,] <- x[i-1,]
      nrej <- nrej+1
    }
  }
  list(sample=x,acc=1-nrej/nrep,prop=y)
}

#Ratio der a posteriori-Verteilungen
r.post <- function(x1,x2,xi,ka2,gam,lam,n,xb,xs)
{
  ## Purpose: Quotient der a posteriori Dichten an den Stellen x1 und x2
  ## ----------------------------------------------------------------------
  ## Arguments: x=(mu,log(sigma)) = Parameter der Normal-Verteilung
  ##            xi,ka2,gam,lam = Parameter der a priori Verteilung
  ##            n=Anzahl Daten, xb=Mittel, xs=Stichprobenvarianz
  ## ----------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 24 Jan 2004, 17:25
  quot <- 0.5*((x2[1]-xi)^2-(x1[1]-xi)^2)/ka2 + (2*gam+n)*(x2[2]-x1[2]) +
    (lam+0.5*n*(xs+(x2[1]-xb)^2))*exp(-2*x2[2]) -
    (lam+0.5*n*(xs+(x1[1]-xb)^2))*exp(-2*x1[2])
  exp(quot)
}


f.post.contour(xi,ka2,gam,lam,n,xb,xs)
points(xb,sqrt(xs),pch=15)
title(main="log a posteriori Dichte")
#Wenige Iterationen fuer didaktische Zwecke
metrop <- f.metrop(x0=c(xb,0.5*log(xs)),scale=c(0.6,0.3),nrep=21,
                   r.target=r.post,xi,ka2,gam,lam,n,xb,xs)
sample <- metrop$sample
sample[,2] <- exp(sample[,2])
prop <- metrop$prop
prop[,2] <- exp(prop[,2])
points(sample[,1],sample[,2],col=2)
arrows(sample[-21,1],sample[-21,2],prop[,1],prop[,2],col=4)
#Resultate von vielen Iterationen
library("ts")
metrop <- f.metrop(x0=c(xb,0.5*log(xs)),scale=c(0.6,0.3),nrep=1000,
                   r.target=r.post,xi,ka2,gam,lam,n,xb,xs)
metrop$acc  #Akzeptierungsrate
param <- metrop$sample
points(param[,1]+0.04*(runif(1000)-0.5),
       exp(param[,2])+0.04*(runif(1000)-0.5),cex=0.5)
#Punkte mit Ueberlappung sichtbar gemacht
par(mfrow=c(2,1))
plot(param[1:500,1],type="l")
plot(param[1:500,2],type="l")
acf(param[,1])
acf(param[,2])