#################################################################
# Uebung zu Schaetzen und Testen I im WiSe 2015/2016
# R-Loesung zur Uebung 2 am 10.12.2016
# Autor: Ludwig Bothmann; Ueberarbeitet: David Ruegamer
#################################################################

################################
# Zufallszahlen generieren
################################

# Wahrer Erwartungswert und Varianz
mu <- 5
sigma <- 1

# seed setzen und Ergebnisse reproduzierbar zu machen
set.seed(1054)

# 10 Zufallszahlen aus einer NV ziehen 
x <- rnorm(10, mu, sigma)

# Auf 2 Nachkommastellen runden (wie auf Uebungsblatt)
round(x, 2)

##################
# Aufgabe 18.d)
##################

nu <- 0
tau2 <- 10000
sigma2 <- 1

# Funktion, die fuer gegebene Stichprobe den Posteriori-EW berechnet
# vgl Teilaufgabe (b) fuer Berechnung des Posteriori-EW
mu.post <- function(x, sigma2=1, nu=0, tau2=10000){
	n <- length(x)
	return((n/sigma2 + 1/tau2)^(-1) * (n*mean(x)/sigma2 + nu/tau2))
}

# Da die Normalverteilung symmetrisch ist, sind Posteriori-EW, -Modus und
#	-Median identisch:

mu.x <- mu.post(x)
mu.x
# Zum Vergleich: ML-Schaetzer
mean(x)

# Funktion, die fuer gegebene Stichprobe die Posteriori-Varianz berechnet
# vgl Teilaufgabe (b) fuer Berechnung der Posteriori-Varianz
sigma2.post <- function(x, sigma2=1, tau2=10000){
	n <- length(x)
	return((n/sigma2 + 1/tau2)^(-1))
}

sigma2.x <- sigma2.post(x)
sigma2.x

# Gleichendiges Kredibilitaetsintervall: 
# 2.5% und 97.5% - Quantil der Posteriori-Verteilung berechnen, fuer alpha=5%
alpha <- 0.05

KI <- qnorm(c(alpha/2, 1-alpha/2), mean=mu.x, sd=sqrt(sigma2.x))
KI

# HPD-Intervall: Identisch, da Normalverteilung symmetrisch und
#	eingipfelig

##################
# Aufgabe 18.e)
##################

x.100 <- rep(x, 10)
x.100

# Punktschätzer: Bleibt (nahezu) gleich
mu.x100 <- mu.post(x.100)
mu.x100
mu.x

# (geschaetztes sigma ist kleiner fuer n=100, als fuer n=10)
sigma2.x100 <- sigma2.post(x.100)
sqrt(sigma2.x100)
sqrt(sigma2.x)

# Intervallschaetzer: Kleiner, da die Daten mehr Information enthalten
#  => Unsicherheit ist geringer
KI.100 <- qnorm(c(alpha/2, 1-alpha/2), mean=mu.x100, sd=sqrt(sigma2.x100))
KI.100
KI

##################
# Aufgabe 18.f)
##################


# Gitter von Werten von mu, fuer die die Plots erstellt werden sollen
mu.grid <- seq(-6, 8, length=1000)

# Funktion zur Berechnung der Likelihood
lik <- function(mu, x, sigma2){
	
	loglik <- sum( -log(sqrt(2*pi*sigma2)) - (1/(2*sigma2)) * (x-mu)^2)
	
	return(exp(loglik))
}


# Plot von Priori-Dichte, Likelihood und Posteriori-Dichte:
#	Funktion, an welche die Priori-Parameter uebergeben werden koennen
plot.func <- function(nu=0, tau2=10000, sigma2=1, alpha=0.05, plot.KI=TRUE){
	
	mu.x <- mu.post(x, nu=nu, tau2=tau2)
	sigma2.x <- sigma2.post(x, tau2=tau2)
	KI <- qnorm(c(alpha/2, 1-alpha/2), mean=mu.x, sd=sqrt(sigma2.x))
	
	mu.x100 <- mu.post(x.100, nu=nu, tau2=tau2)
	sigma2.x100 <- sigma2.post(x.100, tau2=tau2)
	KI.100 <- qnorm(c(alpha/2, 1-alpha/2), mean=mu.x100, sd=sqrt(sigma2.x100))
	
	
	# Plot von Priori-Dichte, Likelihood und Posteriori-Dichte

	# Fuer Stichprobe x
	
	# Priori-Dichte
	dens.mu <- dnorm(mu.grid, mean=nu, sd=sqrt(tau2))
	plot(mu.grid, dens.mu, ylim=c(0, max(dens.mu)),
		  type="l", xlab=expression(mu), ylab=expression(f(mu)), 
		  main=paste("Priori-Dichte mit nu = ", nu, ", tau = ", sqrt(tau2), sep=""))
	
	# Likelihood fuer Stichprobe x (= Dichte einer Multivariaten NV)
	plot(mu.grid, sapply(mu.grid, FUN=lik, x=x, sigma2=sigma2), 
		  type="l", xlab=expression(mu), ylab=expression(L(mu,x)), 
		  main="Likelihood (Stichprobe: x)")
	
	# Posteriori-Dichte 
	plot(mu.grid, dnorm(mu.grid, mean=mu.x, sd=sqrt(sigma2.x)), 
		  type="l", xlab=expression(mu), ylab=expression(f(mu,x)), 
		  main="Posteriori-Dichte")
	
	
	if(plot.KI){
		abline(v=KI, col="blue", lty=2)
		legend("topleft", paste(1-alpha, " - KI",sep=""), col="blue", lty = 2)
	}
	
	# Fuer Stichprobe x.100
	# Priori-Dichte
	dens.mu <- dnorm(mu.grid, mean=nu, sd=sqrt(tau2))
	plot(mu.grid, dens.mu, ylim=c(0, max(dens.mu)),
		  type="l", xlab=expression(mu), ylab=expression(f(mu)), 
		  main=paste("Priori-Dichte mit nu = ", nu, ", tau = ", sqrt(tau2), sep=""))
	
	# Likelihood für Stichprobe x (= Dichte einer Multivariaten NV)
	plot(mu.grid, sapply(mu.grid, FUN=lik, x=x.100, sigma2=sigma2), 
		  type="l", xlab=expression(mu), ylab=expression(L(mu,x)), 
		  main="Likelihood (Stichprobe: x.100)")
	
	# Posteriori-Dichte 
	plot(mu.grid, dnorm(mu.grid, mean=mu.x100, sd=sqrt(sigma2.x100)), 
		  type="l", xlab=expression(mu), ylab=expression(f(mu,x)), 
		  main="Posteriori-Dichte")
	
	if(plot.KI){
		abline(v=KI.100, col="blue", lty=2)
		legend("topleft", paste(1-alpha, " - KI",sep=""), col="blue", lty=2)
	}
}

# Plots als pdf abspeichern

pdf("Plots_Aufgabe18.pdf", width=12, height=24)
par(mfrow=c(6,3), mar=c(5,5,4,2))
	# Parameter wie auf dem Uebungsblatt
	plot.func()
	
	# Kleinere Priori-Varianz
	plot.func(nu=0, tau2=1, alpha=0.05, plot.KI=TRUE)
	# => Auswirkung auf Posteriori bei x groesser als bei x.100
	# => Klar: x.100 hat ja auch mehr Information, 
	#		die Likelihood ist "staerker"
	
	# Stark informative Priori, die nicht zu den Daten passt
	plot.func(nu=-5, tau2=0.01, alpha=0.05, plot.KI=TRUE)
	# => Bei x fallen die Daten kaum noch ins Gewicht, bei x.100 sind
	#		Priori und Likelihood in etwa gleich stark

dev.off()