# Fitting Beta Distributions to Data

I am trying to reproduce some beta distribution parameters found in this published paper. I have two data sets, y1 and y2, that are generated in the following way (NumPy code):

size = 5000
x = np.arange(size)

# Dataset #1
I = np.random.randint(0, size, size=size)
k = I.shape[0]
nnmark = np.zeros(k)
y1 = np.zeros(k)

for i in range(k):
j = I[i]
nnmark[min(i, j)] = nnmark[min(i, j)] + 1
nnmark[max(i, j)] = nnmark[max(i, j)] - 1

y1 = np.cumsum(nnmark)

# Dataset #2
I = np.empty(size, dtype=np.int)
for i in range(I.shape[0]):
I[i] = np.random.randint(i, size)

for i in range(k):
j = I[i]
nnmark[min(i, j)] = nnmark[min(i, j)] + 1
nnmark[max(i, j)] = nnmark[max(i, j)] - 1

y2 = np.cumsum(nnmark)


y1 is shown in black and y2 is shown in red above. According to the paper, both curves can be approximated by a beta distribution. In the original paper, they claim that since y1 has a width that is equal to the size and a height that is equal to 0.5*size then this is a special case of the beta distribution with beta(2, 2, a, c). Finally, their code shows that this any value from this distribution can be computed from:

y = 2*(x)*(size-z)/size


It isn't clear how this is derived. I am only familiar with basic beta distributions with alpha and beta parameters bounded between [0,1] and some references explain that a and c often refer to the location and scale of the distribution. I also noticed that the height of y1 and y2 are quite large and, unfortunately, I don't know how to create, fit, or rescale a beta distribution so as to obtain the correct parameters that would fit either of these curves.

My final goal would be to change the size (currently it is set to 5,000 for demonstration purposes but it can be smaller or larger) and be able to represent that dataset with an approximate beta distribution (albeit, it will likely need to be stretched horizontally and vertically).

I have been looking at the SciPy beta distribution function but the documentation is vague. I've gotten as far as:

a1, b1, c1, d1 = beta.fit(y1, loc=0, scale=size)
a2, b2, c2, d2 = beta.fit(y2, loc=0, scale=size)


But neither of the PDFs look like the original data when plotted next to it.

• The way you are using scale isn't correct. Try setting it to 1. c1 and d1 should be 0 or close to 0. Then try plotting x = np.linspace(0, 1, 100); plt.plot(x, beta.pdf(x, a1, b1)) Jun 23, 2019 at 13:37
• According to the (vague) documentation, it looks like fit should use fscale and not scale. I believe that loc and scale are actually ignored. Additionally, in order to specify floc and fscale together one needs to take care and have floc < x < floc+fscale (so one might need to add/subtract a very small constant 1e-10 if x contains values that land right on the boundary limits). The x_range of my data is between [0, 5000] and the y_range is between [0, 2500] for y1 and it isn't clear if I need to normalize it first before fitting.
– slaw
Jun 23, 2019 at 13:53
• Ah, you are right, it is fscale. Is this a demonstration of what you are trying to do? gist.github.com/CamDavidsonPilon/… Jun 23, 2019 at 19:17
• I think I realize what the problem is. y1 and y2 aren't "data". Instead, they are histograms that approximates some beta pdf. That is, I'm not drawing from a distribution and then fitting a beta distribution from the drawn data. Instead, what I am looking for is a way to fit a beta distribution to each histogram.
– slaw
Jun 24, 2019 at 1:04
• I am wondering if something like this where I can draw data samples from my histogram might be appropriate: stackoverflow.com/a/54535378/2955541
– slaw
Jun 24, 2019 at 1:28

To generate an accurate curve fit, you can sample the histogram via the scipy.stats.rv_histogram function and then average over the fitted parameters from multiple iterations:

n_iter = 1000
params = np.empty((n_iter, 2))
for i in range(n_iter):
n_samples = 1000
hist_dist = scipy.stats.rv_histogram((y2, np.append(np.arange(k), k)))
data = hist_dist.rvs(size=n_samples)
a, b, c, d = scipy.stats.beta.fit(data, floc=0, fscale=k)

params[i, 0] = a
params[i, 1] = b

a_mean = np.round(np.mean(params[:, 0]), 2)
b_mean = np.round(np.mean(params[:, 1]), 2)
print(a_mean, b_mean)

y2_normalized_fit = scipy.stats.beta.pdf(np.arange(k), a_mean, b_mean, loc=0, scale=k)
scaling_factor, _, _, _ = np.linalg.lstsq(y2_normalized_fit.reshape(-1,1), y2, rcond=None)


Then, you multiply the scaling_factor by the fitted curve to get the right curve.

y2_normalized_fit * scaling_factor