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!

Could you post a complete code example? I guess `curve_fit` is from `scipy.optimize`. But what is `v2`? — Dietmar, Jun 04 '18 at 12:42
@Dietmar, I'm sorry, v2 should be changed, I will edit it now. Thanks — Enrique, Jun 04 '18 at 12:44
@Enrique Using your code and removing the last point (x>1) the fit converges without problems. — Christian K., Jun 04 '18 at 12:57
What is `-0.1 ** -0.5` supposed to mean without straying into the realm of imaginary numbers? You even get an error message: `RuntimeWarning: invalid value encountered in power` — Mr. T, Jun 04 '18 at 12:58
oh, yikes, you are right @ChristianK., and I have no clue how is it possible when I do the same in my proper code but it fails. Now I think most likely it is a problem in some other part of my code, maybe double definitions or something like that, but for sure the problem is somewhere else, because as you realized, it literally I only use the few lines in my question, it actually works. I will find out what the problem is but of course now this question is useless. In any case thanks for checking it! — Enrique, Jun 04 '18 at 13:09

DrM · Accepted Answer · 2018-06-04T13:50:43.887

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.

why curve_fit does not converge for a beta function fit?

1 Answers1