# simulation de la distribution normale
# a l'aide de l'algorithme de Metropolis-Hasting 
# 
# parametres:
# n    -- nombre d'iterations de l'algorithme
# step -- longeur du saut maximal pour une marche aleatoire sous-jacente
# x0   -- position initiale (0 par defaut)
#
# retour: un vecteur de longueur n, qui contient la trajectoire d'une
#         chaine de Markov, generee par l'algorithme M-H

mhnorm <- function(n, step=.1, x0=0) {
	x <- vector("numeric", n)
	d <- runif(n, -step, step)
	u  <- runif(n)

	# une fonction proportionnelle a la densite normale
	# (nous n'avons pas besoin de la constante 1/sqrt(2*pi))
	f <- function(x) { exp(-x*x/2) }

	for(i in 1:n) {
		x[i] <- x0;

		# valeur propose
		x1 <- x0 + d[i]
		v0 = f(x0)
		v1 = f(x1)

		# condition d'acceptation
		if(v0*u[i] < v1) x0 <- x1;
	}

	return(x);
}

# exemple 1:
# histogramme de valeurs de mhnorm(n)
# tracee contre la densite de la loi normale

# essayez les commandes suivantes
# et expliquez le resultat:
# 
# exemple1(n=1000, step=.1)
# exemple1(n=1000, step=2)
# exemple1(n=1000, step=3)
# exemple1(n=10000, step=3)
# exemple1(n=10000, step=2)
# exemple1(n=10000, step=.1)

exemple1 <- function(n=10000, step=.1) {
	x <- mhnorm(n, step)
	hist(x, breaks=100, freq=FALSE)
	curve(dnorm, col='red', add=TRUE)
}

# exemple 2:
# deux suites independantes de mhnorm(n)
# tracees comme un chemin en deux dimensions
# (permet d'illustrer la dependence de valuers)
		
# essayez les commandes suivantes
# et expliquez le resultat:
# 
# exemple2(1000, .1)
# exemple2(1000, 3)
# exemple2(1000, 10)
# exemple2(10000, .01)

		
exemple2 <- function(n=1000, step=.1) {
	x <- mhnorm(n, step)
	y <- mhnorm(n, step)
	z <- complex(real=x, imag=y)
	plot(z, type='l')
}