Why is this distribution uniform? We are investigating Bayesian statistical testing, and come across an odd (to me atleast) phenomenon. 
Consider the following case: we are interested in measuring which population, A or B, has a higher conversion rate. For a sanity check, we set $p_A = p_B$, that is, the probability of conversion is equal in both groups. We generate artificial data using a binomial model, e.g. $$n_A \sim \text{Binomial}(N, p_A)$$
We then try to estimate the $p_A, p_B$ using a Bayesian beta-binomial model so we get posteriors for each conversion rate, e.g. $$P_A \sim \text{Beta}(1 + n_A, N - n_A +1 )$$
Our test statistic is computed by calculating $S = P(P_A > P_B\; |\; N, n_A, n_B)$ via monte carlo. 
What surprised me was that if $p_A = p_B$, then $S \sim \text{Uniform(0,1)}$. My thoughts were that it would be centered around 0.5, and even converge to 0.5 as the sample size, $N$, grows.  
My question is, why is  $S \sim \text{Uniform(0,1)}$ when $p_A = p_B$?

Here's some Python code to demonstrate:
%pylab
from scipy.stats import beta
import numpy as np
import pylab as P

a = b = 0.5
N = 10000
samples = [] #collects the values of S
for i in range(5000):
    assert a==b
    A = np.random.binomial(N, a); B = np.random.binomial(N, b)
    S = (beta.rvs(A+1, N-A+1, size=15000) > beta.rvs(B+1, N-B+1, size=15000)).mean() 
    samples.append(S)

P.hist(samples)
P.show()

 A: To get some intuition for what is going on, let us feel free to make $N$ very large and in so doing ignore $O(1/N)$ behavior and exploit asymptotic theorems that state both Beta and Binomial distributions become approximately Normal. (With some trouble, all this can be made rigorous.)  When we do this, the result emerges from a specific relationship among the various parameters.

Because we plan to use Normal approximations we will pay attention to the expectations and variances of the variables:


*

*As Binomial$(N, p)$ variates, $n_A$ and $n_B$ have expectations of $pN$ and variances of $p(1-p)N$.  Consequently $\alpha=n_A/N$ and $\beta=n_B/N$ have expectations of $p$ and variance $p(1-p)/N$.

*As a Beta$(n_A+1, N+1-n_A)$ variate, $P_A$ has an expectation of $(n_A+1)/(N+2)$ and a variance of $(n_A+1)(N+1-n_A) / [(N+2)^2(N+3)]$.  Approximating, we find that $P_A$ has an expectation of 
$$\mathbb{E}(P_A) = \alpha+O(1/N)$$
and a variance of 
$$\text{Var}(P_A) = \alpha(1-\alpha)/N + O(1/N^2),$$
with similar results for $P_B$.
Let us therefore approximate the distributions of $P_A$ and $P_B$ with Normal$(\alpha, \alpha(1-\alpha)/N)$ and Normal$(\beta,\beta(1-\beta)/N)$ distributions (where the second parameter designates the variance).  The distribution of $P_A-P_B$ consequently is approximately Normal; to wit,
$$P_A-P_B \approx \text{Normal}\left(\alpha-\beta, \frac{\alpha(1-\alpha) + \beta(1-\beta)}{N}\right).$$
For very large $N$, the expression $\alpha(1-\alpha) + \beta(1-\beta)$ will not vary appreciably from $p(1-p)+p(1-p)=2p(1-p)$ except with very low probability (another neglected $O(1/N)$ term).  Accordingly, letting $\Phi$ be the standard normal CDF,
$$\Pr(P_A\gt P_B) =\Pr(P_A-P_B\gt 0) \approx \Phi\left(\frac{\alpha-\beta}{\sqrt{2p(1-p)/N}}\right).$$
But since $\alpha-\beta$ has zero mean and variance $2p(1-p)/N,$ $Z=\frac{\alpha-\beta}{\sqrt{2p(1-p)/N}}$ is a standard Normal variate (at least approximately).  $\Phi$ is its probability integral transform; $\Phi(Z)$ is uniform.
A: TL;DR:  Mixtures of normal distributions may look uniform when bin sizes are large.
This answer borrows from @whuber's sample code (which I thought first was an error, but in retrospect was probably a hint).
The underlying proportions in the population are equal: a = b = 0.5.
Each group, A and B has 10000 members: N = 10000.
We are going to conduct 5000 replicates of a simulation: for i in range(5000):.  
Actually, what we are doing is a $\rm simulation_\rm{prime}$ of a $\rm simulation_\rm{underlying}$.  In each of the 5000 iterations $\rm simulation_\rm{prime}$ we will do $\rm simulation_\rm{underlying}$.
In each iteration of $\rm simulation_\rm{prime}$ we will simulate a random number of A and B that are 'successes' (AKA converted) given the equal underlying proportions defined earlier: A = np.random.binomial(N, a); B = np.random.binomial(N, b).  Nominally this will yield A = 5000 and B = 5000, but A and B vary from sim run to sim run and are distributed across the 5000 simulation runs independently and (approximately) normally (we'll be coming back to that).
Let's now step through  $\rm simulation_\rm {underlying}$ for a single iteration of $\rm simulation_\rm{prime}$ in which A and B have taken on an equal number of successes (as will be the average the case).  In each iteration of $\rm simulation_\rm{underlying}$ we will, given A and B, create random variates of the beta distribution for each group.  Then we will compare them and find out if ${\rm Beta}_A > {\rm Beta}_B$, yielding a TRUE or FALSE (1 or 0).  At the end of a run of $\rm simulation_\rm {underlying}$, we have completed 15000 iterations and have 15000 TRUE/FALSE values.  The average of these will yield a single value from the (approximately normal) sampling distribution of the proportion of ${\rm Beta}_A > {\rm Beta}_B$.
Except now $\rm simulation_\rm{prime}$ is going to select 5000 A and B values.  A and B will rarely be exactly equal, but the typical differences in the number of A and B successes are dwarfed by the total sample size of A and B.  Typical As and Bs will yield more pulls from their sampling distribution of proportions of ${\rm Beta}_A > {\rm Beta}_B$, but those on the edges of the A/B distribution will also get pulled.
So, what in essence we pull over many sim runs is a combination of sampling distributions of ${\rm Beta}_A > {\rm Beta}_B$ for combinations of A and B (with more pulls from the sampling distributions made from the common values of A and B than the uncommon values of A and B).  This results in mixtures of normal-ish distributions.  When you combine them over a small bin size (as is the default for the histogram function you used and was specified directly in your original code), you end up with something that looks like a uniform distribution.
Consider:
a = b = 0.5
N = 10
samples = [] #collects the values of S
for i in range(5000):
    assert a==b
    A = np.random.binomial(N, a); B = np.random.binomial(N, b)
    S = (beta.rvs(A+1, N-A+1, size=15000) > beta.rvs(B+1, N-B+1, size=15000)).mean() 
    samples.append(S)

P.hist(samples,1000)
P.show()

