########################################## # Trajectoire dans le cas du jeu Pile ou Face # probabilité de succès p=0.5 # nombre de lancers n=10 # réalisation de n lancers Succes_Echec=rbinom(n,1,p) Succes_Echec # nombre de succès succes=cumsum(Succes_Echec) # la fonction suivante permet de mettre deux graphiques dans la même fenêtre par(mfrow=c(2,1)) # on trace la trajectoire du nombre de succès en fonction du nombre de lancers plot((0:(n+1)),c(0,succes,succes[n]),col="red",type="s",xlab="Nombre de lancers",ylab="Nombre de succès obtenus",main="Nombre de succès en fonction du nombre de lancers",xlim=c(0,n+1),ylim=c(0,n)) # fréquence de succès moy_emp=succes/(1:n) # on trace la trajectoire de la fréquence de succès en fonction du nombre de lancers plot((0:(n+1)),c(0,moy_emp,moy_emp[n]),col="blue",type="b",,xlab="Nombre de lancers",ylab="Moyenne empirique",main="Fréquence de succès en fonction du nombre de lancers") ############################# # loi des grands nombres LGN<-function(n,p) { #un seul graphe dans la fenêtre graphique par(mfrow=c(1,1)) # fréquence de succès sur n lancers où p est la probabilité de succès moy_emp=cumsum(rbinom(n,1,p))/(1:n) # on trace la trajectoire de la fréquence de succès en fonction du nombre de lancers plot(1:n,moy_emp,col="blue",type="l",xlab="Nombre de lancers",ylab=" ", main="Fréquence de succès en fonction du nombre de lancers",xlim=c(1,n),ylim=c(0,1)) # on trace sur le même graphique la droite représentant la probabilité de succès lines(c(1,n),c(p,p)) # légende du graphique legend("topright",c(paste("p=",p),"Moyenne empirique"),col=c("black","blue"),lty=c(1),bty="n") } LGN(50,0.5) LGN(500,0.5) LGN(10000,0.5) # Si onveut enregister les graphes dans des fichiers # la commande pdf('nomfichier.pdf')... dev.off() permet d'enregistrer le graphe dans un fichier pdf. #Utiliser la fonction rechercher de votre ordinateur pour voir où a été sauvegardé le fichier pdf("LGN50.pdf") LGN(50,0.5) dev.off() pdf("LGN100.pdf") LGN(500,0.5) dev.off() pdf("LGN1000.pdf") LGN(10000,0.5) ############################################################ # Vitesse de convergence dans la loi des grands nombres #probabilité de succès p=0.5 # nombre de lancers n=500 # nombre de courbes réalisées Nexp=15 # couleur des différentes courbes c=rainbow(Nexp) # droite représentant la probabilité de succès plot(c(1,n),c(p,p),type="l",xlab="Nombre de lancers",ylab="Moyenne empirique",main="Fréquence de succès en fonction du nombre de lancers",xlim=c(1,n),ylim=c(0,1)) # on réalise les différentes courbes for (i in 1:Nexp){ # on calcule la ième fréquence de succès en fonction de n moy_emp=cumsum(rbinom(n,1,p))/(1:n) # on trace la ième courbe lines((1:n),moy_emp,col=c[i]) } # on trace sur le même graphique les courbes d'équation p+/- 1/sqrt(n) IF1=p-1/sqrt(1:n) IF2=p+1/sqrt(1:n) lines(1:n,IF1,lwd=2) lines(1:n,IF2,lwd=2) # légende du graphique legend(300,0.8,c(expression(p +1/sqrt(n)),expression(p -1/sqrt(n))), col=c("black","white"),lty=c(1),bty="n",lwd=2) ########### # Si on veut tester d'autres vitesses de convergence #pdf("vit-LGN1.pdf") c=rainbow(Nexp) plot(c(1,n),c(p,p),type="l",xlab="Nombre de lancers",ylab="Moyenne empirique",main=" ",xlim=c(1,n),ylim=c(0,1)) for (i in 1:Nexp){ moy_emp=cumsum(rbinom(n,1,p))/(1:n) lines((1:n),moy_emp,col=c[i]) } IF1=p-1/(1:n) IF2=p+1/(1:n) lines(1:n,IF1,lwd=2) lines(1:n,IF2,lwd=2) legend(300,0.8,c(expression(p +1/n),expression(p -1/n)),col=c("black","white"),lty=c(1),bty="n",lwd=2) #dev.off() #pdf("vit-LGN2.pdf") c=rainbow(Nexp) plot(c(1,n),c(p,p),type="l",xlab="Nombre de lancers",ylab="Moyenne empirique",main=" ",xlim=c(1,n),ylim=c(0,1)) for (i in 1:Nexp){ moy_emp=cumsum(rbinom(n,1,p))/(1:n) lines((1:n),moy_emp,col=c[i]) } IF1=p-1/log((2:n)) IF2=p+1/log(2:n) lines(2:n,IF1,lwd=2) lines(2:n,IF2,lwd=2) legend(300,0.8,c(expression(p +1/log(n)),expression(p -1/log(n))),col=c("black","white"),lty=c(1),bty="n",lwd=2) #dev.off() #pdf("vit-LGN3.pdf") c=rainbow(Nexp) plot(c(1,n),c(p,p),type="l",xlab="Nombre de lancers",ylab="Moyenne empirique",main=" ",xlim=c(1,n),ylim=c(0,1)) for (i in 1:Nexp){ moy_emp=cumsum(rbinom(n,1,p))/(1:n) lines((1:n),moy_emp,col=c[i]) } IF1=p-log((2:n))/(2:n) IF2=p+log(2:n)/(2:n) lines(2:n,IF1,lwd=2) lines(2:n,IF2,lwd=2) legend(300,0.8,c(expression(p +log(n)/n),expression(p -log(n)/n)),col=c("black","white"),lty=c(1),bty="n",lwd=2) #dev.off() ###################################################################### # Histogramme de la loi Binomiale B(n,p) pour différentes valeurs de n histo<-function(p) { # on sépare la fenêtre graphique pour y placer 2x3 graphes par(mfrow=c(2,3)) # histogramme de la loi B(1,p) n=1 x=0:n plot(x,dbinom(x,n,p),type="h",col="red",lwd=5,xlab="",ylab="",main="Binomiale B(1,p)") # histogramme de la loi B(5,p) n=5 x=0:n plot(x,dbinom(x,n,p),type="h",col="red",lwd=5,xlab=paste("p=",p),ylab="",main="Binomiale B(5,p)") # histogramme de la loi B(10,p) n=10 x=0:n plot(x,dbinom(x,n,p),type="h",col="red",lwd=5,xlab="",ylab="",main="Binomiale B(10,p)") # histogramme de la loi B(20,p) n=20 x=0:n plot(x,dbinom(x,n,p),type="h",col="red",lwd=5,xlab="",ylab="",main="Binomiale B(20,p)") # histogramme de la loi B(50,p) n=50 x=0:n plot(x,dbinom(x,n,p),type="h",col="red",lwd=5,xlab="",ylab="",main="Binomiale B(50,p)") # histogramme de la loi B(100,p) n=100 x=0:n plot(x,dbinom(x,n,p),type="h",col="red",lwd=5,xlab="",ylab="",main="Binomiale B(100,p)") } histo(p=0.5) histo(p=0.3) histo(p=0.9) ########################################################################## # Histogramme de la loi sqrt(n/p(1-p))(Sn/n-p) pour différentes valeurs de n TCL<-function(p) { # nombre de variables B(n,p) simulées pour tracer l'histogramme Nexp=10000 # séparation de la fenêtre graphique par(mfrow=c(2,3)) n=1 echantillon=sqrt(n/(p*(1-p)))*(rbinom(Nexp,n,p)/n-p) hist(echantillon,border="red",freq=FALSE,main=" ",xlab=paste("n=",n," , p=",p),ylab="") n=5 echantillon=sqrt(n/(p*(1-p)))*(rbinom(Nexp,n,p)/n-p) hist(echantillon,border="red",freq=FALSE,main=expression(over(sqrt(n),sqrt(p(1-p)))(sum(over(x[i], n), i == 1, n)-p)),xlab=paste("n=",n," , p=",p),ylab="") # on trace la densité normale N(0,1) par dessus l'histogramme aux=seq(min(echantillon),max(echantillon),by=0.1) lines(aux,dnorm(aux,0,1)) n=10 echantillon=sqrt(n/(p*(1-p)))*(rbinom(Nexp,n,p)/n-p) hist(echantillon,border="red",freq=FALSE,xlab=paste("n=",n," , p=",p),ylab="",main="") aux=seq(min(echantillon),max(echantillon),by=0.1) lines(aux,dnorm(aux,0,1)) n=50 echantillon=sqrt(n/(p*(1-p)))*(rbinom(Nexp,n,p)/n-p) hist(echantillon,border="red",freq=FALSE,xlab=paste("n=",n," , p=",p),ylab="",main="") aux=seq(min(echantillon),max(echantillon),by=0.1) lines(aux,dnorm(aux,0,1)) n=100 echantillon=sqrt(n/(p*(1-p)))*(rbinom(Nexp,n,p)/n-p) hist(echantillon,border="red",freq=FALSE,xlab=paste("n=",n," , p=",p),ylab="",main="") aux=seq(min(echantillon),max(echantillon),by=0.1) lines(aux,dnorm(aux,0,1)) n=1000 echantillon=sqrt(n/(p*(1-p)))*(rbinom(Nexp,n,p)/n-p) hist(echantillon,border="red",freq=FALSE,xlab=paste("n=",n," , p=",p),ylab="",main="") aux=seq(min(echantillon),max(echantillon),by=0.1) lines(aux,dnorm(aux,0,1)) } TCL(0.5) TCL(0.3) TCL(0.9) ################################################## #Comparaison des intervalles de fluctuations comparaisonIF<-function(p,n,alpha) { # un seul graphe dans la fenêtre graphique par(mfrow=c(1,1)) # Calcul du t_alpha t=qnorm(1-alpha/2) # moyenne empirique pour les paramètres n, p moy_emp=cumsum(rbinom(n,1,p))/(1:n) # On trace l'évolution de la moyenne empirique plot((1:n),moy_emp,type="l",xlab="Nombre de lancers",ylab="Moyenne empirique", main=paste("Comparaison des intervalles de fluctuations pour","alpha=",alpha),xlim=c(1,n),ylim=c(0,1)) # Intervalle de fluctuations niveau seconde IF1=p-1/sqrt((1:n)) IF2=p+1/sqrt((1:n)) lines(1:n,IF1,col="blue") lines(1:n,IF2,col="blue") # Intervalle de fluctuations niveau terminale IF1TCL=p-t*sqrt(p*(1-p))/sqrt((1:n)) IF2TCL=p+t*sqrt(p*(1-p))/sqrt((1:n)) lines(1:n,IF1TCL,col="red") lines(1:n,IF2TCL,col="red") # on trace la droite égale à p lines(c(1,n),c(p,p),col="green") # légende du graphique legend("topright",c(paste("p=",p),expression(p+-1/sqrt(n)), expression(p+-t*sqrt(over(p*(1-p),n)))),col=c("green","blue","red"),lty=c(1),bty="n") } comparaisonIF(0.2,500,5/100) comparaisonIF(0.5,500,5/100) comparaisonIF(0.9,500,5/100) #################################### #Intervalle de confiance comparaisonIC<-function(p,n,alpha) { # Calcul du t_alpha t=qnorm(1-alpha/2) # On trace l'évolution de la moyenne empirique moy_emp=cumsum(rbinom(n,1,p))/(1:n) plot((1:n),moy_emp,type="l",xlab="Nombre de lancers",ylab="Moyenne empirique",main="Comparaison des intervalles de confiance",xlim=c(1,n),ylim=c(0,1)) # Intervalle de confiance simple IC1=moy_emp-1/sqrt((1:n)) IC2=moy_emp+1/sqrt((1:n)) lines(1:n,IC1,col="blue") lines(1:n,IC2,col="blue") # Intervalle de confiance avec le TCL IC1TCL=moy_emp-t*sqrt(moy_emp*(1-moy_emp))/sqrt((1:n)) IC2TCL=moy_emp+t*sqrt(moy_emp*(1-moy_emp))/sqrt((1:n)) lines(1:n,IC1TCL,col="red") lines(1:n,IC2TCL,col="red") # vraie valeur du paramètre p lines(c(1,n),c(p,p),col="green") # légende du graphique legend("topright",c(paste("p=",p),expression(bar(X) +-1/sqrt(n)), expression(bar(X) +-1.96*sqrt(over(bar(X)(1-bar(X)),n)))), col=c("green","blue","red"),lty=c(1),bty="n") #legend("topright",c(paste("p=",p),expression(bar(X) +-1/sqrt(n)),expression(bar(X) +-1.96*sqrt(over(bar(X)(1-bar(X)),n))), expression(bar(X)==sum(over(x[i], n), i == 1, n))),col=c("green","blue","red","white"),lty=c(1),bty="n") } comparaisonIC(p=0.2,n=500,alpha=5/100) comparaisonIC(p=0.5,n=500,alpha=5/100) comparaisonIC(p=0.9,n=500,alpha=5/100) ######################################################################## # Probabilité de sortir de l'intervalle IF pour différentes valeurs de p n=50 Nexp=10000 moy_emp=rep(0,Nexp) proba=rep(0,Nexp) par(mfrow=c(3,1)) p=0.1 IF1=p-1/sqrt(n) IF2=p+1/sqrt(n) moy_emp=rbinom(Nexp,n,p)/n proba=cumsum(moy_empIF2)/(1:Nexp) plot((1:Nexp),proba,type="l",xlab="p=0.1 ",ylab=" ",main="Probabilité que la moyenne empirique soit hors de l'I.F.",col="blue") lines(c(-1,Nexp),c(proba[Nexp],proba[Nexp]),col="red") p=0.5 IF1=p-1/sqrt(n) IF2=p+1/sqrt(n) moy_emp=rbinom(Nexp,n,p)/n proba=cumsum(moy_empIF2)/(1:Nexp) plot((1:Nexp),proba,type="l",xlab="p=0.5",ylab=" ",main=" ",col="blue") lines(c(-1,Nexp),c(proba[Nexp],proba[Nexp]),col="red") p=0.9 IF1=p-1/sqrt(n) IF2=p+1/sqrt(n) moy_emp=rbinom(Nexp,n,p)/n proba=cumsum(moy_empIF2)/(1:Nexp) plot((1:Nexp),proba,type="l",xlab="p=0.9",ylab=" ",main=" ",col="blue") lines(c(-1,Nexp),c(proba[Nexp],proba[Nexp]),col="red") ########################### n=100 p=(1:10)/10 Nexp=10000 proba=rep(0,length(p)) moy_emp=0 for (i in 1:length(p)){ IF1=p[i]-1/sqrt(n) IF2=p[i]+1/sqrt(n) moy_emp=rbinom(Nexp,n,p[i])/n proba[i]=sum(moy_empIF2)/Nexp } par(mfrow=c(1,1)) plot(p,proba,type="l",col="red",xlab="probabilité de succès",ylab="",main=" Probabilité que la moyenne empirique soit hors de l'I.F. en fonction de p")