<<CheminPackLocaux`

Collège du Sud, Bulle

Applications des mathématiques

Simulation

Version pour Mathematica
Edition 2009/2010
Marcel Délèze, Eugène Pasquier
www.collegedusud.ch/applmaths/

§ 1	Nombres pseudo-aléatoires
	à distribution spécifiée

§ 1.1	Introduction

Pour résoudre un problème dans lequel intervient le hasard, il existe trois familles de méthodes:

Méthode théorique

Le calcul des probabilités sera étudié plus tard, en 4-ème année, dans le cadre du cours de mathématiques. Nous nous contentons ici d'en faire une simple mention.
Lorsque le problème est assez simple, il peut être résolu d'une manière exacte et complète par le calcul des probabilités. Mais, pour des problèmes plus complexes, il est généralement difficile, voire même impossible, de donner une solution explicite. Il est utile de connaître d'autres approches.

Méthode expérimentale

Il est parfois possible d'organiser une série d'expériences.  Par exemple, lancer un dé et noter les résultats obtenus est facile. Mais, pour obtenir une certaine précision, il faut faire un grand nombre d'expériences ce qui est coûteux.

Simulation (ou méthode de Monte Carlo)

Il existe des algorithmes capables de produire des "nombres pseudo-aléatoires" c'est-à-dire des nombres qu'on peut utiliser à la place des nombres donnés par le dé. L'ordinateur étant capable d'effectuer un grand nombre d'expériences en peu de temps, on peut obtenir des résultats plus précis et à meilleur compte qu'avec la méthode d'expérimentation.

Alors que les nombres aléatoires sont engendrés par un phénomène aléatoire (par exemple par tirage à la roulette), les nombres pseudo-aléatoires sont calculés au moyen d'un algorithme; ils offrent néammoins les mêmes caractéristiques que les nombres aléatoires et peuvent avantageusement les remplacer dans les calculs.

On appelle "méthode de Monte-Carlo" toute méthode de calcul utilisant des nombres pseudo-aléatoires. Cette dénomination se réfère à la roulette du célèbre casino de Monte-Carlo. Parmi les applications des méthodes de Monte-Carlo, citons
 *	la programmation de jeux de hasard (dés, cartes, lotto, roulette, ...);
 *	la simulation de phénomènes aléatoires (
			cheminements d'un neutron au sein d'une matière fissile,
			fluctuations de charge d'un central téléphonique,
			longueur de la file d'attente à un guichet, ...).

§ 1.2	Nombres pseudo-aléatoires à distribution discrète

Distribution de Bernoulli

n=40;

?BernoulliDistribution

Par exemple, une expérience aléatoire consistant à lancer un dé, le nombre de six est

six=BernoulliDistribution[(1/6)];
RandomInteger[six,n]

Distribution binomiale

?BinomialDistribution

Par exemple, une expérience aléatoire consistant à lancer 10 dés, le nombre de six est

six10=BinomialDistribution[10,(1/6)];
RandomInteger[six10,n]

Distribution discrète uniforme

?DiscreteUniformDistribution

Par exemple, une expérience aléatoire consistant à lancer 1 dé, le nombre de points est

de=DiscreteUniformDistribution[{1,6}];
RandomInteger[de,n]

Distribution géométrique

?GeometricDistribution

Par exemple, une expérience aléatoire consistant à lancer un dé jusqu'à ce que l'on obtienne un six, le nombe de lancers avant d'obtenir le permier six est

jets=GeometricDistribution[(1/6)];
RandomInteger[jets,n]

Si on désire que le dernier jet, qui est un six, soit compté,

1+RandomInteger[jets,n]

§ 1.3	Nombres pseudo-aléatoires à distribution continue

n=40;

Distribution uniforme

?UniformDistribution

Par exemple, pour créer un échantillon de n mesures principales d'angles, sous l'hypothèse que la distribution est uniforme:

angle=UniformDistribution[{-π,π}];
RandomReal[angle,n]

Distribution normale

?NormalDistribution

Par exemple, pour créer un échantillon des tailles de n personnes, de taille moyenne 1.75 m d'écart-type 0.2 m, sous l'hypothèse que la distribution des tailles est normale:

taille=NormalDistribution[1.75,0.2];
RandomReal[taille,n]

§ 1.4	Estimation de l'erreur

Pour de grandes valeurs de n, la moyenne empirique de n tirages indépendants a une distribution approximativement normale. En notant m la moyenne empirique, s l'écart-type empirique et (Overscript[(Subscript[σ, M]), ∧])=(s/(Sqrt[n-1])) l'estimation de l'écart-type de la moyenne M, nous aurons approximativement
	la probabilité que l'intervalle [m-(Overscript[(Subscript[σ, M]), ∧]),m+(Overscript[(Subscript[σ, M]), ∧])]  contienne μ est de 68.3 % environ;
	la probabilité que l'intervalle [m-2(Overscript[(Subscript[σ, M]), ∧]),m+2(Overscript[(Subscript[σ, M]), ∧])] contienne μ est de 95.5 % environ;
	la probabilité que l'intervalle [m-3(Overscript[(Subscript[σ, M]), ∧]),m+3(Overscript[(Subscript[σ, M]), ∧])] contienne μ est de 99.7 % environ.

Justification:
Voir Probabilités II, § 2.6

§ 1.5	Gestion du temps de calcul

Recherche de l'efficacité

Pour générer des tirages aléatoires, utilisez si possible les distributions de Mathematica et la commande RandomArray[distr, n].

Si la variable est à valeurs {0, 1}, calculez l'écart-type au moyen de la formule     s=(Sqrt[m (1-m)]).

Si la variable est discrète et ne comporte qu'un nombre fini de modalités, regroupez les données puis calculez la moyenne et l'écart-type sur les données groupées.

Contrôle de la durée

La précision s'améliore avec la taille de l'échantillon. Mais la durée des calculs ne peut pas dépasser la limite du supportable, définie par l'utilisateur. Nous adopterons la règle suivante: la taille de l'échantillon est choisie de telle manière que la durée de la simulation est de 1 minute environ.

Pour déterminer le nombre d'expériences, nous effectuons un calcul préalable qui consiste à chronométrer la durée d'un "certain nombre" d'expériences. A partir de là, on calcule combien d'expériences on peut réaliser en 1 minute.

Exercice 1 - 1

Dans le but d'effectuer la simulation suivante : " Tirer 20 nombres réels pseudo-aléatoires uniformément répartis entre 0 et 1 et afficher ceux qui sont compris entre 0.4 et 0.6", un étudiant a écrit le programme suivant :

echant=RandomReal[{0,1},20];

critQ[x_]:=0.4`≤RandomReal[{0,1}]<0.6`

Select[echant,critQ]

Pourquoi le programme ne donne-t-il pas le résultat escompté ?
Comment le corriger pour qu'il accomplisse la tâche désirée ?

Exercice 1 - 2                    (facultatif)

Ecrivez un programme qui distribue aléatoirement les 36 cartes du jeu de jass à quatre joueurs.

Indications

a)	Le début du programme qui construit la liste ordonnée des cartes vous est donné:
		accédez au site  http://www.collegedusud.ch/app/applmaths/
		suivez les liens  Documents Mathematica / Annexes: Simulation
		et téléchargez le cahier  1-2_donnees_exercice.nb

b)	Pour battre les cartes, on peut utiliser l'idée suivante:
	Charger le module de combinatoire:

<<Combinatorica`

La fonction RandomPermutation[k] produit une permutation aléatoire de la liste
	{1,2,3,...,k}. Par exemple,

RandomPermutation[10]

Il s'agit d'appliquer la permutation RandomPermutation[36] à la liste des cartes
	pour obtenir un jeu brassé.

c)	Pour le joueur numéro 1, extraire du jeu les cartes numéros 1, 2, 3, 13, 14, 15, 25, 26, 27;
	pour le joueur numéro 2, extraire du jeu les cartes numéros 4, 5, 6, 16, 17, 18, 28, 29, 30;
	etc.

§ 2	Exemples de simulations

§ 2.1	Deux six consécutifs

Combien de fois faut-il lancer le dé, en moyenne, jusqu'à ce que l'on obtienne 2 six consécutifs?

Comme question préalable, lançons d'abord le dé jusqu'à l'apparition du premier six.

?GeometricDistribution

n6=GeometricDistribution[(1/6)];

Le nombre de lancers d'un dé pour obtenir un six est donné par

1+RandomInteger[n6,20]

Si, dans la liste précédente, on veut déterminer le nombre de lancers pour obtenir deux six consécutifs (pour la première fois), on peut faire comme suit:
	*	il faut sommer le nombre de coups : s=11+3+3+... 
	*	il faut interrompre la sommation au premier 1 rencontré, à partir du deuxième rang :
		  -	le premier terme doit toujours être pris, même si c'est 1;
		  -	le dernier terme est 1 et fait partie de la somme;
	*	le résultat est s=11+3+3+6+1.

Le générateur de nombres pseudo-aléatoires qui simule "le nombre de fois qu'il faut lancer le dé pour obtenir deux six consécutifs" est

deuxSix[]:=Module[{s,t},
		                       s=1+Random[n6];t=1+Random[n6];s+=t;
		                        While[t>1,t=1+Random[n6];s+=t];
		                         Return[s]]

Mesurons d'abord la durée d'une "petite" simulation:

duree=Timing[
		                           n=5000;
		                           x=Table[deuxSix[],{n}]
	              ][[1]]

Accordons maintenant une minute à l'exécution de la simulation principale:

n=n*Round[((60 )/duree)]

x=Table[deuxSix[],{n}];

m=N[Mean[x]]

StandardDeviationMLE[data_]:=StandardDeviation[data]Sqrt[(Length[data]-1)/Length[data]]

s=StandardDeviationMLE[x];
sm=N[(s/(Sqrt[n-1]))]

Interprétation : avec une probabilité de 95 %, le nombre cherché μ  est situé dans l'intervalle

{m-2 sm, m+2 sm}

Ce problème peut être résolu par le calcul des probabilités. Le résultat exact est μ=42.

§ 2.2	Aiguille de Buffon

On lance une aiguille sur un plancher. La longueur de l'aiguille est 1 et la largeur d'une planche est 2.
Quelle est la probabilité que l'aiguille touche 2 planches ?

aiguille[]:=Module[{x,y,z,a},x=20 RandomReal[];y=12 RandomReal[];a=RandomReal[{-(π/2),(π/2)}];Line[{{x-0.5` Cos[a],y-0.5` Sin[a]},{x+0.5` Cos[a],y+0.5` Sin[a]}}]]

Show[Graphics[{Table[Line[{{0,i},{20,i}}],{i,1,12,2}],Table[aiguille[],{i,1,16}]},AspectRatio->Automatic]]

Critère de réussite

Dans un premier temps, organisons une seule expérience.
Notons
	d = distance entre le centre de gravité de l'aiguille et la droite la plus proche;
	θ = angle non orienté entre l'aiguille et la perpendiculaire aux droites.

Show[Graphics[{Circle[{0,0},0.14`,{1-π,-(π/2)}],Circle[{0,0},0.16`,{1-π,-(π/2)}],Text["θ",0.19` {Cos[0.5`-((3 π)/4)],Sin[0.5`-((3 π)/4)]}],Line[{{-0.5` Cos[1],-0.5` Sin[1]},{0.5` Cos[1],0.5` Sin[1]}}],Line[{{-0.6`,-0.2`},{1.3`,-0.2`}}],Arrowheads[{-.02,.02}],Arrow[{{-0.3`,0},{-0.3`,-0.2`}}],Text["d",{-0.35`,-0.1`}],Arrow[{{0.3`,0},{0.3`,-0.5` Sin[1]}}],Text["\!\(\*FractionBox[\(1\), \(2\)]\)cos(θ)",{0.5`,-0.35` Sin[1]}],{Dashing[{0.01`,0.01`}],Line[{{0,-0.6`},{0,0}}],Line[{{-0.6`,0},{1,0}}],Line[{{-0.5` Cos[1],-0.5` Sin[1]},{1,-0.5` Sin[1]}}]}},AspectRatio->Automatic,ImageSize->{500,250}]]

Show[Graphics[{Circle[{0,0},0.1`,{0.3`-π,-(π/2)}],Circle[{0,0},0.11`,{0.3`-π,-(π/2)}],Text["θ",0.15` {Cos[0.15`-((3 π)/4)],Sin[0.15`-((3 π)/4)]}],Line[{{-0.5` Cos[0.3`],-0.5` Sin[0.3`]},{0.5` Cos[0.3`],0.5` Sin[0.3`]}}],Line[{{-0.8`,-0.2`},{1.3`,-0.2`}}],Arrowheads[{-.02,.02}],Arrow[{{-0.7`,0},{-0.7`,-0.2`}}],Text["d",{-0.75`,-0.1`}],Arrow[{{0.3`,0},{0.3`,-0.5` Sin[0.3`]}}],Text["\!\(\*FractionBox[\(1\), \(2\)]\)cos(θ)",{0.5`,-0.3` Sin[0.3`]}],{Dashing[{0.01`,0.01`}],Line[{{0,-0.3`},{0,0}}],Line[{{-0.8`,0},{1,0}}],Line[{{-0.5` Cos[0.3`],-0.5` Sin[0.3`]},{1,-0.5` Sin[0.3`]}}]}},AspectRatio->Automatic,ImageSize->{500,150}]]

On voit sur les deux graphiques précédents que
	l'aiguille coupe un interstice entre deux planches si et seulement si   (1/2) cos(θ)>d

dist=UniformDistribution[{0,1}];

angle=UniformDistribution[{0,(π/2)}];

Mesurons d'abord la durée d'une "petite" simulation:

duree=Timing[
		                n=10000;
		                d=RandomReal[dist,n];
		                θ=RandomReal[angle,n];
		                succes=Count[Table[0.5 Cos[θ[[i]]]>d[[i]],{i,1,n}],True]
	              ][[1]]

Accordons maintenant une minute à l'exécution de la simulation principale:

n=n*Round[(60/duree)]

d=RandomReal[dist,n];
θ=RandomReal[angle,n];

succes=Count[Table[0.5 Cos[θ[[i]]]>d[[i]],{i,1,n}],True]

m=N[(succes/n)]

s=(Sqrt[m(1-m)]);
sm=N[(s/(Sqrt[n-1]))]

Interprétation : avec une probabilité de 95 %, le nombre cherché est situé dans l'intervalle

{m-2 sm, m+2 sm}

Ce problème peut être résolu par le calcul des probabilités. Le résultat exact est μ=(1/π)≃0.31831.

§ 2.3	Calcul d'un volume par la méthode de Monte-Carlo

Considérons le corps   C = { (x, y, z) | 0 < x < (3/2),  0 < y < 1 et 0 < z < 4 - (x^2) - (y^2) }. 
Quel est son volume ?

Plot3D[4-(x^2)-(y^2),{x,0,1.5`},{y,0,1},PlotRange->{0,4},BoxRatios->{1,1,1}]

Une expérience élémentaire consiste à placer un point (x, y, z)  au hasard dans le parallélépipède avec

0≤x<1
0≤y<1.5
0≤z<4

L'expérience est réussie si le point est situé au-dessous de la surface z = f(x,y) c'est-à-dire si

z<4-(x^2)-(y^2)

La fréquence de réussite m est égale au rapport entre le volume "V" cherché et le volume "boite" du parallélépipède:

m=(V/boite)       ⟺        V=m*boite

absc=UniformDistribution[{0,1.5}];

ordon=UniformDistribution[{0,1}];

cote=UniformDistribution[{0,4}];

Mesurons d'abord la durée d'une "petite" simulation:

duree=Timing[
		                n=10000;
		                x=RandomReal[absc,n];
		                y=RandomReal[ordon,n];
		                z=RandomReal[cote,n];
		                succes=Count[Table[z[[i]]<4-((x[[i]])^2)-((y[[i]])^2),{i,1,n}],True]
	              ][[1]]

Accordons maintenant une minute à l'exécution de la simulation principale:

n=n*Round[(60/duree)]

x=RandomReal[absc,n];
y=RandomReal[ordon,n];
z=RandomReal[cote,n];

succes=Count[Table[z[[i]]<4-((x[[i]])^2)-((y[[i]])^2),{i,1,n}],True]

m=N[(succes/n)]

s=(Sqrt[m(1-m)]);
sm=N[(s/(Sqrt[n-1]))]

boite=1*1.5*4;
vol=m*boite

incert=sm*boite

Interprétation : avec une probabilité de 95 %, le nombre cherché est situé dans l'intervalle

{vol-2 incert, vol+2 incert}

Ce problème peut être résolu par calcul intégral. Le résultat exact est 4.375.

Calcul des probabilités exactes

Le calcul des probabilités exactes des exemples des paragraphes 2.1, 2.2 et 2.3 sort des objectifs de ce cours. Les intéressés (en particulier le professeur) peuvent trouver ces calculs dans les annexes:
sur le site  http://www.collegedusud.ch/app/applmaths/
	suivez les liens  Documents Mathematica / Annexes: Simulation
	et téléchargez le cahier  2S-Probabilites_exactes.nb

Exercices du § 2

Exercice 2 - 1

Perpendiculairement à un vaste grillage, on lance une balle au hasard.
La balle a un rayon de 0.6 et les mailles du filet sont rectangulaires de 4 sur 3.

Show[Graphics[{Table[Line[{{i,0},{i,16}}],{i,0,25,4}],Table[Line[{{0,i},{25,i}}],{i,0,16,3}],Disk[{6,12},0.6`],Disk[{15,7},0.6`]},AspectRatio->Automatic]]

Quelle est la probabilité que la balle traverse le grillage sans le toucher?
a)	Calculez la probabilité exacte μ de cet événement.
	Indication: utilisez une méthode géométrique).
b)	Estimez la fréquence de succès m au moyen d'une simulation.

Exercice 2 - 2

A l'aide d'une méthode de Monte-Carlo, estimez le volume d'une boule de rayon 1, centrée à l'origine:
	B = { (x, y, z) | (x^2)  + (y^2) + (z^2) ≤  1}.

Commentaire:
(x^2)  + (y^2) + (z^2) représente le carré de la distance entre le point P(x, y, z) et l'origine O(0, 0, 0):
			(OP^2) = (x^2)  + (y^2) + (z^2) .
B est donc l'ensemble des points P tels que (OP^2) ≤ 1.

Exercice 2 - 3	(Permutations sans point fixe)

k hommes, avec leurs k épouses respectives, participent à une soirée dansante. L'animateur de la soirée forme des couples au hasard.
Estimez la probabilité qu'aucun homme ne danse avec son épouse.

Indications:
1°	Numérotons les couples et représentons la situation par une liste v = {3, 47, ..., 	14}
	qui signifie, dans le cas particulier,
		l'homme numéro 1 danse avec la femme numéro 3,
		l'homme numéro 2 danse avec la femme numéro 47,
		...,
		l'homme numéro k danse avec la femme numéro 14.

2°	Initialement, la liste est (Subscript[v, 0])= {1, 2, ..., k}.
	Effectuez une permutation au hasard selon la méthode décrite dans  l'exercice 1-2 pour battre les cartes.

3°	Une épreuve est ratée s'il existe un numéro j tel que v[[j]] = j c'est-à-dire 
	s'il existe une coïncidence entre les listes v et (Subscript[v, 0]). En Mathematica, le test peut s'écrire 
		If[MemberQ[v - (Subscript[v, 0]), 0], 0, 1]

4°	Il faut organiser une série de n épreuves.
	Pour de grandes valeurs de k, l'espérance mathématique μ tend vers (1/e) ≈ 0.368

Exercice 2 - 4

Combien de fois faut-il lancer un dé pour qu'apparaisse l'événement "tirer un six" ?

Exercice 2 - 5

Combien de fois faut-il lancer trois dés pour que la somme des 3 dés dépasse 15 ?

Indications: structurez votre programme en sous-programmes, chaque sous-programme correspondant à une fonction:
 *	une expérience consiste à lancer 3 dés et à tester si leur somme dépasse 15;
 *	la variable aléatoire est le nombre d'expériences nécessaires pour obtenir un expérience réussie;
 *	une simulation consiste à évaluer un grand nombre de fois la variable aléatoire.

Exercice 2 - 6

On lance le dé jusqu'à ce que l'on obtienne deux issues consécutives égales puis on considère la somme des points obtenus.
Par exemple : 1 + 5 + 3 + 5 + 4 + 4  = 22.
Estimez la moyenne et l'écart-type des valeurs ainsi obtenues.
Faites un histogramme des résultats obtenus.