# 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()

-
Note that $S$ cannot be exactly uniform, because it is a discrete variable. You are therefore asking about asymptotic behavior. Moreover, for small $N$ (less than $100/\min(p,1-p)$, approximately, with $p=p_A=p_B$) the distribution isn't even remotely close to uniform. –  whuber May 12 at 15:04
@whuber S is not discrete, it is a probability that can fall between 0 and 1. Also, even for low N, I am observing uniform behaviour. –  Cam.Davidson.Pilon May 12 at 15:22
I must be misunderstanding your setup, then. As far as I can tell, for any given values of $N,n_A,n_B,$ the value of $S$ is a number. Therefore, accepting that $N, p_A,$ and $p_B$ are fixed for the moment (as they are in your code), $S$ is a function of $(n_A,n_B)$. But the latter, being realizations of two Binomial distributions, can attain only a discrete set of values. When I reproduce your code in R, I get decidedly non-uniform histograms for small $N$. –  whuber May 12 at 15:28
Although indeed your $S$ has values between $0$ and $1$, do not confuse that with non-discrete: it can have at most $N^2$ distinct values (and actually has fewer than that). This might not be perfectly clear to you because your simulation generates estimates of $S$ rather than its correct values and the estimates essentially do have a continuous distribution. –  whuber May 12 at 15:35
@whuber yes, you are correct, excellent observation. I'm still stuck on why it looks uniform then. –  Cam.Davidson.Pilon May 12 at 15:37

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()

-
So there is a difference between mine and your code. I sample A and B in each loop, you sample it once and calculate S 5000 times. –  Cam.Davidson.Pilon May 12 at 16:08
The discrepancy lies in your calls to rbinom, which returns a vector. The subsequent call to rbeta inside replicate is vectorized, so the inner (internal) loop is using a different $A$ and $B$ for each of the 15000 random variables generated (wrapping around for the final 5000 since your NSIM = 10000). See ?rbeta for more. This differs from @Cam's code with has a single fixed $A$ and $B$ used in all 15000 random-variate calls for each of the 5000 sampling (replicate) loops. –  cardinal May 12 at 17:57
here's the output for those curious: imgur.com/ryvWbJO –  Cam.Davidson.Pilon May 12 at 18:37
The only things I'm aware of that are potentially pertinent at a conceptual level are that a) the expected distribution of results are symmetric, b) a bin size of 1 is always uniform, c) a bin size of 2 for a symmetric distribution will also always appear uniform, d) the number of possible sampling distributions that can be drawn from increases with N, e) values of S can't stack up on 0 or 1 alone because beta is undefined when there are 0 successes in either group, and f) the samples are restricted between 0 and 1. –  rpierce May 12 at 19:21
As a matter of observation alone we can see that the distances between the centroids of the sampling distributions get reduced as the centroids of sampling distributions moves away from .5 (probably related to point f above). This effect tends to counteract the tendency for the high frequencies of observations for the more common nearly equal successes in group A and group B case. However, giving a mathematical solution as to why that is or why it should yield normal distributions for certain bin sizes isn't anywhere near my territory. –  rpierce May 12 at 19:29

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.

-
I'mn with you up until $P_A - P_B \approx Normal$... then you go off another direction that I didn't quite follow. Is $\Phi$ defined twice, once as the standard normal CDF and then as the probability integral transform? I'm hoping you can expand your description around these steps and relate them to the initial code/problem. Maybe loop back around and restate which specific parameters produce the uniform result. –  rpierce May 12 at 20:40
@rpierce (1) The difference $P_A-P_B$ is approximately normal because $P_A$ and $P_B$ are independent and each is approximately normal. The mean is the difference of the means and the variance is the sum of the variances. (2) The probability integral transform is the CDF: it is the case for any random variable $X$ with continuous distribution $F$, that $F(X)$ is uniform. –  whuber May 12 at 21:12
Oh I got 1, it was the stuff after it where I got lost. This will be mindbogglingly dumb, but why is $Pr(P_A>P_B)$ the same as the CDF? –  rpierce May 12 at 21:25
@rpierce That follows rather directly from the definition, but there's a slight twist in which the symmetry of the Normal distribution is invoked. We're dealing with a Normal variate $X = P_A-P_B$ assumed to have an expectation of $\mu=\alpha-\beta$ and variance $\sigma^2 = 2p(1-p)/N$. Standardizing $X$, it is natural to rewrite the probability as $$\Pr(X\gt 0) = \Pr((X-\mu)/\sigma \gt (0-\mu)/\sigma) = 1-\Phi(-\mu/\sigma) = \Phi(\mu/\sigma).$$ –  whuber May 12 at 21:30
@whuber this is pretty amazing. You are a wonderful teacher. I appreciate both yours and rpierce's answer, I'll still give him credit as it did solve our problem, and you have shown why the behaviour occurs. Ty! –  Cam.Davidson.Pilon May 12 at 22:21