library(mvtnorm)

# Gibbs sampler, bivariate Normalverteilung, Cov=Cor=bekannt, non-informative
# priori \propto const
varx1 <- 1
varx2 <- 1

# Korrelation höher -> Sampler braucht länger
corrx1x2 <- 0.90
covmat <- matrix( nrow=2, ncol=2, data=c(varx1, corrx1x2, corrx1x2, varx2))
print(covmat)

# Daten erzeugen
n <- 100
x <- rmvnorm(n=n, mean=c(3,10), sigma=covmat )

# Jetzt mu_1 und mu_2 mittels Gibbs-Sampler schätzen
# summary measures
x1bar <- mean( x[,1] )
x2bar <- mean( x[,2] )

numBurnIn <- 100
numSamples <- 300
mu <- matrix( nrow=(numBurnIn+numSamples), ncol=2, data=0 )

# Startwerte
mu[1,2] <- -7

# Full-conditional mu_1 | mu_2, daten, START
mu[1,1] <- rnorm( n=1, mean=x1bar + corrx1x2 * (mu[1,2]-x2bar), sd=sqrt( (1-corrx1x2^2)/n) )

# jetzt gehts erst richtig los ...

for ( i in 2:(numBurnIn+numSamples) ) {
  # Full-conditional mu_2 | mu_1, daten
  mu[i,2] <- rnorm( n=1, mean=x2bar + corrx1x2 * (mu[i-1,1]-x1bar), sd=sqrt( (1-corrx1x2^2)/n))
  # Full-conditional mu_1 | mu_2, daten, START
  mu[i,1] <- rnorm( n=1, mean=x1bar + corrx1x2 * (mu[i,2]-x2bar), sd=sqrt( (1-corrx1x2^2)/n))
  matplot( mu[1:i, ], type="l")
}

print("Geschätzte a posteriori-Erwartungswerte von mu_1 und mu_2 (ohne Burnin):")
print( colMeans(mu[(numBurnIn+1):(numBurnIn+numSamples),] )  )
print("Im Vergleich dazu ML-Schätzer:")
print(colMeans(x))