##########################################
#fiche de TD n°10 - méthode de Monte Carlo
#L. Carraro - O. Roustant
#Telecom Saint-Etienne - decembre 2010
##########################################

######################
#1° - methode du rejet
######################

#on ajoute comme parametre le cote du carré dans lequel on tire les v.a. U et V
MC1<-function(n,cote=1)
{
	X<-vector(length=n)		#contient les simulations de la v.a. X
	M<-vector(length=n+1)	#contient les estimations de pi au fil des simulations
	S<-vector(length=n+1)	#contient les sommes S au numérateur de M
	M[1]<-0
	S[1]<-0
#boucle de Monte Carlo
	for (i in 1:n)	
	{
#simulation de la v.a. X à partir des v.a. U et V uniformes dans le carre [-cote,cote]*[-cote,cote]
		U<-runif(1,min=-cote,max=cote)
		V<-runif(1,min=-cote,max=cote)
		if (U^2+V^2 <= 1) X[i]<-1 else X[i]<-0	
		S[i+1]<-S[i]+X[i]					#on met à jour la somme S
		M[i+1]<-(2*cote)^2*S[i+1]/i		#on met à jour M
	}
#on dessine le graphique des simulations au cours du temps (attention M[i] est le résultat pour la simulation i+1)
	plot(M[2:n+1],pch=19,cex=0.5,col="blue",xlab="nb simulations",ylab="valeur estimée",main=paste("estimation de pi par rejet avec un coté du carré égal à",eval(cote)))
	abline(a=pi,b=0)
#on retourne toutes les valeurs prises par X et M pour des indices entre 1 et n
	result<-rbind(X,M[2:(n+1)])
}



#####################
#2° - seconde methode
#####################

MC2<-function(n)
{
	Y<-vector(length=n)		#contient les simulations de la v.a. Y
	M<-vector(length=n+1)	#contient les estimations de pi au fil des simulations
	S<-vector(length=n+1)	#contient les sommes S au numérateur de M
	M[1]<-0
	S[1]<-0
#boucle de Monte Carlo
	for (i in 1:n)	
	{
#simulatin de la v.a. Y à partir de la v.a. U uniforme dans [0,1]
		U<-runif(1)
		Y[i]<-sqrt(1-U^2)	
		S[i+1]<-S[i]+Y[i]		#on met à jour la somme S
		M[i+1]<-4*S[i+1]/i		#on met à jour M
	}
#on dessine le graphique des simulations au cours du temps (attention M[i] est le résultat pour la simulation i+1)
	plot(M[2:n+1],pch=19,cex=0.5,col="blue",xlab="nb simulations",ylab="valeur estimée",main="estimation de pi par la seconde méthode")
	abline(a=pi,b=0)
#on retourne toutes les valeurs prises par M pour des indices entre 1 et n
	result<-rbind(Y,M[2:(n+1)])
}

########################################
#comparaison des deux premieres methodes
########################################
n<-1000
M1<-MC1(n)
M2<-MC2(n)
plot(M1[2,],pch=19,cex=0.5,col="blue",xlab="nb simulations",ylab="valeur estimée",main="comparaison des méthodes 1 et 2",ylim=c(2.5,3.5))
abline(a=pi,b=0)
points(M2[2,],pch=19,cex=0.5,col="red")
legend("topright",legend=c("méthode 1","méthode 2"),pch=c(19,19),col=c("blue","red"))

#######################
#3° - troisieme methode
#######################

#fonction devant etre modifiee

MC3<-function(n)
{
	Z<-vector(length=n)		#contient les simulations de la v.a. Z
	M<-vector(length=n+1)	#contient les estimations de pi au fil des simulations
	S<-vector(length=n+1)	#contient les sommes S au numérateur de M
	M[1]<-0
	S[1]<-0
#boucle de Monte Carlo
	for (i in 1:n)	
	{
#simulatin de la v.a. Z à partir de la v.a. U uniforme dans [0,1]
		U<-runif(1)
#		Z[i]<-########??????????###########	
		S[i+1]<-S[i]+Z[i]			#on met à jour la somme S
#		M[i+1]<-#############????????###########*S[i+1]/i		#on met à jour M
	}
#on dessine le graphique des simulations au cours du temps (attention M[i] est le résultat pour la simulation i+1)
	plot(M[2:n+1],pch=19,cex=0.5,col="blue",xlab="nb simulations",ylab="valeur estimée",main="estimation de pi par la seconde méthode")
	abline(a=pi,b=0)
#on retourne toutes les valeurs prises par Z et M pour des indices entre 1 et n
	result<-rbind(Z,M[2:(n+1)])
}