################################################################################# ##### ##### Correction TD4 - Simulation de variables aleatoires ##### ################################################################################## rm(list=ls()) ############################################ # Avec les fonctions usuelles ############################################ # 1. help(rnorm) # rnorm simule n variables gaussiennes de moyenne et d'ecart type donne, par defaut, respectivement 0 et 1. # dnorm donne la densite gaussienne en un vecteur x de moyenne et d'ecart type donne, par defaut, respectivement 0 et 1. # pnorm donne de la meme maniere la fonction de repartition (on peut changer une option pour que ce soit la fonction de queue). # qnorm donne les quantiles. Une option permet aussi de transformer p en log(p). # De maniere generale, toutes les fonctions qui commencent par r donneront des variables aleatoires, une fonction qui commence par d donnera la densite pour les variables continues ou la probabilite d'avoir telle valeur pour les lois discretes, une fonction qui commence par p la fonction de repartition et par q les quantiles. # Ex : rpois, dpois, ppois, qpois pour la loi de Poisson. # 2. rbinom(1,4,0.2) # Donne une variable binomiale avec n=4 et p=0.2. rbinom(1,1,0.2) # On voit la loi de Bernoulli comme un cas particulier de la binomiale de parametre n=1 et p=0.2. runif(1) # Donne une varaible uniforme sur [0,1], on peut aussi specifier les bornes dans les options, la c'est ce qui se passe par defaut. rexp(1,2) # Donne une variable exponentielle de parametre 2. rpois(1,3) # Donne une variable de Poisson de parametre 3. # NB : le premier 1 dans les commandes est la pour dire que l'on veut une seule variable aleatoire. # On pourrait en demander plus. Dans ce cas R renverrait un vecteur de variables aleatoires. # 3. X= rnorm(1000, mean=2, sd=2) # NB : la fonction rnorm attend l'ecart-type et non la variance (donc ici 2 et non 4). # 4. # Variable quantitative continue, on trace donc l'histogramme des frequences. hist(X, freq=FALSE, main='Simulation N(2,4)', xlab='variables', ylab='Densite') x= seq(-5, 10, 0.1) lines(x, dnorm(x,mean=2,sd=2), col='red') legend('left', c('simulation','theorie'), lty=c(1,1), col=c('black','red')) # Attention a bien renormaliser l'histogramme sinon la comparaison histogramme/densite n'a aucun sens. ############################################ # A partir d'une variable uniforme ############################################ # 1. # La fonction de repartition d'une variable exponentielle de parametre a est t-> F(t)=1-exp(-a*t) pour t>0, et 0 sinon. # Son inverse est u-> -log(1-u)/a pour 00) # Si on a trouve au moins un indice negatif, alors... { print('Theta doit etre positif ou nul.') vartest= 1 } if(abs(sum(theta)-1)>10^(-15)) # On regarde ici si la somme des theta[i] ne fait pas 1 (a la precision machine pres). # De maniere generale faire attention aux tests d'egalite sur des reels et les remplacer par des tests de difference plus petite que la precision machine. { print('Theta doit etre un vecteur de somme 1.') vartest= 1 } if(vartest==0) # Si on a bien passe tous les tests. { U= runif(n) # Simulation des n variables uniformes sur [0,1]. X= 0*(U=theta[1])&(U<(theta[1]+theta[2])))+2*(U>=(theta[1]+theta[2])) # (U 1/n (X_1+....+X_n) que l'on souhaite dessiner. On en prend plusieurs pour observer les fluctuations. # On calcule ce qu'il faut. Attention de n'utiliser que p et nmax, et non directement les valeurs que l'on a mise dedans. nbp= length(p) # Evite de tout changer si on change le vecteur des p. par(mfrow=c(1,nbp)) # On va faire automatiquement autant de dessins que le nombre de p que l'on veut tester. for(i in 1:nbp) # Boucle sur la probabilite p. { # Trace la ligne vers laquelle les trajectoires doivent converger. plot(c(1,nmax), c(p[i],p[i]), type='l', col='red', ylim=c(0,1), xlim=c(1,nmax), main=paste('p=',p[i]), xlab='n', ylab='1/n(X1+...+X_n)') for(j in 1:Ntraj) # Boucle sur les trajectoires. { X= rbinom(nmax, 1, p[i]) # Simulation d'un grand vecteur contenant des variables de Bernoulli de parametre p[i]. traj= cumsum(X)/(1:nmax) # Calcul de la trajectoire complete. lines(1:nmax, traj) # Trace de cette trajectoire sur le graphe. } # Pour la legende c'est un peu complique car ca depend de p. Ce n'est pas grave si elle n'y est pas dans un premier temps. # En voici une qui fonctionne a tous les coups : if(p[i]<0.5) # Si p est petit, on met la legende en haut. Sinon, en bas. { legend(0.1*nmax, 0.9,c('trajectoires','valeur limite'), lty=c(1,1), col=c('black','red')) }else{ legend(0.1*nmax, 0.2,c('trajectoires','valeur limite'), lty=c(1,1), col=c('black','red')) } } # On voit bien que les trajectoires ne sont pas les memes, que ca fluctue mais que pour n grand on est toujours proche de la limite. # Remarquez que les trajectoires sont generalement plus dispersees vers p=0.5, car c'est la ou 1/n (X1+...+Xn) a la plus grande variance. # 2. # Ici la LGN dit que 1/n (X_1+....+X_n) tend vers 1/lambda car X_1 suit une loi exponentielle de parametre lambda. lambda= c(0.5,1,3) # On se donne un ensemble de lambda possibles. nmax= 1000 # Le nombre maximal pour n. Ntraj= 3 # Le nombre de trajectoires n-> 1/n (X_1+....+X_n) que l'on souhaite dessiner. On en prend plusieurs pour observer les fluctuations. nbl= length(lambda) par(mfrow=c(1,nbl)) for(i in 1:nbl) # Boucle sur le parametre lambda. { # Trace la ligne vers laquelle les trajectoires doivent converger. plot(c(1,nmax), c(1/lambda[i],1/lambda[i]), type='l', col='red', ylim=c(0,3), xlim=c(1,nmax), main=paste('lambda=',lambda[i]), xlab='n', ylab='1/n(X1+...+X_n)') for(j in 1:Ntraj) # Boucle sur les trajectoires. { X= rexp(nmax,lambda[i]) # Simulation d'un grand vecteur contenant des variables exponentielles de parametre lambda[i]. traj= cumsum(X)/(1:nmax) # Calcul de la trajectoire complete. lines(1:nmax,traj) # Trace de cette trajectoire sur le graphe. } legend('topright',c('trajectoires','valeur limite'),lty=c(1,1),col=c('black','red')) # je me mets presque au bout (7/10 de nmax) et vu les parametres que j'ai pris je serai toujours sous 2.5 } # Meme commentaire qualitatif, la variance est plus grande pour les lambda petits. # 3. # la loi de n*barX est celle d'une binomiale de parametre n et p. Avec pbinom on a donc acces a # la fonction x-> P(n*barX <= x). # Mais P(n*barX <= x) = P(sqrt(n)*(barX-p) <= x/sqrt(n)-sqrt(n)*p) donc on peut finalement deduire de pbinom la fonction de repartition cherchee. n= c(10,100,1000) # On se donne un vecteur de n. p= c(0.3,0.5,0.8) # On se donne un vecteur de p. nbp= length(p) nbn= length(n) par(mfrow=c(nbp,nbn)) # On va faire un trace n par p. for(i in 1:nbp) # Boucle sur les p { for(j in 1:nbn) # Boucle sur les n { x_binom= 0:n[j] # Les abscisses pour n*barX (on met toutes les valeurs potentielles de n*barX). y= pbinom(x_binom, n[j], p[i]) # Les ordonnees (qui elles ne bougent pas) de la fonction de repartition de n*barX, binomilae de parametre n[j] et p[i]. x_renorm= x_binom/sqrt(n[j])-sqrt(n[j])*p[i] # Les abscisses renormalisees pour sqrt(n)*(barX-p). x= seq(-3,3,0.1) # Les abscisses pour la fonction de repartion de la loi normale. plot(x, pnorm(x,mean=0,sd=sqrt(p[i]*(1-p[i]))), type='l', col='red', main=paste('TCL avec p=', p[i],' et n=', n[j]), xlab='sqrt(n)(bar X-p)', ylab='fonction de repartition') # On trace la fonction de repartition de la loi limite (loi normale centree et de variance p*(1-p)). lines(x_renorm, y, type='s') # On superpose ensuite la fonction de repartition de sqrt(n)(bar X-p). Attention variable qui ne prend qu'un nombre fini de valeurs, comme n barX qui est une binomiale, d'ou le type 's' pour stairs (escalier en anglais). # legend(-3, 0.8, c('fonction de repartition','gaussienne limite'), lty=c(1,1), col=c('black','red')) # On enleve la legende, sinon elle prend toute la place... } } # On observe bien la convergence de la fonction de repartition vers celle de la loi normale lorsque n grandit. # 4. # On va faire la meme chose avec un diagramme en baton. n= c(10,100,1000) # On se donne un vecteur de n. lambda= c(1,2,4) # On se donne un vecteur de lambda. Nsimu= 1000 nbl= length(lambda) # Extraction des longueurs qu'il faut pour les plots. nbn= length(n) par(mfrow=c(nbl,nbn)) # Trace n par lambda. for(i in 1:nbl) # Boucle sur les lambda. { for(j in 1:nbn) # Boucle sur les n. { x_limite= qpois(1-10^(-4), lambda[i]) # On recupere le quantile d'ordre 1-10^(-4) de la variable poisson limite. Si l'approximation est bonne, tout ce qui se passe apres est negligeable. x_poiss= 0:x_limite # Les valeurs possible de la loi de Poisson jusqu'a ce quantile. p= lambda[i]/n[j] # On prend ce p par hypothese. nbarX= rbinom(Nsimu,n[j],p) # On simule Nsimu fois la variable binomiale. tx= table(factor(nbarX,levels=0:x_limite))/Nsimu # On la transforme en facteur avec les bons niveaux qui nous interessent. b= barplot(tx, ylim=c(0,1), main=paste('n=',n[j],' et p=',lambda[i],'/n'), xlab='valeurs', ylab='probabilite') # Affichage du diagramme en baton pour les donnees simulees. points(b, dbinom(x_poiss,n[j],p), pch=24, col='cyan', bg='cyan') # Les vraies valeurs des probabilites pour la loi binomiale. points(b, dpois(x_poiss,lambda[i]), pch=20, col='red', bg='red') # Les valeurs des probas vers lesquelles on doit converger (loi de Poisson). legend('topright', c('simulation','vraie valeur','Poisson limite'), pch=c(22,24,20), col=c('black','cyan','red'), pt.bg=c('grey','cyan','red')) # NB : on a fait en sorte que le symbole du 'cyan' soit plus gros que le symbole du 'red'. Meme si le point rouge tombe dessus, on voit les deux ! } } # On observe bien une convergence de la loi binomiale vers la loi de Poisson lorsque n grandit. # 5. # On va approcher la fonction de repartition theorique a n fixe par la courbe des frequences cumulees sur un echantillon de taille Nsimu. n= c(10,100,1000) # On se donne un vecteur de n. a= c(1,3) # On se donne un vecteur de a. Nsimu= 1000 nba= length(a) # Extraction des longueurs qu'il faut pour les plots. nbn= length(n) par(mfrow=c(nba,nbn)) for(i in 1:nba) # Boucle sur les a. { for(j in 1:nbn) # Boucle sur les n. { simu= matrix(data=rexp((Nsimu*n[j]), rate=a[i]), ncol=Nsimu, nrow=n[j]) # Creation d'une grande matrice remplie de variables exponentielles iid avec les dimensions voulues. Dans la colonne s je vais avoir les 'X_i' pour i allant de 1 a n pour la simu numero s. # La commande colMeans donne les moyennes colonne par colonnes, donc les variables qui nous interessent sont : var= sqrt(n[j])*(colMeans(simu) - 1/a[i]) x= seq(-3,3,0.1) plot(x, pnorm(x,mean=0,sd=1/a[i]), type='l', main=paste('a=', a[i],' et n=',n[j]), xlab='sqrt(n)(barX-1/a)', ylab='fonction de repartition', col='red') # On trace la fonction de repartition de la gaussienne limite. h= hist(var,plot=FALSE) # On recupere les classes et leurs effectifs. lines(h$breaks,c(0,cumsum(h$counts)/Nsimu)) # On trace la courbe des frequences cumulees comme au td 2. # Ne pas oublier le 0 au debut du vecteur des frequences sinon les objets n'ont pas la meme taille et le plot ne fonctionne pas. legend('bottomright', c('fonction de repartition','gaussienne limite'), lty=c(1,1), col=c('black','red')) } } ################################### ### Methode du rejet ################################### # 1. # On commence par coder la fonction g de l'enonce. g_gauss<-function(x) { C= sqrt(pi/(2*exp(1))) # constante C, ici egale a 1/a fx= exp(-x^2/2)/sqrt(2*pi) f0x= exp(-abs(x))/2 gx= C*fx/f0x return(gx) } monrejet_gauss <-function(n) { var= rep(0,n) # On initialise le vecteur des resultats. Il est de taille n et ne contient que des 0 pour l'instant. for(i in 1:n) { Z= rexp(1,1) # Simule la variable Z qui est une exponenitelle de parametre 1. epsilon= 2*rbinom(1,1,0.5)-1 # Cree une variable qui prend deux valeurs -1 et 1 avec proba 0.5, a partir d'une bernoulli qui prend deux valeurs 0 et 1 avec proba 0.5. Y= epsilon*Z # Calcul du Y correspondant. U= runif(1) # Simule la variable U qui est une uniforme sur [0,1]. gY= g_gauss(Y) # Calcul de g(Y). while(U>gY) # Tant que la condtion U<=g(Y) n'est pas verifiee, on recommence jusqu'a ce qu'elle le soit. { Z= rexp(1,1) epsilon= 2*rbinom(1,1,0.5)-1 Y= epsilon*Z U= runif(1) gY= g_gauss(Y) } # Ici on est sur d'avoir U et Y tels que U<=g(Y). On stocke alors ce Y a la position i. var[i]=Y } return(var) } monrejet_gauss(10) # On regarde rapidement si ca fonctionne. # 2. # On procede comme d'habitude : comparaison histogramme/densite. # On simule 1000 variables avec la methode du rejet. X= monrejet_gauss(1000) hist(X, freq=FALSE, main='Methode du rejet pour la loi N(0,1)', xlab='x', ylab='densite', ylim=c(0,dnorm(0)*1.1)) x= seq(-3,3,0.1) lines(x, dnorm(x), col='red') legend('topleft', c('simulation','theorie'), lty=c(1,1), col=c('black','red')) # On retrouve bien la loi normale centree reduite. ######################################################## ## L'algorithme de Box-Muller ######################################################### # 1. mon_BoxMuller<-function(n) { U1= runif(n) # On tire les variables aleatoires uniformes. U2= runif(n) X1= sqrt(-2*log(U1))*cos(2*pi*U2) # Les formules demandees. X2= sqrt(-2*log(U1))*sin(2*pi*U2) return(rbind(X1,X2)) # on renvoie une matrice a deux lignes et n colonnes } mon_BoxMuller(10) # On verifie que ca fonctionne. # 2. # Comme avant, je prepare pour la question #3 les 3 plots X= mon_BoxMuller(1000) par(mfrow=c(1,3)) hist(X[1,], freq=FALSE, main='BoxMuller X1', xlab='x', ylab='densite', ylim=c(0,dnorm(0)*1.1)) x=seq(-3, 3, 0.1) lines(x, dnorm(x), col='red') legend('topleft', c('simulation','theorie'), lty=c(1,1), col=c('black','red')) hist(X[2,], freq=FALSE, main='BoxMuller X2', xlab='x', ylab='densite', ylim=c(0,dnorm(0)*1.1)) x=seq(-3, 3, 0.1) lines(x, dnorm(x), col='red') legend('topleft', c('simulation','theorie'), lty=c(1,1), col=c('black','red')) # 3. hist(X[1,]+X[2,], freq=FALSE, main='BoxMuller X1+X2', xlab='x', ylab='densite', ylim=c(0,dnorm(0)*1.1)) x=seq(-5, 5, 0.1) # NB : on change les x car variance plus grande donc distribution plus evasee lines(x, dnorm(x,mean=0,sd=sqrt(2)), col='red') legend('top', c('simulation','theorie'), lty=c(1,1), col=c('black','red')) # X1+X2 semble bien etre de loi N(0,2).