#################################################################
# Uebung zu Schaetzen und Testen 1 im WiSe 15/16 
# R-Loesung zur Uebung 12, Aufgabe 1
# Funktionen
# Autor: Ludwig Bothmann; Aenderungen: David Ruegamer
#################################################################

# Pakete laden
library(mvtnorm)
library(MCMCpack)
library(coda)

################################################
# Metropolis-Algorithmus für das Logit-Modell
################################################

logit.sample <- function(D=1000,
								 y,
								 X,
								 B.0,
								 beta.0,
								 beta.start,
								 B.prop, 
								 do.for=TRUE){
	
	# Anzahl Kovariablen
	p <- ncol(X)
	
	# Matrix, in die spaeter die Samples gespeichert werden
	beta.sample <- matrix(ncol=p, nrow=D+1)
	beta.sample[1,] <- beta.start
	
	# Der erste "alte" Wert ist der Startwert
	beta.alt <- beta.start
	
	# Eigentlicher Algorithmus
	if(do.for){
		
		# For-Schleife: Richtiger MH-Algorithmus
		
		# Anzahl der akzeptierten Vorschläge
		d <- 0
		
		for(i in 1:D){
		
			# 1.) Ziehung aus Vorschlagsdichte
			beta.neu <- rmvnorm(n=1, mean=beta.alt, sigma=B.prop)
			
			# 2.) Akzeptanzrate bestimmen
			akzeptanzrate <- alpha.rate(beta.alt, 
												 beta.neu, 
												 y=y, 
												 X=X, 
												 beta.0=beta.0, 
												 B.0=B.0)
			
			# 3.) Zufaellig akzeptieren mit W' akzeptanz
			akzeptieren <- rbinom(n=1,size=1, prob=akzeptanzrate)
			
			# Falls akzeptiert wird
			if(akzeptieren){
				
				# Vorschlag abspeichern
				beta.sample[i+1,] <- beta.neu
				
				# Der akzeptierte Vorschlag ist ab jetzt der "alte" Wert
				beta.alt <- beta.neu			
								
				# Anzahl der akzeptierten Werte
				d <- d+1	
				
			}else{
				# Falls nicht akzeptiert wird: Vorherigen Wert noch mal 
				#	abspeichern
				
				# Alten Vorschlag abspeichern
				beta.sample[i+1,] <- beta.alt
			}
		
		}
		out <- list(beta.sample = beta.sample, d=d)
		
	}else{
	
		# While-Schleife: In der Uebung vorgestellter, 
		#	aber falscher MH-Algorithmus
		
		d <- 1
		i <- 0
		
		# Vektor, in den spaeter gespeichert wird, welche Vorschlaege akzeptiert
		#	wurden
		d.vec <- vector(length=d)
		
		while(d<=D){
			
			# Wie oft musste gezogen werden?
			i <- i+1
			
			# 1.) Ziehung aus Vorschlagsdichte
			beta.neu <- rmvnorm(n=1, mean=beta.alt, sigma=B.prop)
			
			# 2.) Akzeptanzrate bestimmen
			akzeptanzrate <- alpha.rate(beta.alt, 
										  beta.neu, 
										  y=y, 
										  X=X, 
										  beta.0=beta.0, 
										  B.0=B.0)
			
			# 3.) Zufaellig akzeptieren mit W' akzeptanz
			akzeptieren <- rbinom(n=1,size=1, prob=akzeptanzrate)
			
			# Falls akzeptiert wird: Abspeichern und d um 1 erhoehen
			if(akzeptieren){
				
				# Vorschlag abspeichern
				beta.sample[d+1,] <- beta.neu
				
				# Vorschlag ist ab jetzt der "alte" Wert
				beta.alt <- beta.neu			
				
				# Zeitpunkt der Ziehung
				d.vec[d] <- i
				
				# Naechste Iteration
				d <- d+1	
			}
			# sonst: Wieder von oben		
		}
		
		out <- list(beta.sample = beta.sample,
						d.vec = d.vec)
	}


	
	return(out)
}

################################################
# Hilfsfunktionen
################################################

###########
# 1.) Vorschlagsdichte: N(beta.kq, B.kq)
# => Braucht nicht implementiert zu werden: 
#		einfach dmvnorm und rmvnorm aus Paket mvtnorm benutzen
###########

###########
# 2.) Responsefunktion
###########

h <- function(x){ 
	h <- exp(x)/(1+exp(x)) 
	return(h)
}


###########
# 3.) Posteriori-Dichte (bis auf Proportionalitaet)
###########

post <- function(beta, y, X, beta.0, B.0){
	
	# Likelihood berechnen
	lik <- 1
	
	for(i in 1:length(y)){
		lik <- lik * (h(t(X[i,]) %*% beta))^y[i] * (1-h(t(X[i,]) %*% beta))^(1-y[i]) 
	}
	
	# Priori berechnen
	prio <- exp(-1/2 * t(beta - beta.0) %*% solve(B.0) %*% (beta - beta.0))
	
	# Posteriori berechnen
	post <- lik * prio
	
	return(post)
	
}


###########
# 4.)  Akzeptanzrate
###########

alpha.rate <- function(beta.alt, 
						beta.neu, 
# 						B.prop, # nicht nötig
						y,
						X,
						beta.0,
						B.0
						){
	
	# 	ratio.vorschlag <- dmvnorm(x=beta.alt, mean=beta.neu, sigma=B.prop) / 
	# 		dmvnorm(x=beta.neu, mean=beta.alt, sigma=B.prop)	
	# => Immer = 1, da symmetrisch
	
	ratio.post <- post(t(beta.neu), y, X, beta.0, B.0) / 
		post(t(beta.alt), y, X, beta.0, B.0)
	
	# 	a <- min(1, ratio.vorschlag * ratio.post)
	a <- min(1, ratio.post)
	
}