Exercice 1¶

Les questions 1 à 3 se résolvent simplement en calculant des puissances de la matrice de transition :

In [1]:
# Matrice de transition
P_ = [
    # D    J     M      (états)
    [0.9, 0.05, 0.05], # Dormir
    [0.8, 0.2,  0.0],  # Jouer
    [0.7, 0.3,  0.0]   # Manger
]

import numpy as np
import numpy.linalg as li
P = np.matrix(P_) #pour calculer P^n

Le futur ne dépendant que de l'état courant, la probabilité que Gribouille joue à $t=8$ sachant qu'il dort à $t=5$ est la même que celle qu'il joue à $t=3$ sachant qu'il dort à $t=0$. Il suffit de regarder $P^3[1,2]$ (indices à partir de 1).

In [2]:
li.matrix_power(P, 3)[0,1]
Out[2]:
0.07175000000000001
In [3]:
li.matrix_power(P, 12)[2,0]
Out[3]:
0.8839779005525006
In [4]:
li.matrix_power(P, 50)
Out[4]:
matrix([[0.8839779, 0.0718232, 0.0441989],
        [0.8839779, 0.0718232, 0.0441989],
        [0.8839779, 0.0718232, 0.0441989]])

On observe numériquement la distribution stationnaire $w = (0.8839779, 0.0718232, 0.0441989)$. La convergence est en fait très rapide sur cet exemple. On peut aussi calculer $w$ en résolvant l'équation $w P = w$ :

$$\begin{align*} 0.9 w_1 + 0.8 w_2 + 0.7 w_3 &= w_1\\ 0.05 w_1 + 0.2 w_2 + 0.3 w_3 &= w_2\\ 0.05 w_1 &= w_3\\ \end{align*}$$

Le système admet une infinité de solutions car $P$ est de rang 2 : sa dernière colonne se déduit des autres étant donné que la somme par lignes vaut 1. On observe que le choix $w_1 = 0$ mène à $w = 0$ et n'est donc pas envisageable. Fixons alors arbitrairement $w_1 = 1$ : alors $w_3 = 0.05 = \frac{1}{20}$ et $0.8 w_2 = 1 - 0.9 - 0.7 \times 0.05$ donc $w_2 = \frac{0.065}{0.8} = \frac{13}{160}$. On normalise ensuite par $w_1 + w_2 + w_3 = \frac{181}{160}$ pour obtenir :

In [5]:
[160/181, 13/181, 8/181]
Out[5]:
[0.8839779005524862, 0.0718232044198895, 0.04419889502762431]

Pour la dernière question, il faut considérer l'état "manger" comme absorbant. On obtient

$$Q = \begin{pmatrix}0.9 & 0.05\\ 0.8 & 0.2\end{pmatrix} \quad \mbox{et} \quad N = \begin{pmatrix}0.05\\0.0\end{pmatrix}$$
In [6]:
Q = np.matrix([[0.9, 0.05],[0.8, 0.2]])
R = np.matrix([[0.05],[0.0]])
N = li.inv( np.identity(2) - Q )
print(np.sum(N, axis=1))

[[21.25]
 [22.5 ]]
In [7]:
np.identity(2) - Q
Out[7]:
matrix([[ 0.1 , -0.05],
        [-0.8 ,  0.8 ]])

Il faudra donc attendre un peu plus de 21 changements d'états en moyenne pour que Gribouille mange.

On peut aussi faire le calcul à la main dans ce cas simple : $I - Q = \begin{pmatrix}0.1 & -0.05\\ -0.8 & 0.8\end{pmatrix}$ puis $\mbox{det}(I-Q) = 0.08 - 0.04 = 0.04$, et donc $(I-Q)^{-1} = 25 \begin{pmatrix}0.8 & 0.05\\0.8 & 0.1\end{pmatrix} = \begin{pmatrix}20 & 1.25\\20 & 2.5\end{pmatrix}$.

In [8]:
N @ R #colonne de 1 : normal, on finit par manger avec probabilité 1
Out[8]:
matrix([[1.],
        [1.]])

Exercice 2¶

Choisir une stratégie revient à sélectionner une fonction $f$ de $[1, 7]$ dans lui-même, qui à un montant en euros associe l'argent parié. Bien sûr $f(x) \leq x$ : on ne peut pas parier plus que ce qu'on a. On ne définit pas $f(0)$ ou $f(8)$ car dans les états "0" et "8" le jeu est terminé : dans le premier cas on a perdu (ruine) et dans le second on a gagné (sortie de prison).

L'idée la plus simple consiste à toujours miser tout ce qu'on a : $f(x) = x$. Dans ce cas, on passe initialement de 3€ à 0 ou 6€, puis de 6€ à 0 ou 8€ (12 en fait, mais toute somme supérieure ou égale à 8 est ramenée à 8€ puisque 8 suffit). Il n'y a donc que 4 états et la mise en équation est très simple.

In [9]:
# Matrice de transition
P_ = [
    # 3    6    0    8    (états)
    [0.0, 0.4, 0.6, 0.0], # 3
    [0.0, 0.0, 0.6, 0.4], # 6
    [0.0, 0.0, 1.0, 0.0], # 0
    [0.0, 0.0, 0.0, 1.0]  # 8
]
# Remarque : l'ordre dans lequel les états apparaissent n'importe pas,
# mais il est plus facile de regrouper les états absorbants en bas à droite.

P = np.matrix(P_) #pour utiliser numpy

# Matrice fondamentale : N = (I - Q)^{-1})
indices = range(2) #2 = nombre d'états non absorbants (transients)
Q = P[indices,:][:,indices]
N = li.inv( np.identity(2) - Q )

Rappel 1 : np.identity(n) construit une matrice carrée $I$ de taille $n$, contenant des 1 sur la diagonale et des 0 partout ailleurs. Elle vérifie $I X = X I = X$ pour toute matrice carrée $X$ de taille $n$ (d'où le nom "matrice [fonction] identité" : c'est l'élément neutre pour la multiplication).

In [10]:
print(np.identity(3))
print("")
print(np.identity(4))

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]

Rappel 2 : l'inverse d'une matrice $M$, notée $M^{-1}$, est l'unique matrice $A$ vérifiant $A M = M A = I$ (identité).
Attention certaines matrices ne sont pas inversibles.

Calcul pratique : Pour une petite matrice (taille 2 ou 3) on peut appliquer la méthode des cofacteurs :
$M^{-1} = \frac{1}{\mbox{det}(M)} \mbox{t(com}(M))$
Pour une matrice contenant un grand nombre de zéros, on peut chercher à inverser le système $M X = Y$,
l'écrivant sous la forme $A Y = X$ : $A$ est alors l'inverse cherché.
Dans le cas général on peut appliquer l'algorithme du pivot de Gauss.

Ce calcul ne sera pas demandé à l'examen (je vous donnerai directement la matrice fondamentale $N$), mais comme toujours, c'est mieux si vous savez comment obtenir l'inverse d'une matrice à la main.

In [11]:
# Suite de l'exercice : matrice R, puis N x R
R = P[indices,:][:,range(2, P.shape[1])]
absorb_probs = N @ R
absorb_probs[0, 1] #partant de l'état 0 = 3€, probabilité de sortir de prison (état absorbant 1)
Out[11]:
0.16000000000000003

Le calcul devrait donner exactement 0.16, mais on observe une erreur d'arrondi. Explication courte : 0.4 = 4/10 = 2/5 n'est pas représentable de manière exacte en écriture binaire. Explication longue : voir par exemple ce document.

In [12]:
print(0.4 * 0.4) #2/5 * 2/5 : erreur
print(0.5 * 0.5) #1/2 * 1/2 : OK

0.16000000000000003
0.25

Voir aussi https://stackoverflow.com/a/16444778/12660887.

Bref, sur ce premier scénario simple, probabilité de sortir de prison = 0.16.
C'est ce qu'on obtient immédiatement en dessinant l'arbre de probabilités, mais dans un cas plus complexe l'écriture matricielle est préférable.

Remarque : miser 6 lorsqu'on a 6€ est moins intéressant que de miser 2 seulement (puisque tout montant de 8€ au moins suffit à sortir). En cas de défaite on peut alors rejouer à partir de 4€ au lieu d'avoir perdu.

In [13]:
# Modélisation modifiée : on mise 2 dans l'état "6€", et 4 dans l'état "4€" :
P_ = [
    # 3    4    6    0    8    (états)
    [0.0, 0.0, 0.4, 0.6, 0.0], # 3
    [0.0, 0.0, 0.0, 0.6, 0.4], # 4
    [0.0, 0.6, 0.0, 0.0, 0.4], # 6
    [0.0, 0.0, 0.0, 1.0, 0.0], # 0
    [0.0, 0.0, 0.0, 0.0, 1.0]  # 8
]

# Résumé des opérations réalisées, dans une fonction :
def getAbsorbProbs(P):
    nbTransients = P.shape[0] - 2
    indices = range(nbTransients)
    Q = P[indices,:][:,indices]
    N = li.inv( np.identity(nbTransients) - Q )
    R = P[indices,:][:,range(nbTransients, P.shape[1])]
    return N @ R

getAbsorbProbs(np.matrix(P_))[0, 1]
Out[13]:
0.256

C'est mieux. On obtient le même résultat avec un arbre : 0.4 * 0.4 + 0.4 * 0.6 * 0.4

Que dire de la stratégie "la plus prudente" consistant à toujours miser exactement 1€ ?

In [14]:
P_ = [
    # 1    2    3    4    5    6    7    8    0    (états)
    [0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6], # 1
    [0.6, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], # 2
    [0.0, 0.6, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0], # 3
    [0.0, 0.0, 0.6, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0], # 4
    [0.0, 0.0, 0.0, 0.6, 0.0, 0.4, 0.0, 0.0, 0.0], # 5
    [0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.4, 0.0, 0.0], # 6
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.4, 0.0], # 7
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], # 8
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]  # 0
]

# 7 états transients. On part du 3eme = indice 2.
getAbsorbProbs(np.matrix(P_))[2, 0]
Out[14]:
0.09643140364789855

Pas terrible : moins bien que la stratégie naïve "tout miser".

Pour faire le calcul à l'aide d'un arbre, il faudrait envisager "tous" les chemins. Le souci est qu'ils sont ici en nombre infini : 1 -> 2 -> ... -> 8 en direct, mais aussi 1 -> 2 -> 1 -> 2 -> 3 etc. On ne s'en sort pas : le calcul matriciel est nécessaire.

Intuitivement on sent que moins il y a de tours mieux c'est, car on a toujours plus de chances de perdre (60%) que de gagner (40%). Cette stratégie "prudente" nécessite de gagner au moins 7 fois (1 -> 2, 2 -> 3 etc), et n'est donc pas du tout prudente. Qu'en est-il alors d'une stratégie intermédiaire ? Dans l'état "$n$€", disons qu'on mise min(max($1, n-1$),$8-n$), pour qu'une défaite ne soit pas nécessairement fatale et qu'on ne dépasse pas 8 :

In [15]:
# Fonction calculant P étant donné le paramètre de gain p (0.4 jusqu'ici), et la stratégie S
def buildP(p, S):
    P = np.zeros((9, 9)) #9 états = cas le plus général - sur cet exercice
    for state in range(1, 8):
        mise = S(state)
        P[state, state - mise] = 1 - p #on perd
        P[state, state + mise] = p #on gagne
    # On complète avec les états absorbants
    P[0, 0] = 1
    P[8, 8] = 1
    # réagencement de P : état 0 à la fin https://stackoverflow.com/a/34443002/12660887
    # L'état 7 sera en premier. Pas très grave : 7 1 2 3 4 5 6 donc. "3€" en 4eme position
    P[[0, 7], :] = P[[7, 0], :]
    P[:, [0, 7]] = P[:, [7, 0]]
    return (P)

def strategy(n):
    return min( max(1, n-1), 8-n )

P = buildP(0.4, strategy)
getAbsorbProbs(P)[3, 1]
Out[15]:
0.21408450704225357

Pas mieux que la stratégie "miser tout". Resterait à prouver qu'elle est optimale... (Je pense que oui).

In [16]:
# Tests : pour p = 0.5, on doit trouver 0.5 quelle que soit la stratégie, à condition de partir de 4€.
#         les probabilités d'absorption ne doivent alors pas dépendre de la stratégie.
def S1(n):
    return 1

def S2(n):
    return min(n, 8-n) #optimal, semble-t-il

print( getAbsorbProbs( buildP(0.5, S1) )[:, 1] )
print( getAbsorbProbs( buildP(0.5, S2) )[:, 1] )
print( getAbsorbProbs( buildP(0.5, strategy) )[:, 1] )

[0.875 0.125 0.25  0.375 0.5   0.625 0.75 ]
[0.875 0.125 0.25  0.375 0.5   0.625 0.75 ]
[0.875 0.125 0.25  0.375 0.5   0.625 0.75 ]
In [17]:
# Pour p = 0.6 on devrait retrouver exactement les résultats précédents
# en partant cette fois de 5€ et en cherchant à perdre (situation symétrique) :
print( getAbsorbProbs( buildP(0.6, S1) )[5, 0] )
print( getAbsorbProbs( buildP(0.6, S2) )[5, 0] )
print( getAbsorbProbs( buildP(0.6, strategy) )[5, 0] ) #erreur d'arrondi, le retour ? TODO...

0.09643140364789848
0.256
0.21694915254237299