#Feuille de TD n°4

#avec les fonctions usuelles

#1)
# rnorm(n,m,s) permet de générer n observations d'une loi normale
					# d'espérance m et d'écart-type s
# dnorm(x,m,s) permet d'évaluer la fonction de densité d'une loi normale
					#d'espérance m et d'écart-type s aux points
					#du vecteur x
# qnorm(p,m,s) permet d'évaluer la fonction quantile d'une loi normale
					#d'espérance m et d'écart-type s aux points
					#du vecteur p (chaque point etant entre 0 et 1)
# pnorm(x,m,s) permet d'évaluer la fonction de répartition d'une loi normale
					#d'espérance m et d'écart-type s aux points
					#du vecteur x

#2)
#pour simuler une
# - binomiale de paramètres N et p et réaliser n observations : rbinom(n,N,p)
# - bernoulli de paramètre p et réaliser n observations : rbinom(n,1,p)
# en effet une bernoulli de paramètre pest une binomiale de paramètre 1 et p
# - uniforme sur [a;b] et réaliser n observations : runif(n,a,b)
# - exponentielle de paramètre lambda et réaliser n observations : rexp(n,lambda)
# - poisson de paramètre lambda et réaliser n observations : rpois(n,lambda)

#3)
A=rnorm(1000,3,sqrt(10))

#4)
H=hist(A,plot=FALSE) 		#récupérer les éléments de la construction de
					#l'histogramme sans le tracer
limits=H$breaks
xmin=min(limits)			#plus petite limite de classe
xmax=max(limits)			#plus grande limite de classe
x=seq(xmin,xmax,0.01)		#vecteur sur lequel on va évaluer la densité théorique
y=dnorm(x,3,sqrt(10))		#évaluation de la densité
ymax=max(H$density,y)

hist(A,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x,y,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique

#A partir d'une variable uniforme
#1)
#par la méthode d'inversion, on sait que si F est la fonction de répartition 
#d'une variable continue X, alors F^{-1}(U), où U est une variable uniforme sur
#[0,1], a même loi que X

#ici la fonction de répartition d'une loi exponentielle de paramètre lambda est
# F(x)=1-exp(-lambda.x) si x>=0 et 0 sinon
#si on restrait F à [0,+infini[, elle admet une inverse car F est sur [0,+nfini[
#une fonction continue et strictement croissante

#on résout F(x)=y de façon à trouver x=
#la solution est x=-1/lambda.ln(1-y)
#donc la variable -1/lambda.ln(1-U) est de loi exponentielle de paramètre lambda

simuexp<-function(n,a)
{
u=runif(n)		#simulation de n observations d'une loi uniforme sur [0,1]
y=-1/a*log(1-u)	#observations d'une loi exponentielle de paramètre a
simuexp=y
}

#2)
B=simuexp(1000,2)
H=hist(B,plot=FALSE,breaks=15) 		#récupérer les éléments de la construction de
					#l'histogramme sans le tracer
limits=H$breaks
xmin=min(limits)			#plus petite limite de classe
xmax=max(limits)			#plus grande limite de classe
x=seq(xmin,xmax,0.01)		#vecteur sur lequel on va évaluer la densité théorique
y=dexp(x,2)		#évaluation de la densité
ymax=max(H$density,y)

hist(B,freq=FALSE,breaks=15,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x,y,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique



#avec 10 observations, pas suffisamment de classes pour créer un histogramme représentatif

#3)
#l'idée est de découper l'intervalle [0,1] en trois intervalles de longueur
#respective theta_0,theta_1 et theta_2, et de dire si la réalisation de la loi
#uniforme sur [0,1] tombe dans l'intervalle de longuur theta_0, on associe la
#valeur 0, si la réalisation de la loi uniforme sur [0,1] tombe dans 
#l'intervalle de longuur theta_1, on associe la valeur 1 et sinon la valeur 2.

simu012<-function(n,theta0,theta1,theta2)
{
u=runif(n)		#n observations de la loi uniforme sur [0,1]
y=0*(u<=theta0)+1*(u>theta0)*(u<=theta0+theta1)+2*(u>theta0+theta1)
simu012=y
}

#4)
C=simu012(1000,0.2,0.3,0.5)
counts=table(C)
barplot(counts/1000)

#la hauteur des bâtons est proche des probabilités entrées.


#Illustration des grands théorèmes

LGNber<-function(n,k,p)
{
M=matrix(data=rbinom(n*k,1,p),ncol=n)		#simualtion de k jeux de données
								#de longueur n

N=t(apply(M,1,cumsum))		#N est la matrice avec sur chaque ligne les sommes partielles
				#autrement dit si on note x1,x2,..,xn les éléments de la
				#première ligne de M, la première ligne de N est
				#x1,x1+x2,x1+x2+x3,...,x1+x2+x3+..+xn
R=N%*%diag(1/(1:n))		#cela divise tous les éléments de la colonne i par i
				#ainsi on obitient les moyennes empriques à savoir en colonne j
				#on a la moyenne empirique des j premières colonnes
ymin=min(R)
ymax=max(R)
for (i in 1:k)		#tracé des k trajectoires de la moyenne empirique
{
plot(1:n,R[i,],type='l',xlim=c(1,n),ylim=c(ymin,ymax),col=palette()[i])
par(new=TRUE)
}
}


LGNber(700,4,0.2)

#3)
#d'après le théorème central limit la variable considérée converge vers une 
#loi normale d'espérance 0 et de variance p.(1-p)

#la loi de n.Xbar est une binomiale de paramètre n et p car n.Xbar=X1+X2+...+Xn
#et par ailleurs P(sqrt{n}(Xbar-p)<=t)=P(n.Xbar<=sqrt{n}.t+np)

TCLber<-function(n,p)
{
t=seq(-10,10,0.05)
yt=pnorm(t,0,sqrt(p*(1-p)))		#évaluation fonction de répartition limite
yb=pbinom(sqrt(n)*t+n*p,n,p)		#évaluation fonction de répartition de sqrt{n}(Xbar-p)
plot(t,yt,type='l',col='red',xlim=c(-10,10),ylim=c(0,1))
par(new=TRUE)
plot(t,yb,type='s',xlim=c(-10,10),ylim=c(0,1))
}

TCLber(5,0.2)
TCLber(20,0.2)
TCLber(50,0.2)
TCLber(3,0.7)
TCLber(10,0.7)
TCLber(50,0.7)

#4)
#il faut simuler un grand nombre de réalisation de n.Xbar
limitpos<-function(n,k,lambda)
{
M=matrix(data=rbinom(n*k,1,lambda/n),ncol=n)	#simulation des bernoulli
N=apply(M,1,mean)						#k simulations de Xbar
R=n*N								#calcul de n.Xbar
a=max(R)							#calcul de la plus grande valeur observee
Rb=factor(R,levels=0:a)					#transformation en facteur en disant que
								#l'on aurait pu observer des valeurs
								#entre 0 et a
counts=table(Rb)						#tableau de contingence
ct=dpois(0:a,lambda)					#probabilité théorique de la poisson
barplot(height=rbind(x=counts/k,y=ct),beside=TRUE)	#bâtons avec juxtaposé les fréquences
									#observées et les proba théoriques
}

limitpos(500,10000,2)
limitpos(1000,10000,2)

#5)
#le théorème central limit dit que sqrt{n}(Xbar-1/lambda) converge vers une
#loi normale d'espérance 0 et de variance 1/lambda^2

#il faut faire plein de réalisation de sqrt{n}(Xbar-1/lambda)

tclexpo<-function(n,k,lambda)
{
M=matrix(data=rexp(n*k,lambda),ncol=n)	#simulations des exponentielles
N=apply(M,1,mean)					#on obtient k réalisations de Xbar
T=sqrt(n)*(N-1/lambda)				#k réalisations de sqrt{n}(Xbar-1/lambda)
Ts=sort(T)						#on range ces réalisations par ordre croissant
xmin=min(T)
xmax=max(T)
x=seq(xmin,xmax,0.01)				#grille pour évaluer la fucntion de
							#répartition limite
y=pnorm(x,0,1/lambda)				#évaluation de la fonction
							#de répartition limite
yemp=1/k*(1:k)					#valeurs de la fonction de répartition
							#empirique en les points de T
plot(x,y,xlim=c(xmin,xmax),ylim=c(0,1),type='l',col='red')
par(new=TRUE)
plot(Ts,yemp,type='s',xlim=c(xmin,xmax),ylim=c(0,1))
}

tclexpo(2,1000,2)
tclexpo(10,1000,2)
tclexpo(50,1000,2)

#méthode du rejet
#1)
rejet<-function(n)
{
x=c()
l=0
while(l<n)				#boucle jusqu'à obtenir un vecteur de taille n
{
Z=rexp(1,1)
epsilon=2*rbinom(1,1,0.5)-1	#simule une variable qui prend les valeurs -1 et 1
					#avec proba 1/2
Y=epsilon*Z				#1 observation de la variable Y
u=runif(1)				#1 observation de la variable U
gY=1/sqrt(2*pi)*exp(-Y^2/2)/(sqrt(2*exp(1)/pi)*exp(-abs(Y))/2	#calcul de g(Y)
if(gY>=u)				#incrémentation de x si la condition g(Y)>=U vérifiée
{x=c(x,Y)
l=length(x)
}
}
rejet=x
}

A=rejet(1000)
H=hist(A,plot=FALSE,breaks=15) 		#récupérer les éléments de la construction de
					#l'histogramme sans le tracer
limits=H$breaks
xmin=min(limits)			#plus petite limite de classe
xmax=max(limits)			#plus grande limite de classe
x=seq(xmin,xmax,0.01)		#vecteur sur lequel on va évaluer la densité théorique
y=dnorm(x)			#évaluation de la densité
ymax=max(H$density,y)

hist(A,freq=FALSE,breaks=15,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x,y,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique



#Box Muller
#1)
boxmuller<-function(n)
{
u1=runif(n)
u2=runif(n)
X1=sqrt(-2*log(u1))*cos(2*pi*u2)
X2=sqrt(-2*log(u1))*sin(2*pi*u2)
M=matrix(data=0,ncol=2,nrow=n)
M[,1]=X1
M[,2]=X2
boxmuller=M
}

#2)
M=boxmuller(100000)
X1=M[,1]
X2=M[,2]

H1=hist(X1,plot=FALSE) 		#récupérer les éléments de la construction de
					#l'histogramme sans le tracer
limits=H1$breaks
xmin=min(limits)			#plus petite limite de classe
xmax=max(limits)			#plus grande limite de classe
x1=seq(xmin,xmax,0.01)		#vecteur sur lequel on va évaluer la densité théorique
y1=dnorm(x1)				#évaluation de la densité
ymax=max(H1$density,y1)

hist(X1,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x1,y1,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique

H2=hist(X2,plot=FALSE) 		#récupérer les éléments de la construction de
					#l'histogramme sans le tracer
limits=H2$breaks
xmin=min(limits)			#plus petite limite de classe
xmax=max(limits)			#plus grande limite de classe
x2=seq(xmin,xmax,0.01)		#vecteur sur lequel on va évaluer la densité théorique
y2=dnorm(x2)				#évaluation de la densité
ymax=max(H1$density,y2)

hist(X2,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x2,y2,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique


#3)
X3=X1+X2
H3=hist(X3,plot=FALSE) 		#récupérer les éléments de la construction de
					#l'histogramme sans le tracer
limits=H3$breaks
xmin=min(limits)			#plus petite limite de classe
xmax=max(limits)			#plus grande limite de classe
x3=seq(xmin,xmax,0.01)		#vecteur sur lequel on va évaluer la densité théorique
y3=dnorm(x3,0,sqrt(2))				#évaluation de la densité
ymax=max(H1$density,y3)

hist(X3,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x3,y3,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique

par(mfrow=c(1,3))
hist(X1,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x1,y1,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique


hist(X2,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x2,y2,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique

hist(X3,freq=FALSE,xlim=c(xmin,xmax),ylim=c(0,ymax))	#tracé de l'histogramme
par(new=TRUE)		#permet la superposition mais il faut bien définir le jeu d'axe
plot(x3,y3,type='l',col='red',xlim=c(xmin,xmax),ylim=c(0,ymax))
				#tracé de la fonction de densité théorique