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.
 A: 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 
\begin{equation}
\mu_x \mid \textrm{data} \sim \mathrm{Beta}(\alpha+y_x,\beta+n_x-y_x).
\end{equation}
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).
