Implementing Bates distribution

Question

I've been trying to plot the Bates distribution curve, The Bates distribution is the distribution of the mean of n independent standard uniform variates (from 0 to 1).

(I worked on the interval [-1;1], I made a simple change of variable).

The curve destabilizes after such number of n, which prevents me from moving forward. In order to consider that the variable x is continuous, I sampled interval in 10**6 samples. Here are some examples for different n:

Bates distribution for some different n

But for n greater than 29, the curve diverges, and the greater n, the closer the deformation caused by the divergence is to the (mean) center of the curve:

Divergence of the distribution

The Bates distribution of probability is defined as follows:

Bates distribution

My code:

samples=10**6

def combinaison(n,k):   # combination of K out of N
  cnk=fac(n)/(fac(k)*fac(abs(n-k))) # fac is factoriel 
  return cnk


def dens_probas(a,b,n):
  x=np.linspace(a, b, num=samples)
  y=(x-a)/(b-a)
  F=list()
  for i in range(0,len(y)):
    g=0
    for k in range(0,int(n*y[i]+1)):
      g=g+pow(-1,k)*combinaison(n,k)*pow(y[i]-k/n,n-1)
    d=(n**n/fac(n-1))*g
    F.append(d)         
  return F

Any idea to correct the divergence for larger n?

user_12046409 · Answer

The main problem is that the formula with alternating sums is extremely prone to numerical accuracy issues.

One trick to avoid the problems on the right side, is to assume the distribution is symmetric and only calculate half of it.

A straightforward accuracy optimization is to replace the factorials in the formula for combinaison by a call to scipy.special.comb. This avoids that very large numbers need to be divided.

A smaller accuracy optimization is to calculate g for even and odd numbers together. But at first sight the formula can not be reduced much, so replacing:

        for k in range(0, int(floor(n * y[i] + 1))):
            g += pow(-1, k) * combinaison(n, k) * pow(y[i] - k / n, n - 1)

By:

        last_k = int(floor(n * y[i]))
        for k in range(0, last_k + 1, 2): # note that k increments in steps of 2
            if k == last_k:
                g += combinaison(n, k) * (pow(y[i] - k / n, n - 1))
            else:
                g += combinaison(n, k) * (pow(y[i] - k / n, n - 1) - pow(y[i] - (k + 1)/ n, n - 1) * (n - k) / (k + 1))

Some other remarks:

The variable samples is only used to tell the division in the xaxis. A much smaller number will suffice. (In the code below I renamed the variable to xaxis_steps).
Using append for F will be extremely slow. It is better to create a numpy array of the correct size and then fill it in. (This also makes the copying of the halves easier.)

from matplotlib import pyplot as plt
import numpy as np
from scipy.special import comb
from math import factorial as fac
from math import floor

xaxis_steps = 500

def combinaison(n, k):  # combination of K out of N
    return comb(n, k)

def dens_probas(a, b, n):
    x = np.linspace(a, b, num=xaxis_steps)
    y = (x - a) / (b - a)
    F = np.zeros_like(y)
    for i in range(0, (len(y)+1) // 2):
        g = 0
        for k in range(0, int(floor(n * y[i] + 1))):
            g += pow(-1, k) * combinaison(n, k) * pow(y[i] - k / n, n - 1)
        F[i] = (n ** n / fac(n - 1)) * g
        F[-i-1] = F[i]  # symmetric graph
    plt.plot(x, F, label=f'n={n}')
    return F

for n in (5, 30, 50, 80, 90):
    dens_probas(-1, 1, n)
plt.legend()
plt.show()

All these optimizations together move the accuracy problem from n=30 to around n=80:

A completely different approach would be to generate a lot of uniform samples and take the means. From those samples a kde plot can be generated. The smoothness of such curves depends on the number of samples. A kde can be plotted directly via seaborn's kdeplot. You can also separately calculate the kde function, then apply it to a given x range and plot it via standard matplotlib.

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kde

num_samples = 10 ** 5

def dens_probas(a, b, n):
    samples = np.random.uniform(a, b, size=(num_samples, n)).mean(axis=1)
    samples = np.hstack([samples, a + b - samples])  # force symmetry; this is not strictly necessary
    return gaussian_kde(samples)

for n in (5, 30, 50, 80, 90, 200):
    kde = dens_probas(-1, 1, n)
    xs = np.linspace(-1, 1, 1000)
    F = kde(xs)
    plt.plot(xs, F, label=f'n={n}')
plt.legend()
plt.show()

Implementing Bates distribution

1 Answers