# Understanding default Bayesian posterior probability shape parameters in multi-armed bandit

Multi-Armed Bandit indicates a best "arm" based on a simple test of proportions and Bayesian posterior probability shape parameters, alpha and beta (some implementations may take a different approach, but I want to understand this one).

A simple example could be 3 slot machines I've played thousands of times, with the following results:

Machine A:   441 successes;   73,737 arm pulls
Machine B:   285 successes;   88,398 arm pulls
Machine C: 1,068 successes;  205,234 arm pulls


Machine A clearly has the highest proportion (0.005981) of successes. Given that I have no other information about the machines, I would expect a Multi-Armed Bandit to recommend that I play Machine A.

Surprisingly (to me) the Bandit strongly recommends Machine C (which has an estimated success-to-trial proportion of 0.005204 based on my sample of plays).

The Bayesian posterior probabilities that each arm is the best binomial bandit are:

Machine A: 1.799 E-15
Machine B: 2.378 E-15
Machine C: 2.675133 E-12


This surprises me. These results are computed with shape parameters alpha and beta both = 1. Based on these notes, I would've thought that with the shape parameters set to unity that the Bandit results would be the same as the proportions (successes divided by trials)?

For those unfamiliar with Bayesian posterior probabilities or Multi-Arm Bandits, I have provided links to both topics.

• How did you obtain the Bayesian posterior probabilities that each arm isthe best (right before the last paragraph)? How come they don't sum to 1? Apr 15, 2016 at 16:03
• @JuhoKokkala I will append the code to the question. To my knowledge, they aren't supposed to sum to 1 (?). Apr 15, 2016 at 16:42
• I think you should explain the algorithm and/or give a reference. The R code is not very informative as it just calls some package. And if the question is "how to interpret the output of this software package", it would be off-topic here. About the probabilities summing to 1: how is it possible that "A is not the best", "B is not the best" and "C is not the best" are simultaneously true? Or are these the probabilities of something else than "x is the best of these three"? Apr 15, 2016 at 16:54
• @JuhoKokkala This had nothing to do with R until you requested to see the calculation. This question is really intended for those who are already knowledgeable about Bayesian statistics (and to a lesser extent Multi-Arm Bandit), but I will include additional references to Bayesian statistics and MAB. Apr 15, 2016 at 17:02
• Sure, the means of the posteriors over win probabilities are approximately 0.0059, 0.0032, 0.0052. However we know more than that. The posterior standard deviations of those probabilities are approximately 0.00028, 0.00019, 0.00015, so we're much more certain about C than A. Hence, taking into account those different levels of uncertainty leads to the result you find surprising. This is sec 2.3 of the Scott paper you linked to, I believe. Apr 15, 2016 at 19:39

Based on the edit history, I assume the probabilities were computed using the R package bandit available on CRAN, as no other source for the numbers has been given. The results reported in the question are most likely due to numeric issues in this R package. Based on a quick look at the source code of the R package, it is performing some numeric integration to get these probabilities, I conjecture that this integration does not perform well when the posterior distributions of the win probabilities are too narrow.

In this answer, I explain how to get the posterior distribution of the win probabilities, how to simulate the probabilities (using Python), and finally compare to R bandit results (with smaller data, the approaches agree, which is to show that my results are likely not due to my misunderstandings).

### Model

This is the binomial bandit model, see Scott, S. L. (2010). A modern Bayesian look at the multi‐armed bandit. Applied Stochastic Models in Business and Industry, 26(6), 639-658.

Each arm has a win probability $\mu_x$ which has a beta prior $\mu_x \sim \mathrm{Beta}(\alpha,\beta)$. Conditional on $\mu_x$, the arm pulls are independent Bernoulli draws with probability $\mu_x$. The $\mu_x$-parameters of different arms are assumed to be independent.

This is the well-known beta-binomial conjugate case, so the posterior distribution of each win probability is $$\mu_x \mid \textrm{data} \sim \mathrm{Beta}(\alpha+y_x,\beta+n_x-y_x).$$ Furthermore, since the parameters are a priori independent and the observation models are independent, they are also independent a posteriori.

### Getting the probabilities that arm x has the highest win probability

So, we are dealing with independent Beta distributed parameters. I suspect an analytic expression for the integral does not exist, but a solution is to simply simulate by drawing a beta random variable from each arm's posterior distribution, check which is best, repeat this $N$ times and count the frequencies. Python code:

def best_probs(data,alpha,beta,N):
winprob_samples = np.zeros((data.shape[0],N))

for i in range(data.shape[0]):
winprob_samples[i,:] = np.random.beta(alpha+data[i,0],beta+data[i,1]-data[i,0],size=N)

best_samples = np.argmax(winprob_samples,axis=0)
best_probs = np.bincount(best_samples)/N


With the parameters in the question:

import numpy as np
np.random.seed(1)

data = np.array([[441,73737],[285,88398],[1068,205234]])
alpha = 1
beta = 1

bp = best_probs(data,alpha,beta,1000000)
print(bp)


I get about $0.993$ for arm A and $0.007$ for arm $B$. So, A is indeed very likely to be the best, as would be intuitive from the data.

### Check against the R package bandit

To check that I have understood correctly what the R package is aiming to do and that my Python implementation works as intended, let us check with smaller data (so that the R package does not run into numeric problems). With the data:

• Arm 1: 125 pulls, 50 wins
• Arm 2: 300 pulls, 100 wins
• Arm 3: 400 pulls, 150 wins

I get win probabilities $0.68357, 0.03158, 0.28485$ using the R package. Meanwhile, my Python:

np.random.seed(1)
data = np.array([[50,125],[100,300],[150,400]])
alpha = 1
beta = 1
bp = best_probs(data,alpha,beta,1000000)
print(bp)
[ 0.684182  0.031538  0.28428 ]


So, based on this test, the implementations agree (when the posteriors are not too narrow for the numeric integration the R package performs).