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)

Plotted Data

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.

  • $\begingroup$ 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)) $\endgroup$ – Cam.Davidson.Pilon Jun 23 '19 at 13:37
  • $\begingroup$ 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. $\endgroup$ – slaw Jun 23 '19 at 13:53
  • $\begingroup$ Ah, you are right, it is fscale. Is this a demonstration of what you are trying to do? gist.github.com/CamDavidsonPilon/… $\endgroup$ – Cam.Davidson.Pilon Jun 23 '19 at 19:17
  • $\begingroup$ 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. $\endgroup$ – slaw Jun 24 '19 at 1:04
  • $\begingroup$ I am wondering if something like this where I can draw data samples from my histogram might be appropriate: stackoverflow.com/a/54535378/2955541 $\endgroup$ – slaw Jun 24 '19 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

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.