why curve_fit does not converge for a beta function fit?

Question

I have a small problem with my Python code when I try to fit a beta function to a few points. The problem is that either the solution does not converge (and the result coefficients are nans), or it does nothing (at the results stay the same as my initial guess), or it apparently fits, but then the fit is not similar to the data points at all. I have been reading similar posts about beta function and about curve_fit, because both are discussed questions in the stackoverflow literature, but I have not been able to find a solution to the specific problem I have, so I was wondering if you could give me some ideas.

I have a set of points:

x = np.array([0.1, 0.3, 0.5, 0.7, 0.9, 1.1])
y = np.array([0.45112234, 0.56934313, 0.3996803 , 0.28982859, 0.19682153,
   0.]

and then I try to fit them with the gamma function using curve_fit as follows:

from scipy.optimize import curve_fit
from scipy.special import gamma as gamma
def betafunc(x,a,b,cst):
    return cst*gamma(a+b) * (x**(a-1)) * ((1-x)**(b-1))  / ( gamma(a)*gamma(b) )
popt2,pcov2 = curve_fit(betafunc,x,y,p0=(0.5,1.5,0.5))

And there is where my problem arises, because depending on my initial guess, I get either popt2=[nan,nan,nan], or popt2=p0, or a few times values that, when plotted, do not mimic the data at all.

I also know the beta function is for 0 < x < 1, so I have tried re-scaling the points, or just removing the last point of my data, but that does not well either. Adding errors to curve fit or, as mentioned, changing the initial parameters, also does not help. Also I thought it could be just the fact that I have 3 free parameters and 4 or 5 points, but, as shown in the figure... ,

...I have already fitted another profile (that also uses three free parameters), and there is not problem with it, so I do not understand why this other beta profile does not work. Any guidance is warmly appreciated!

Answer 1

Your implementation is for the probability density function of the beta distribution (rather than the beta function). It is defined over the interval 0 <= x <= 1. Therefore, your independent data (x coordinate) has to be entirely within this interval. How you ensure that, depends on the larger context of what you are trying to do. Possible approaches would include pruning or mapping onto some portion of the interval.

Following your example, the following (pruning) runs without error:

import numpy as np
from scipy.optimize import curve_fit
from scipy.special import gamma as gamma
def betafunc(x,a,b,cst):
    return cst*gamma(a+b) * (x**(a-1)) * ((1-x)**(b-1))  / ( gamma(a)*gamma(b) )

x = np.array( [0.1, 0.3, 0.5, 0.7, 0.9, 1.1])
y = np.array( [0.45112234, 0.56934313, 0.3996803 , 0.28982859, 0.19682153, 0.] )


popt2,pcov2 = curve_fit(betafunc,x[:-1],y[:-1],p0=(0.5,1.5,0.5))

print popt2
print pcov2

and produces the result:

[ 1.22624727  1.74192827  0.37084996]
[[ 0.03758865  0.04888083 -0.00132468]
 [ 0.04888083  0.09142608 -0.00309165]
 [-0.00132468 -0.00309165  0.00094766]]

And this also (rescaling x), runs without error:

import numpy as np
from scipy.optimize import curve_fit
from scipy.special import gamma as gamma
def betafunc(x,a,b,cst,scale):
    x = x / scale
    return cst*gamma(a+b) * (x**(a-1)) * ((1-x)**(b-1))  / ( gamma(a)*gamma(b) )

x = np.array( [0.1, 0.3, 0.5, 0.7, 0.9, 1.1])
y = np.array( [0.45112234, 0.56934313, 0.3996803 , 0.28982859, 0.19682153, 0.] )

popt2,pcov2 = curve_fit(betafunc,x,y,p0=(0.5,1.5,0.5,1.1))

print popt2
print pcov2

and produces the result:

[ 1.37100253  2.36832069  0.32337175  1.16052822]
[[ 0.04972377  0.15943756 -0.00792804  0.02550767]
 [ 0.15943756  0.71001918 -0.04180131  0.14426687]
 [-0.00792804 -0.04180131  0.00312037 -0.00983075]
 [ 0.02550767  0.14426687 -0.00983075  0.0373759 ]]

Note that in the second example, the x range scaling is one of the fit variables also. But it could just as well be constant determined by your problem. It all depends on your context.

Again, which approach is the correct one for you to use, depends on the details of where the data comes from. The approach that you chose should make physical sense in the context of the data that you are trying to fit and what you hope to accomplish.