Problème du bandit à $K$ bras¶

*Master Class Lycéennes - 3 avril 2025*

Margaux Brégère

Une statisticienne rentre dans un casino et se retrouve face à plusieurs machines à sous. Il existe des machines qui, lorsqu'elles sont jouées, rapportent souvent, et des machines sur lesquelles il est beaucoup plus rare de gagner. La statisticienne dispose d'un certain nombre de coups à jouer et a pour objectif de maximiser ses gains. Quelles sont les bonnes stratégies ? Nous formaliserons ce problème mathématique du bandit à K bras. Nous aborderons ensuite certaines stratégies optimales pour le résoudre. Ces dernières sont, par exemple, utilisées dans les algorithmes de recommandation ou encore pour optimiser les tests cliniques.

Simulation du problème du bandit à $K$ bras¶

Formalisation mathématique :¶

Il y a $K \in \mathbb{N}^\star$ machines à sous dans le casino. La statisticienne dispose de $T\in \mathbb{N}^\star$ coups et a pour objectif de maximiser ses gains.

Pour $t = 1,\dots, T$

$\qquad \bullet$ Elle choisit une machine à sous $I_t \in \{1,\dots, K\}$

$\qquad \bullet$ Elle reçoit le gain $g_t$ issu de la machine $I_t$

Mathématiquement, une loi de probabilité est associée à chaque machine. Ici, par exemple, le gain de chaque machine à sous $k$ suit une loi de probabilité de Bernoulli de paramètre $p_k$. Ainsi, si la statisticienne choisit la machine $k$, avec $100 \times p_k \%$ de chance, elle recevra un gain de $1$ et avec $100 \times (1 - p_k) \%$ de chance, $0$. Si la statisticienne connait les paramètres $p_1,p_2,\dots,p_K$ des machines à sous, pour maximiser ses gains, il suffit de jouer à tous les coups la machine avec le paramètre le plus grand. Mais elle ne connait pas ces paramètres. Il lui faut "explorer" les différentes machines tout en cherchant à maximiser ses gains. Nous allons l'aider à trouver des stratégies.

Simulation du problème :¶

Nous allons simuler le problème dans le cas où la statisticienne choisit à chaque coup la machine de façon complètement aléatoire.

In [1]:
import numpy as np
import random
import math
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.image as mpimg
from IPython import display
from IPython.display import HTML
random.seed(0)
# Choix du nombre machine à sous
K = 3

# Génération des paramètres des lois de Bernoulli
# p = np.random.uniform(0,1,K)

# Ou définition du problème de bandit à K bras que l'on souhaire résoudre
p = np.array([0.5,0.2,0.6])
print(p)

# Nombre de coups à jouer 
T = 100

# Matrice qui contiendra les gains de chaque machine au cours des coups joués
G = np.zeros([K, T])

# Vecteur qui contiendra les numéros des machines jouées 
# (par convention python la première est machine est numérotéé 0)
It = -1*np.ones(T)

# pour chaque coup
for t in range(T):
    # la statisticienne choisit la machine à sous it de manière aléatoire 
    it = math.floor(random.uniform(0,K))
    # it est stocké dans le vecteur It
    It[t] = it 
    # elle reçoit le gain gt : il est tiré selon la loi de Bernoulli de paramètre p[it]
    gt = random.binomialvariate(1,p[it])
    # pour chaque machine à sous k
    for k in range(K):
        # si k est la machiné jouée, alors on met à jour le gain associé à la machine k
        if k == it:
            G[k,t] = G[k,t-1] + gt
        # si k n'est pas la machiné jouée, le gain associé à la machine k reste le même
        else:
            G[k,t] = G[k,t-1] 

[0.5 0.2 0.6]

Illustration du problème avec une animation:¶

In [2]:
# Code pour générer l'animation

let = mpimg.imread('let.png')
pick = mpimg.imread('pick.png')
white = mpimg.imread('white.png')
end = mpimg.imread('end.png')

def JouerBandit(frame): 
    if frame == T:
        for k in range(K):
            im = white
            axs[k].imshow(im)
            axs[k].axis("off")
            axs[k].set_title("", y=-K/10)
        axs[math.floor(K/2)].imshow(end)
        axs[math.floor(K/2)].axis("off")
        axs[math.floor(K/2)].set_title("Gain = %i" %np.sum(G[:, T-1]), y=-K/10, 
                                       fontsize=20, color= "#01A9B1")
    else:
        it = It[frame]
        for k in range(K):
            if(k == it):
                im = pick
            else:
                im = let
            axs[k].imshow(im)
            axs[k].axis("off")
            axs[k].set_title("%i" %G[k, frame], y=-K/10)
   
f, axs = plt.subplots(1,K)    
anim_created = animation.FuncAnimation(f,JouerBandit, frames=T+1, interval=300)
video = anim_created.to_jshtml()
html = display.HTML(video)
display.display(html) 
plt.close()
No description has been provided for this image

Stratégie Explore then commit¶

La staticienne explore de façon déterministe pendant $K \times T_{\mathrm{explore}}$ coups : elle joue chaque machine $T_{\mathrm{explore}}$ fois (phase Explore). Elle approxime ensuite $p_k$ pour chaque machine à sous $k$, grâce à la moyenne empirique. Le reste du temps, elle va jouer la machine avec la moyenne empirique la plus grande, c'est-à-dire celle qui lui a rapporté le plus lors de la phase d'exploration (phase Commit).

In [3]:
# Explore then commit 

# Matrice de taille K x T qui contiendra les gains de chaque machine au cours des coups joués
G_etc = np.zeros([K, T])

# Vecteur de taille T qui contiendra les numéros des machines jouées 
It_etc  = -1*np.ones(T)

# Vecteur de taille K qui contiendra le nombre de fois que l'on a joué chaque machine 
N = np.zeros([K])

# proportion du temps où l'on explorera 
prop_explore = 0.2

# calcul de T_explore en fonction 
T_explore = math.floor(prop_explore*T/(K))

# Algorithme Explore then commit : 
for t in range(T):
    # explore
    if t<(K*T_explore):
        it = math.floor(t/T_explore)

    #commit 
    else:
        it = np.argmax(G_etc[:,K*T_explore-1])
    gt = random.binomialvariate(1,p[it])
    It_etc[t] = it 
    for k in range(K):
        if k == it:
            G_etc[k,t] = G_etc[k,t-1] + gt
            N[k] = N[k] + 1
        else:
            G_etc[k,t] = G_etc[k,t-1] 
            
def JouerBanditETC(frame): 
    if frame == T:
        for k in range(K):
            im = white
            axs[k].imshow(im)
            axs[k].axis("off")
            axs[k].set_title("", y=-K/10)
        axs[math.floor(K/2)].imshow(end)
        axs[math.floor(K/2)].axis("off")
        axs[math.floor(K/2)].set_title("Gain = %i" %np.sum(G_etc[:, T-1]), y=-K/10, 
                                       fontsize=20, color= "#01A9B1")
    else:
        it = It_etc[frame]
        for k in range(K):
            if(k == it):
                im = pick
            else:
                im = let
            axs[k].imshow(im)
            axs[k].axis("off")
            axs[k].set_title("%i" %G_etc[k, frame], y=-K/10)

f, axs = plt.subplots(1,K)    
anim_created = animation.FuncAnimation(f,JouerBanditETC, frames=T+1, interval=300)
video = anim_created.to_jshtml()
html = display.HTML(video)
display.display(html) 
plt.close()
No description has been provided for this image

Stratégie Upper confidence bound¶

L'algorithme Upper confidence bound (UCB, Auer et al., 2002) repose sur un principe d'optimisme devant l'incertidude: pour chaque machine à sous $k$, la statisticienne approxime le paramètre $p_k$ (toujours avec la moyenne empirique - notons $\widehat{p}_{k,t}$ l'estimation de $p_k$ au coup $t$). Mais elle va plus loin: à chaque coup $t$, elle construit aussi un intervalle de confiance $\alpha_{k,t}$ autour de $p_k$. Avec très grande probabilité, le vrai paramètre $p_k$ est dans l'intervalle $[\widehat{p}_{k,t} - \alpha_{k,t} ,\widehat{p}_{k,t} + \alpha_{k,t} ]$. C'est là qu'intervient l'optimisme: elle imagine que pour les $p_k$ appartiennent effectivement aux intervalles de confiance, mais surtout que toutes les machines autant gagnantes qu'il l'est très probablement possible; autrement dit, elle imagine que pour toutes les machines à sous, $p_k = \widehat{p}_{k,t} + \alpha_{k,t}$ et choisit donc $I_{t+1} \in \underset{k \in \{1,\dots,K\}}{\mathrm{argmax}}\{ \widehat{p}_{k,t} + \alpha_{k,t} \}$.

L'enjeu réside à présent dans le choix des $\alpha_{k,t}$. La théorie (essentiellement l'inégalité de concentration d'Hoeffding) donne l'expression explicite: $$I_{t+1} \in \underset{k \in \{1,\dots,K\}}{\mathrm{argmax}} \Bigg\{ \widehat{p}_{k,t} + \sqrt{\frac{2 \log t}{N_{k,t}}} \Bigg\} $$ où $N_{k,t}$ le nombre de fois que la statisticienne à jouer la machine $k$ après le coup $t$. Cela semble un peu compliqué mais retenons deux choses essentielles:

$\qquad \bullet$ $\alpha_{k,t}$ décroit quand $N_{k,t}$ augmente car plus la machine $k$ est jouée, plus nous sommes confiantes sur l'estimation de $\widehat{p}_{k,t} $

$\qquad \bullet$ $\alpha_{k,t}$ augmente tout doucement (via le logarithme) avec le nombre de coups $t$, cela permet de revenir jouer sur une machine à sous qui n'a pas été choisie depuis longtemps et te garantir un bon niveau d'exploration tout au long de la partie.

In [4]:
G_ucb = np.zeros([K, T])
It_ucb = -1*np.ones(T)
ucb = np.zeros([K, T])
N = np.zeros([K,T])

for t in range(T):
    if t<K:
        it = t
    else:
        it = np.argmax(ucb[:,t-1])
    gt = random.binomialvariate(1,p[it])
    It_ucb[t] = it 
    for k in range(K):
        if k == it:
            G_ucb[k,t] = G_ucb[k,t-1] + gt
            if t == 1:
                N[k,t] = 1*(k == it)
            else:
                N[k,t] = N[k,t-1] + 1
        else:
            G_ucb[k,t] = G_ucb[k,t-1] 
            N[k,t] = N[k,t-1] 
        ucb[k,t] = G_ucb[k,t]/max(N[k, t],1) + math.sqrt(2*math.log(t+1)/max(N[k,t],1))

def PlotUCB(t): 
    ylm = -0.3
    ylM = 3
    plt.cla()   
    plt.clf()
    plt.ylim((ylm,ylM))
    plt.xlim((-0.4,K-0.6))
    it = It_ucb[t]
    for k in range(K):
        if k == it:
            c = "#01A9B1"
        else:
            c = '#22509D'
        plt.axvline(x= k , 
                ymin= (G_ucb[k,t]/max(N[k, t],1) - math.sqrt(2*math.log(max(t+1,2))/max(N[k,t],1)))/(ylM-ylm)-ylm,
                ymax= (G_ucb[k,t]/max(N[k, t],1) + math.sqrt(2*math.log(max(t+1,2))/max(N[k,t],1)))/(ylM-ylm)-ylm,
                color=c, label='Machine 1',
                marker = 'x', markersize=7)
        plt.plot(k,(G_ucb[k,t]/max(N[k, t],1)),c,marker='o') 
        plt.plot(k,p[k],'k',marker='*', markersize=7) 
    plt.tick_params(
    axis='x',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom=False,      # ticks along the bottom edge are off
    top=False,         # ticks along the top edge are off
    labelbottom=False)


f, ax1 = plt.subplots(1,1)
anim_created = animation.FuncAnimation(f,PlotUCB, frames=T, interval=300)
video = anim_created.to_jshtml()
html = display.HTML(video)
display.display(html) 
plt.close()
No description has been provided for this image

Dans l'animation précédante, chaque segment correspond à l'intervalle de confiance $[\widehat{p}_{k,t} - \alpha_{k,t} ,\widehat{p}_{k,t} + \alpha_{k,t} ]$ pour chaque des machines à sous $k = 1,\dots,K$, les étoiles noires correspondent aux vraies valeurs $p_k$ et les ronds aux valeurs estimées $\widehat{p}_{k,t}$. La machine jouée (celle qui a l'intervalle de confiance le plus haut) s'allume en bleu clair comme lors des premières animations.

In [5]:
def JouerBanditUCB(frame): 
    if frame == T:
        for k in range(K):
            im = white
            axs[k].imshow(im)
            axs[k].axis("off")
            axs[k].set_title("", y=-K/10)
        axs[math.floor(K/2)].imshow(end)
        axs[math.floor(K/2)].axis("off")
        axs[math.floor(K/2)].set_title("Gain = %i" %np.sum(G_ucb[:, T-1]), y=-K/10, 
                                       fontsize=20, color= "#01A9B1")
    else:
        it = It_ucb[frame]
        for k in range(K):
            if(k == it):
                im = pick
            else:
                im = let
            axs[k].imshow(im)
            axs[k].axis("off")
            axs[k].set_title("%i" %G_ucb[k, frame], y=-K/10)

f, axs = plt.subplots(1,K)    
anim_created = animation.FuncAnimation(f,JouerBanditUCB, frames=T+1, interval=300)
video = anim_created.to_jshtml()
html = display.HTML(video)
display.display(html) 
plt.close()
No description has been provided for this image

Comparaison entre les trois stratégies (Aléatoire, Explore then commit, Upper confidence bound)¶

In [6]:
n = 200
T_explore = math.floor(prop_explore*T/(K))

res_rd = np.zeros([n, T])
res_etc = np.zeros([n, T])
res_ucb = np.zeros([n, T])

# pour chaque simulation i 
for i in range(n):
    
    # Simulation des gains pour chaque machine
    G = np.zeros([K, T])
    for k in range(K):
        for t in range(T):
            G[k,t] = random.binomialvariate(1,p[k])
            
    # Comparaison des gains de chaque stratégie
    G_rd = np.zeros([K, T])
    G_etc = np.zeros([K, T])
    G_ucb = np.zeros([K, T])
    ucb = np.zeros([K, T])
    N_ucb = np.zeros([K])
    for t in range(T):
        
        # Choix de la machine 
        # Aléatoire
        it_rd = math.floor(random.uniform(0,K))
        gt_rd = random.binomialvariate(1,p[it_rd])
        # Explore then commit 
        if t<(K*T_explore):
            it_etc = math.floor(t/T_explore)
        else:
            it_etc = np.argmax(G_etc[:,K*T_explore-1])
        # UCB
        if t<K:
            it_ucb = t
        else:
            it_ucb = np.argmax(ucb[:,t-1])
        
        # Mise à jour des gains de chaque stratégie    
        for k in range(K):
            if k == it_rd:
                G_rd[k,t] = G_rd[k,t-1] + G[k,t]
            else:
                G_rd[k,t] = G_rd[k,t-1] 
            if k == it_etc:
                G_etc[k,t] = G_etc[k,t-1] + G[k,t]
            else:
                G_etc[k,t] = G_etc[k,t-1] 
            if k == it_ucb:
                G_ucb[k,t] = G_ucb[k,t-1] + G[k,t]
                N_ucb[k] = N_ucb[k] + 1
            else:
                G_ucb[k,t] = G_ucb[k,t-1] 
            ucb[k,t] = G_ucb[k,t]/max(N_ucb[k],1) + math.sqrt(2*math.log(t+1)/max(N_ucb[k],1))
    res_rd[i, :] = np.sum(G_rd, axis=0)
    res_etc[i, :] = np.sum(G_etc, axis=0)
    res_ucb[i, :] = np.sum(G_ucb, axis=0)

x = np.linspace(0, T, T)
y_rd = np.mean(res_rd, axis=0)
y_etc = np.mean(res_etc, axis=0)
y_ucb = np.mean(res_ucb, axis=0)

sd_rd = np.std(res_rd, axis=0)
sd_etc = np.std(res_etc, axis=0)
sd_ucb = np.std(res_ucb, axis=0)

fig, ax = plt.subplots()
ax.plot(x, y_rd, label="Aléatoire",color= "#223F6A")
ax.fill_between(x, y_rd + sd_rd, y_rd - sd_rd, color= "#223F6A", alpha = 0.5)
ax.plot(x, y_etc, label="Explore then commit",color= "#5FCC9E")
ax.fill_between(x, y_etc + sd_etc, y_etc - sd_etc, color= "#5FCC9E", alpha = 0.5)
ax.plot(x, y_ucb, label="Upper Confidence Bound",color= "#326DC0")
ax.fill_between(x, y_ucb + sd_ucb, y_ucb - sd_ucb, color= "#326DC0", alpha = 0.5)
ax.set_xlabel ('Itérations')
ax.set_ylabel ('Gain cumulé')
ax.legend()
plt.show()
No description has been provided for this image