RICM4: Probabilité et Simulation

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.

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).

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)\).

• Q1:
  • Vous vous êtes trompé, votre fonction fun1 renvoie 0 (match nul) ou 1 (gagné) mais jamais -1 (perdu). On souhaitait l'espérance du gain, pas la probabilité de gagner. Du coup, vous passez à coté d'un certain nombre d'aspects intéressants de ce DM.
  • D'autre part, avec une telle façon de programmer, qui mélange le codage des probabilités et les règles du jeu dans des tests imbriqués, ça donne un code très difficile à relire et à vérifier. Du coup, vous êtes obligés de coder par copier/coller.
  • Vous en déduisez que la stratégie optimale est obtenue pour x=.8 et y=.1. Vous auriez pu vous douter qu'il fallait aller un peu plus loin et que l'optimal était en fait obtenu pour x=1 et y=0. Et oui, "B doit jouer pierre autant que possible" signifie que B doit jouer systématiquement Pierre…
  • D'autant plus qu'après, vous calculez correctement la probabilité de gagner comme étant 1/4+x/4 mais vous vous obstinez à me dire que le maximum est obtenu pour x=.9 alors que x est dans [0,1].
• Q2:
  • Bonne idée d'inverser les probabilités de jeu pour ne pas jouer de façon trop systématique. Mais pourquoi est-ce que ça marche ? Et est-ce que ça marche bien ? Peut-on vraiment conclure quoi que ce soit avec 20 parties ?…

• Q1: C'est quoi cette histoire de manches ? Je n'ai jamais parlé de manches. Sinon, le code est propre (commenté, variables bien nommées) mais vous mélangez les règles du jeu avec le codage des probas, ce qui oblige à du copier coller. :(
  1. Avec deux joueurs identiques, le joueur B a gagné 17.5 (ce qui est quand même "assez" différent de 0) et vous concluez que personne n'est avantagé ? Que personne n'est avantagé quand c'est symétrique est une évidence, pas une déduction de votre observation…
  2. Vous obtenez 0.6, ce qui, vous dites est très peu différent de ce que vous avez obtenu en question 1, c'est à dire 17.5… Quelle est votre échelle de comparaison ? Ça veut dire quoi très différent pour vous ? En fait, comme vous ne normalisez pas en divisant par N, vous ne calculez pas un estimateur de l'espérance et il est difficile de savoir à quoi s'attendre. Et non, vous ne pouvez clairement pas en déduire que "la différence de probabilité n’est pas assez significative pour changer drastiquement l’espérance de B". Vous êtes juste tombé sur une autre configuration où l'espérance est nulle….
  3. Ce n'est pas une simulation mais une calcul exact… D'autre part, vous vous êtes trombés dans le calcul de pb car vous avez oublié de normaliser x et y…
  4. Par chance, vous tombez sur le bon résultat mais le calcul est complètement faux.
  5. Le calcul est tout à fait correct et tout est bien rédigé.
• Q2:
  • L'idée de toujours jouer aléatoire mais d'adapter "continuement" les probabilités est intéressante
  • On ne peut pas gagner contre quelqu'un qui jouerai avec équi-probabilité. Par contre, comme vous l'avez remarqué, quelqu'un qui jouerai cycliquement aurait les bonnes fréquences mais serait facile à batter. La seule solution pour gérer ce genre de situation est de considérer les paires de coups, d'apprendre les paires effectivement jouées, et de jouer contre ces paires. Votre code est trop compliqué (mélange de la génération des coups du joueur A avec votre propre stratégie) pour le coder simplement…

Le code est très propre, l'analyse est nickel, vous vous êtes posé de bonnes questions en terme de stratégie et vous êtes allez bien plus loin que ce qui était demandé. Seule remarque, votre régression est très bizarre. La bonne régression aurait plutôt été lm(esp_B ~ dataP + dataF + dataC, data=data_frame).

• Code propre, dans le sens où vous séparez bien en petites fonctions bien identifiées. Dommag que vous utilisiez si peu l'expressivité de R par contre (par exemple quand vous testez chaque valeur d'un tableau avec une boucle).
• Q1:
  • Attention, vous confondez la probabilité de gagner et l'espérance du gain.
• Q2: très bien expliqué.

• Q1:
  • Un code propre (plusieurs fonctions bien expliquées) mais un peu compliqué (par exemple passage de plusieurs valeurs dans un tableau dans pouvoir indiquer clairement leur sémantique alors que ces valeurs sont liées les unes aux autres et qu'on pourrait tout recalculer à partir de deux valeurs au lieu de cinq).
  • Bonne idée d'avoir calculer le gain moyen sur plusieurs séries de parties pour bien illustrer la variabilité de la chose.
  • Pour votre troisième plot, il aurait mieux valu mettre la proba de jouer pierre ou de jouer feuille sur l'axe des x que l'index.
  • Bonnes analyses dans l'ensemble.
• Q2:
  • Je ne suis pas d'accord avec vous quand vous écrivez: "De plus, il n’est pas très efficace contre un joueur biaisé.". C'est au contraire l'algo le plus efficace contre un joueur biaisé. Vous vous êtes trompé dans votre implémentation (dans votre if(jB == ((jA[i]+1) %% 3)), vous n'avez pas considéré le match nul), l'espérance de gain ne devrait pas être nulle mais au contraire maximale.
  • Vous auriez mieux fait de faire une fonction à part pour calculer les victoires.
  • Bot0 ("déterministe"): très bien
  • Bot1 (probas "inversées")
  • Bot2 (probas croissantes). Argh le mélange de la détermination de qui gagne avec comment vous mettez vos probas à jour complique bien les choses. Ce bot devrait rapidement avoir la même stratégie que le Bot0 car les probas vont monter de .1 en .1 vers un stratégie pure (1,0,0).
  • Bonne idée d'avoir cherché à apprendre des séquences de coups.

• Vous aimez bien le copier coller, hein ? :) Je ne peux que vous souhaiter de faire souvent des erreurs à la con au tout début et de ne vous en rendre compte que tardivement afin d'apprendre à coder autrement. ;)
• Et afficher plein plein de trucs à l'écran, vous aimez bien aussi, visiblement. :) Il faut sauvegarder ces valeurs et faire des représentations graphiques pour bien montrer les tendances et la variabilité.
• Bonne remarque sur l'humain qui alternerait strictement Ciseau et Pierre.
• L'amélioration que vous proposez consiste à utiliser une estimation glissante des probabilités. Pourquoi pas mais ça ne résoud aucun des problèmes que vous aviez identifié (en particulier, le cas d'un humain qui alterne strictement entre Ciseau et Pierre).

• Q1: Très bien, mais dans votre code, il n'y a aucune simulation, que des calculs exacts de probabilités et d'espérance…
• Q2:
  • Vous proposez un algo qui adapte vos probabilités en fonction du dernier coup joué. C'est une drôle d'idée mais pourquoi pas? Dommage que vous n'ayez visualisé que vb et pas aussi pb pour vérifier/réaliser que vous aviez évolué vers une stratégie pure. D'ailleurs, il me semble que rien n'empèche dans votre code que les p[i] deviennent négatifs… :( Ah et pb est une variable globale ?!? Berk! Dangereux, ça.
  • Stratégie 2 (stratégie "inverse"): vous aimez bien le copier coller, hein ? :) Une fonction commune pour déterminer les gains aurait été adaptée.
• Q3: extension
  • Il y a une erreur dans votre génération (pour Lézard, vous avez oublié p[3] dans votre somme). Vous devriez regarder la fonction sample, ou bien la fonction cumsum.
  • J'ai comme l'impression que quelque chose ne marche pas dans votre fonction, non ? Probablement ce que je vous disais avant, certaines probabilités vont devenir négatives et votre code va faire n'importe quoi.

Quelques remarques mineures sur le code puisque vous y avez été particulièrement attentifs.

• pour votre fonction buildBounds, plutôt que de faire un stop aussi strictement, la véritable chose à tester, c'est si aucune valeur n'est strictement négative. Si ce n'est pas le cas, vous pouvez facilement normaliser en divisant tout par la somme. Ensuite, les bounds s'obtiennent bien avec la fonction cumsum:

    buildBounds = function(pv) {
        if(sum(pv<0)>0 || sum(pv)<=0) { stop("Invalid probability vector"); }
        pv = pv / sum(pv);
        cumsum(pv);
    }
    buildBounds(c(1/4,1/2,1/2))
    buildBounds(c(1/4,1/4,1/4))
    buildBounds(c(-1/4,1/2,1/2))
    buildBounds(c(0,0,0))
  

  [1] 0.2 0.6 1.0
  [1] 0.3333333 0.6666667 1.0000000
  Error in buildBounds(c(-1/4, 1/2, 1/2)) : Invalid probability vector
  Error in buildBounds(c(0, 0, 0)) : Invalid probability vector
  

• Pour votre fonction fromR01toPFC, une version un peu plus dans l'esprit de R:

    fromR01toPFC = function(x, bounds, set=c("P","F","C")) {
        if(x<0 || x>1) { stop("Invalid probability") }
        if(length(bounds)!=length(set)) { stop("Invalid state space") }
        range=1:length(bounds);
        set[min(range[x<bounds])]
    }
    fromR01toPFC(.5,buildBounds(c(1/4,1/2,1/2)))
    fromR01toPFC(1.5,buildBounds(c(1/4,1/2,1/2)))
    fromR01toPFC(.5,buildBounds(c(1/2,1/2)))
  

  [1] "F"
  Error in fromR01toPFC(1.5, buildBounds(c(1/4, 1/2, 1/2))) : 
    Invalid probability
  Error in fromR01toPFC(0.5, buildBounds(c(1/2, 1/2))) : 
    Invalid state space
  

  Sinon, ne pas hésiter à utiliser la fonction sample qui est là pur ça…

• Très bonne utilisation de la fonction sapply et de la fonction mapply. Merci d'avoir indiqué les liens stackoverflow utilisés.
• Petite suggestion: ne pas faire une fonction (allPlotsForXYin01) qui simule et affiche. Faire plutôt une fonction qui simule tout, l'appeler et sauvegarder son résultat et visualiser de plusieurs façons possibles (de façon ad hoc, cette fois-ci) ce résultat.

• Dans la fonction buildDynamicProbabilityVector, plus simple/élégant:

    buildDynamicProbabilityVector = function(playsVector, currentIndex, windowSize) {
        playsVector = playsVector[(max(1,currentIndex - windowSize)):(currentIndex - 1)] # Merci le passage par valeur qui ne détruit par les arguments...
        freqP = length(playsVector[playsVector == "P"]) / length(playsVector);
        freqF = length(playsVector[playsVector == "F"]) / length(playsVector);
        return(c(freqP, freqF, 1 - (freqP + freqF)));
    }
  

Sinon, rien à dire sur le fond, c'est très bien. L'amélioration proposée de votre machine apprenante est "étrange". J'ai l'impression que vous apprenez puis jouez les même probabilités que son adversaire au lieu d'essayer de les contrer. Je pense que c'est pour ça que vous ne voyez aucun gain…

Bonne idée d'avoir tenté une autre stratégie "non aléatoire" (qui fait quand même appel à random) mais effectivement, il n'y a aucune chance pour que ça puisse gagner contre un joeur sans mémoire.

• Super les dessins d'arbres de possibilités pour expliquer ce qu'il se passe.
• Niveau programmation:
  • très bien les break dans les boucles. :) Pour le reste, vous êtes partis sur des codes assez simples au début mais à la fin (partie apprentissage), ce sont des gros paquets de codes, ce qui montre que vous n'avez pas les idées claires.
  • Votre gestion des données avec des tableaux indicés par des listes, par contre, c'est pas terrible.
  • Si vous trouvez (0.9,0,1) comme meilleur vecteur, c'est parceque vous avez itéré sur 1:10 et pas sur =0:
• Pourquoi calculez-vous la médiane à chaque fois ?
• Une moyenne à 3.9, ça ne vous choque pas ? Dommage que vous n'ayez pas trouvé d'où viennent vos problèmes.
• Mais pourquoi donc avoir tenté d'évaluer votre algorithme contre un joueur jouant équiprobable, puisqu'il est alors impossible de gagner?… Vous faites un certain nombre de remarques qui montrent que vous n'avez pas vraiment compris ce qu'il se passait. En particulier la toute dernière.

• Très soigné et bien rédigé.
• Pour votre Algorithm. C'est lent pour deux raisons. D'une part, vous avec une boucle alors que cmp[1]=length(past[past=='R']) aurait été bien plus rapide. Et d'autre part, parcequ'il n'y a pas besoin de tout recalculer depuis le début à chaque coup joué. Il suffirait de mettre cmp à jour. Le véritable problème de cette stratégie, c'est que comme elle est déterministe, un humain va vite la repérer.

• Que les fonctions Ab et Anb sont spécifiques… :( Et s'il n'y a pas de différence entre Ab et Bb, pourquoi ne pas simplement écrire Bb=Ab plutôt que de copier coller ?
• Argh, c'est quoi cette fonction Bxy qui discrétise ? Tout ça pour utiliser la fonction sample à la fin en plus et faire N appels. Faudra réviser ça.
• Argh, c'est quoi ce rown et coln ? La prochaine fois, je demanderai de faire ça avec un pas de 0.01 histoire que vous ne puissiez pas me définir ça à la main.
• Original, la fonction supp qui décrémente x de 1 puis l'incrémente immédiatement après.

• J'attendais des simulations dans la question 1, pas directement des calculs de probabilités.
• À partir du moment où vous décider d'utililiser des chaînes de caractères, mais pourquoi donc les appeler A, B, C, plutôt que P, F, C ?
• Beurk! La fonction play est vraiment inutile mais le codage de la transformation et du gain en global, ça vous gène pas ?

• Ah, enfin un binôme qui a réalisé que ça se calculait bien modulo 0… Vous auriez même pu simplifier encore un peu plus en regardant la différence entre adv et joueur modulo 3.
• Tout le code et les explications sont très clairs et bien écrits. Juste un petit truc étrange, votre première stratégie d'apprentissage obtient un gain (apparemment stabilisé) de 0.185 contre un joueur (1/4,1/4,1/2) alors qu'il devrait avoir un gain optimal de 0.25.

• Votre fonction choix_signe est biaisée. Quel besoin de passer par ce 100 ? Démonstration avec le test du chi2… ;)

    choix_signe <- function(x,y,z){
      rand = sample(100,1, TRUE);
      if(rand<x*100) return(1);
      if(rand<(x+y)*100) return(2);
      return(3);
    }
    N = 1000000;
    X = c();
    for(i in 1:N) { X[i] = choix_signe(1/3,1/3,1/3) }
    chisq.test(table(as.factor(X))) # biaisé...
    chisq.test(table(as.factor(sample(1:3,N,replace=T, p=c(1/3,1/3,1/3))))) # pas biaisé...
  

  
            Chi-squared test for given probabilities
  
    data:  table(as.factor(X))
    X-squared = 208.73, df = 2, p-value < 2.2e-16
  
            Chi-squared test for given probabilities
  
    data:  table(as.factor(sample(1:3, N, replace = T, p = c(1/3, 1/3, 1/3))))
    X-squared = 0.10381, df = 2, p-value = 0.9494
  

• Votre code aurait été un peu plus concis et plus lisible en passant des tableaux de taille 3 plutôt que 3 variables.
• Vous avez oublié d'étudier le cas où le joueur A joue équi-probable.

Très bien. Quelques remarques mineures:x

• Étrange cette première formule avec une somme infinie sur \(i\).
• Attention à votre façon de générer selon une probabilité donnée. Multiplier le nombre d'objets dans le sac n'est pas la bonne façon de faire.
• Le codage par 'P', 'F', 'C', n'aide pas à factoriser le code. L'idéal est d'avoir une fonction de conversion entre cette représentation et une représentation numérique.
• Évitez l'utilisation de variables globales (probA et probB), passez les plutôt en argument à votre fonction Tirages.

• Très bien l'utilisation vectorisée de sample, de sum, de choixB1 == 1, …
• Ceci dit, vous confondez l'espérance du gain et la probabilité de gagner. :( Du coup, vous vous posez des questions étranges et vous passez à côté de l'analyse du jeu. En particulier, vous espériez que la probabilité de gagner soit de 1/3 dans le cas où les deux joueurs osnt biaisés. Il se trouve qu'elle ne l'est pas. Une façon de s'en convaincre, c'est de regarder la situation extrème où deux joueurs jouent avec la même stratégie biaisée (.99,005,.005). Vous voyez bien qu'il y aura quasiment tout le temps des ex-aequo et que, même si chaque joueur va gagner aussi souvent qu'il va perdre, la probabilité de gagner n'est absolument pas d'1/3.
• Attention, vous confondez "aléatoire" et "équiprobable": quand vous écrivez "le fait qu’un joueur choisisse de jouer de manière totalement aléatoire" vous vouliez dire "le fait qu'un joueur joue de façon équiprobable". Et votre remarque sur le fait que le joueur serait favorisé parce que la valeur obtenue cette fois-ci est un peu plus élevée que celle d'avant est … fausse. C'est une observation donc vous ne pouvez rien déduire car il les valeurs sont proches et qu'il y a trop de variabilité pour conclure quoi que ce soit.
• Tiens, res = sample(0, 1000, TRUE), intéressante façon d'initialiser un tableau. :) Je vous suggère plutôt rep.int(0,1000). Notez, que j'avais demandé d'avoir x et y qui varient régulièrement, pas qu'ils soient tirés aléatoirement.

• Il faut commenter un peu plus vos résultats.

• Dans pccprob (quel nom explicite et judicieux d'ailleurs… ;), vos tirages et transformations de nombres aléatoiries sont incorrects. Multiplier par 100 ne sert à rien et comme vous tirez entre 1 et 100 (i.e., sur un intervalle de longueur 99), vos renormalisations font n'importe quoi. Démonstration:

    votretirage = function(ppA) {
      j1 = runif(1,1,100);
      if(j1 <= ppA*100) { 
          return("pp");
      } else {
        return("Autre chose")
      }
    }
    X = c()
    for(i in 1:1000000) {
       X[i] = votretirage(.001);
    }
    length(X[X=="pp"]); ### Alors, sur 1 million de tirages, combien de fois avons-nous obtenu "pp" ?
    length(X[X=="Autre chose"]);
  

  [1] 0
  [1] 1000000
  

• Pas besoin de la fonction maxtab quand la fonction max et which.max sont déjà définies.
• Balancer tout le code comme ça sans explications et l'utiliser après, n'est pas très pédagogique. Est-ce vraiment comme ça que vous avez travaillé ?
• Comme un certain nombre d'autres, vous confondez espérance de gain et probabilité de gagner. :( Du coup, vous passez à côté du sujet.

• Il ne faut pas donner le même label à tous vos blocs… C'est pour ça que vous ne pouviez pas exporter.
• Vous confondez probabilité de gagner et espérance de gain.

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! ;)

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…

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