RICM4: Probabilité et Simulation

Table of Contents

<div id="site-map"><h2>Sitemap</h2><div id="text-site-map"> <?php include ("lib/PHPLIB.php"); // taken from PHPLib include ("lib/layersmenu-common.inc.php"); include ("lib/treemenu.inc.php"); if ( @ is_file ("../../menu-structure-file.txt") ) { $mid = new TreeMenu(); $mid->setMenuStructureFile("../../menu-structure-file.txt"); $mid->setPrependedUrl("../../"); $mid->parseStructureForMenu("treemenu1"); print $mid->newTreeMenu("treemenu1"); } ?> </div></div>

Informations Générales

Florence Perronnin est chargée des cours. Arnaud Legrand et Florence Perronnin s'occupent des TDs.

Le planning avec les salles de cours est disponible ici.

Voici les annales des quicks des années précédentes.

Programme du cours

Machines perso obligatoires en TD.

  • 6 Septembre 2017 (8:00 - 9:30): Cours (FP)
    • Documents: Présentation des objectifs du cours. Activité autour du paradoxe des aniversaires:
    • À faire pour vendredi:
      1. Installer Rstudio sur votre portable. Pour cela, voir la section sur Rstudio un peu plus bas.
      2. Préparez le TD en suivant ce tutorial. Il y en a pour 40-45 minutes de lecture mais c'est très progressif et il y a 10 minutes de petits exercices pour vérifier que vous avez bien compris.
  • 8 Septembre 2017 (8:00 - 11:15): TD (AL) Au programme:
    • Évaluation de votre autonomie (Rstudio installé, tuto essayé, …)
    • Prise en main de R et de RStudio, premiers documents avec Rmarkdown.
    • Poursuite du document sur les anniversaires afin de simuler ce scénario et étudier l'impact de la taille de la classe sur la fréquence de "colision". Voici le résultat final (sources)
  • 15 Septembre 2017 (9:45 - 13:00): TD (FP)
    • Préparation du TD:
      • Préparez le TD en suivant ce tutorial. Il y en a pour 40-45 minutes de lecture. Il présente les principaux types et structures de données en R et permettre de démystifier un certain nombre de choses. C'est très progressif mais ne brulez pas les étapes pour bien tout comprendre. Les exercices sont faciles.
      • Objectifs de la séance: Modéliser et simuler des situations simples. Faire le lien entre la modélisation mathématique et le code de la simulation.
  • 25 septembre 2017 (8:00-9:30): Cours (FP)
    • Notions abordées:
      • Probabilités conditionnelles, formule de Bayes.
    • Documents:
      • Une belle illustration et application de ces concepts aux maladies rares:
      • Correction avec l'approche 1 «statistique » (celle faite par l’article) et 1 « modélisation » avec définition plus systématique d’événements et manipulation de probabilités conditionnelles, loi de Bayes et loi des probabilités totales. Les deux approches sont en fait basées sur le même raisonnement. La 3e méthode c’est la simulation mais à vous de le coder. Voici notre correction.
      • Pour finir, un exemple du paradoxe de Yule-Simpson, sur les calculs rénaux (ambiance médicale !! ☺)
    • À faire pour la semaine prochaine: Modéliser l'exemple sur les calculs rénaux avec des probas conditionnelles
  • 29 Septembre 2017 (9:45 - 13:00): TD (AL) Correction des points pas clairs de l'évaluation de rentrée. Monty hall (Are Birds Smarter Than Mathematicians?) si le temps le permet.
    • Quelques petits exemples de code:
      • Probabilité d'avoir observé au moins une fois un 6 en lançant un dé 4 fois de suite:

              N=10000;
              R1=sample(1:6,N,replace = T)
              R2=sample(1:6,N,replace = T)
              R3=sample(1:6,N,replace = T)
              R4=sample(1:6,N,replace = T)
              sum(ifelse((R1==6) | (R2==6) |(R3==6) | (R4==6),10, -10))
        
      • Probabilité d'avoir un chemin dans le circuit:

              N=1000000
              a = sample(c(FALSE,TRUE), N, replace = T, prob = c(0.1, 0.9))
              # hist(ifelse(a,1,0))
              b = sample(c(FALSE,TRUE), N, replace = T, prob = c(0.1, 0.9))
              c = sample(c(FALSE,TRUE), N, replace = T, prob = c(0.1, 0.9))
              d = sample(c(FALSE,TRUE), N, replace = T, prob = c(0.4, 0.6))
              e = sample(c(FALSE,TRUE), N, replace = T, prob = c(0.25, 0.75))
              f = sample(c(FALSE,TRUE), N, replace = T, prob = c(0.25, 0.75))
              sum((a & b & c) | (d) | (e & f))/N
        
      • Monty Hall. Three different versions:
        1. Fichue fonction sample (lisez la doc et regardez la sémantique de la variable x).

                      resample <- function(x, ...) x[sample.int(length(x), ...)]
                      N = 1000 # Number of trials
                      S1 = 0   # Number of cases where the first strategy leads to a win
                      S2 = 0   # Number of cases where the second strategy leads to a win
                      doors = 1:3
                      for(trial in 1:N) {
                        win = sample(doors, 1)
                        strat1 = sample(doors,1)
                        reveal =  resample(doors[-c(strat1,win)],1) # Damn, sample does not work here and the help page of sample explains why. The resample function is directly stolen from there.
                        # Bon, mais ce qui compte ici, c'est qu'on voit bien là qu'avec la stratégie 2, on a une chance sur deux alors qu'on n'a qu'une chance sur trois avec la stratégie 1
                        strat2 = doors[-c(strat1,reveal)]
                        S1 = S1 + (strat1 == win)
                        S2 = S2 + (strat2 == win)
                      }
                      S1/N
                      S2/N
          
        2. Si on code les noms des portes avec des lettres, le problème ne se pose plus et ça s'écrirait:

                      N = 1000 # Number of trials
                      S1 = 0   # Number of cases where the first strategy leads to a win
                      S2 = 0   # Number of cases where the second strategy leads to a win
                      doors = c("A","B","C")
                      for(trial in 1:N) {
                        win = sample(doors, 1)
                        strat1 = sample(doors,1)
                        reveal =  sample(doors[!(doors %in% c(strat1,win))],1) 
                        strat2 = doors[!(doors%in%c(strat1,reveal))]
                        S1 = S1 + (strat1 == win)
                        S2 = S2 + (strat2 == win)
                      }
                      S1/N
                      S2/N
          
        3. Bon, enfin, pourrait-on faire une version vraiment vectorielle sans cette boucle for ? Oui, mais il faut vraiment se tordre la cervelle et utiliser apply. Berk, oubliez, c'est inutilement moche.
  • 3 Octobre 2017 (8:00-9:30): Cours (AL): le hasard, c'est quoi au fait ? Comment en générer ? Voir https://tinyurl.com/RICM4PS17.
  • 6 Octobre 2017 (9:45 - 13:00): TD: Pas de TD (dev. P.)
  • 10 Octobre 2017 (8:00-9:30): Cours (FP): Variable aléatoire discrète (uniforme, bernouilli, binomiale, géométrique).
  • 13 Octobre 2017 (9:45 - 13:00): TD: Pas de TD (dev. P.)
  • 17 Octobre 2017 (8:00-9:30): Cours (FP): Notion de générateur pseudo-aléatoire, de graine, de période. Passage d'une loi uniforme (sur \([0,1]\) ou sur \(\{0,1\}\)) vers une autre loi uniforme (sur \([a,b]\) ou sur \(\{0,1,...,n\}\)) avec rejet si besoin.
  • 20 Octobre 2017 (9:45 - 13:00): TD (FP): Génération de nombres aléatoires et Simulation de lois uniformes.
  • 24 Octobre 2017 (8:00-9:30): Cours (FP): Générateurs de loi uniforme et de lois discrètes (tabulation, fonction de la fonction de répartition, rejet).
  • 20 Octobre 2017 (9:45 - 13:00): TD (AL): Génération de structures
    • Énumération des permutations

      Perm_list = list()
      enum_permutation = function(perm,remain) {
          if(length(remain)==0) {
              print(perm);
              Perm_list[[length(Perm_list)+1]] <<- perm 
          }
          for(i in remain) {
              enum_permutation(c(perm,i),remain[remain!=i])
          }
      }
      N=3
      enum_permutation(c(),1:N)
      
      [1] 1 2 3
      [1] 1 3 2
      [1] 2 1 3
      [1] 2 3 1
      [1] 3 1 2
      [1] 3 2 1
      
    • Génération des permutations (4 méthodes)

      # Method 0 :)
      sample(1:N) 
      # Method 1
      Perm_list[[sample(1:(length(Perm_list)),1)]]
      # Method 2
      gen_permutation = function(val, N, perm) {
          possibilities = 1:N;
          possibilities = possibilities[!(possibilities %in% perm)];
          if(length(perm)==(N-1)) { perm[N] = possibilities[1]; return(perm); }
          num_perm = factorial(N-length(perm)-1);
          idx = floor(val / num_perm) + 1;
          perm[length(perm)+1] <- possibilities[idx];
          return (gen_permutation(val %% num_perm, N, perm));
      }
      gen_permutation(floor(runif(1,max=factorial(N))), N, c());
      # Method 3
      perm = 1:N
      for(i in 1:(N-1)) {
          idx = sample(i:N, 1)
          tmp = perm[i];
          perm[i] = perm[idx];
          perm[idx] = tmp;
      }
      perm
      # Method 3 bis
      set = 1:N
      perm = c();
      for(i in 1:(N-1)) {
          perm[i] = sample(set[!set %in% perm], 1);
      }
      perm[N] = set[!set %in% perm];
      perm
      #Method 4 (Knuth)
      perm = 1:N
      for(k in 1:N) {
          i = sample(i:N, 1)
          j = sample(i:N, 1)
          tmp = perm[i];
          perm[i] = perm[j];
          perm[j] = tmp;
      }
      perm
      
      [1] 3 1 2
      [1] 1 2 3
      [1] 2 3 1
      [1] 2 3 1
      [1] 1 3 2
      [1] 1 2 3
      
    • Énumération des sous-ensembles à \(k\) éléments
    • Génération des sous-ensembles à \(k\) éléments (5 méthodes)
    • Génération d'arbres binaires à \(n\) sommets.

          tree_t create_tree(int n) {
            tree_t t1, t2;
            if (n==0) 
              return empty_tree;
            else {
              q=floor((double)(n-1)*rand()/RAND_MAX);
              a1=create_tree(q);
              a2=create_tree(n-q-1);
              return create_tree(a1,a2);
            }
          }
      
      • Tous les arbres à \(n\) noeuds internes peuvent-ils être créés ?
      • Le générateur est il de loi uniforme ?
      • Soit \(X_n\) le nombre de noeuds internes sur le chemin le plus à gauche dans l’arbre généré ayant \(n\) noeuds.
        • Calculer la loi de \(X_0\) , \(X_1\) , \(X_2\) et \(X_3\). Vérifier que \(E[X_3] = 1 + \frac{1}{3}(E[X_0] + E[X_1] + E[X_2] )\).
        • Montrer que \(E[X_n] = 1 + \frac{1}{n}\sum_{i=0}^{n} E[X_i]\). En déduire une expression de \(E[X_n]\) en fonction de \(E[X_{n−1}]\). Calculer \(E[X_n]\) et donner son équivalent.
      • Énumérer les arbres à \(n\) noeuds et en déduire comment débiaiser l'algorithme de génération précédent.
  • 17 Novembre 2017 (9:45 - 13:00): TD (AL) Génération selon la loi normale
    • Rappel de ce que c'est qu'une variable aléatoire continue, fonction de répartition, une densité de probabilité (http://fr.wikipedia.org/wiki/Fonction_de_répartition), l'inverse de la fonction de répartition…
    • Une "correction" du TP: http://rpubs.com/alegrand/10381.
    • Voici aussi un lien qui explique comment superposer une densité sur un histogramme.
    • Un lien avec des slides sur ggplot2: https://github.com/alegrand/SMPE/blob/master/lectures/lecture_R_crash_course.pdf
    • Au final, on a vu trois méthode pour générer une loi normale:
      • En utilisant l'inverse de la fonction de répartition (qnorm(runif())). Mais cette méthode demande de savoir inverser la fonction de répartition, ce que l'on ne sait faire que via des approximations numériques;
      • En sommant d'autres variables aléatoires (par exemple une douzaine de lois uniformes sur \([0,1]\)). Cette méthode a le mérite d'être simple mais demande beaucoup d'appels à la fonction random. Elle n'est pas parfaite mais donne une bonne approximation;
      • Avec la méthode de Box-Muller: en générant un angle uniforme dans \([0,2\pi]\) et un rayon au carré selon une loi exponentielle, on obtient un point dont chacune des coordonnées sont indépendantes et générées selon une loi normale! C'est une méthode très très élégante, qui produit deux nombre pour deux appels à random, mais qui utilise des fonctions mathématiques un peu coûteuses (log, cos, sin, pi).

DM

DM1: Pierre Feuille Ciseaux

Règles du jeu:

Dans ce DM on s'intéresse au jeu de Pierre Papier Ciseaux et on cherche à savoir quelle est la bonne stratégie à avoir. On se concentre sur la version (non violente!) dite "à somme nulle" où le gagnant empoche 1 alors que le perdant perd 1. En cas d'égalité, personne ne perd ni ne gagne rien.

  • Rendu
    • À faire en binôme. Je vous suggère de vous mélanger autant que possible entre les deux groupes de TD.
    • Vous rendrez votre devoir http://rpubs.com à l’aide de Rstudio et enverrez l'url à arnaud.legrand@imag.fr avant le 23 novembre à minuit en indiquant la chaîne [RICM4-PS] dans le sujet du message. La deadline est stricte car nous corrigerons en TD le 24 novembre.
    • Vous indiquerez bien vos noms et prénoms dans l'entête du document Rmd afin qu'il n'y ait pas d'ambiguïté
    • Ce qui m'intéresse n'est pas tant le code que ce que vous concluez de vos simulations.
  • Plagiat
    • Les URLs des DMs sont publiques. Si je vois trop de ressemblances entre deux rendus, ça va vraiment m'agacer. Vous avez le droit de discuter d'un binôme à l'autre si vous rencontrez des difficultés mais il faut alors l'indiquer et l'expliquer clairement dans votre DM.
    • De même si vous vous inspirez de sources (wikipedia, stackoverflow, …), il convient de l'indiquer clairement.

Q1

Les poissons chirurgiens, à la différence des humains, sont capable de parfaitement jouer "au hasard" sans tenir compte de leur historique. Voyons si ce trouble répandu de la mémoire à court terme leur confère un avantage au jeu de Pierre Feuille Ciseaux.

Un joueur sans mémoire \(A\) est parfaitement défini par son vecteur de probabilité \(P^{(A)}\) indiquant avec quelle probabilité il joue chaque symbole. Un joueur \(A\) ayant un \(P^{(A)}=(1/4,1/4,1/2)\) jouera Pierre et Papier avec probabilité \(1/4\) et Ciseau avec probabilité 1/2.

On cherche d'abord à comprendre ce qu'il se passe quand deux joueurs sans mémoire s'affrontent.

Vous réaliserez une simulation vous permettant d'estimer les espérances de gains (en faisant suffisamment de parties, le "suffisamment" étant à estimer par vos soins mais soyez raisonnables).

  • Le joueur biaisé

    On suppose que le joueur \(A\) joue avec probabilités \(P^{(A)}=(1/4,1/4,1/2)\).

    1. Estimer l'espérance de gain du joueur \(B\) contre le joueur \(A\) lorsque \(P^{(B)}=(1/4,1/4,1/2)\).
    2. Même question lorsque \(P^{(B)}=(1/3,1/3,1/3)\).
    3. Estimer l'espérance de gain du joueur \(B\) tel que \(P^{(B)}=(x,y,1-x-y)\) lorsque \(x\) et \(y\) varient entre \(0\) et \(1\) par pas de \(0.1\) ?
    4. Déduisez-en la stratégie optimale pour le joueur \(B\).
    5. Retrouver et vérifier les résultats précédents à l'aide d'un calcul exact de ces espérances (i.e, en obtenant une formule dépendant de \(x\) et de \(y\)).
  • Le joueur non biaisé

    Mêmes questions que précédemment mais lorsque le joueur \(A\) est non biaisé, i.e., lorsque \(P^{(A)}=(1/3,1/3,1/3)\).

    Un joueur non biaisé peut-il perdre ? Peut-il gagner ?

Q2: Apprentissage

Notre objectif est maintenant de créer un programme qui arrive à gagner contre un joueur humain. Nous savons que les humains jouent mal au hasard et peuvent en particulier avoir du mal à maintenir l'uniformité.

  1. Proposez un algorithme qui "apprenne" les fréquences de jeu d'un adversaire et s'y adapte optimalement en permanence.
  2. Évaluer le gain de votre algorithme contre un joueur biaisé qui jouerait avec \(P^{(A)}=(1/4,1/4,1/2)\). N'hésitez pas à visualiser l'évolution du jeu au cours du temps.
  3. Votre algorithme a-t-il une chance de battre un humain pas trop bête? Comment pourriez-vous simplement l'améliorer?
  4. Vous évaluerez la performance de votre nouvel algorithme contre un joueur biaisé qui jouerait avec \(P^{(A)}=(1/4,1/4,1/2)\).

Q3 (bonus au choix):

Technical references

Installing R and Rstudio

Virtual Machine

Les anneés précédentes, je préparais une VM avec une debian récente, toute bien installée mais ça n'a pas eu l'air d'aider tant que ça les gens donc débrouillez vous, et installez le en natif, ça sera formateur! ;)

Mac OSX

A few years ago, a nice RICM4 student, Remi Gattaz, has taken the time to explain how to install a bunch of useful stuff. Here it is. In particular he gave many tips for MacOSX…

Linux

Here is how to proceed on debian-based distributions:

sudo apt-get install r-base r-cran-ggplot2 r-cran-reshape

Make sure you have a recent (>= 3.2.0) version or R. For example, here is what I have on my machine:

R --version
R version 3.2.0 (2015-04-16) -- "Full of Ingredients"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under the terms of the
GNU General Public License versions 2 or 3.
For more information about these matters see
http://www.gnu.org/licenses/.

If it's not the case, it may be because you're running debian stable or a LTD ubuntu. In such case, you may want to include testing packages… Ask your local linux guru or run a VM (see previous section) if you're affraid to break your OS. For the braves, let's keep going!

Rstudio and knitr are unfortunately not packaged within debian so the easiest is to download the corresponding debian package on the Rstudio webpage and then to install it manually (depending on when you do this, you can obviously change the version number).

wget https://download1.rstudio.org/rstudio-1.0.153-amd64.deb
sudo dpkg -i rstudio-1.0.153-amd64.deb
sudo apt-get -f install # to fix possibly missing dependencies

You will also need to install knitr. To this end, you should simply run R (or Rstudio) and use the following command.

install.packages("knitr")

If r-cran-ggplot2 or r-cran-reshape could not be installed for some reason, you can also install it through R by doing:

install.packages("ggplot2")
install.packages("reshape")

You may have trouble when installing some R packages. If so, try to install these ones:

sudo apt-get install libcurl4-openssl-dev libssl-dev

Bibliographie