scipy.optimize.curvefit: Asymmetric error in fit

Question

I try to fit a function to my data using scipy.optimize.curvefit.

Q=optimization.curve_fit(func,X,Y, x0,ERR)

and it works well.

However, now I am trying to use an asymmetric error and I have no idea how to do that - or even if it is possible.

By asymmetric error I mean that the error is not for example: 3+-0.5 but 3 +0.6 -0.2. So that ERR is an array with two columns.

It would be great if somebody had an idea how to do that - or could me point to a different Python routine which might be able to do it.

That a snippet of the code I am using - but I am not sure it makes it clearer:

A=numpy.genfromtxt('WF.dat')
cc=A[:,4]
def func(A,a1,b1,c1):
    N=numpy.zeros(len(x))
    for i in range(len(x)):
        N[i]=1.0*erf(a1*(A[i,1]-c1*A[i,0]**b1))

return N


x0   = numpy.array([2.5  , -0.07 ,-5.0])
Q=optimization.curve_fit(func,A,cc, x0, Error)

And Error=[ErP,ErM] (2 columns)

Can you post some (preferably runnable) code so that we can see the problem in more detail? — YXD, Oct 01 '13 at 12:53
not really. That would defeat the purpose of having such precise errors. — user2622479, Oct 01 '13 at 13:54
This problem might be more fun for people to tackle if you did as Mr E suggested and include a simple, runnable example that highlights the problem. Instead of loading data from a file, generate an asymmetric gaussian and add some random noise (np.random.random). Then solutions can try to fit this and get an asymmetric error for the mean. — DanHickstein, Oct 01 '13 at 15:38

score 1 · Answer 1 · edited Jul 30 '17 at 10:51

In the current version, I am afraid it is not doable. curve_fit is a wrap around the popular Fortran library minipack. Check the source code of \scipy_install_path\optimize\minipack.py, you will see: (line 498-509):

if sigma is None:
    func = _general_function
else:
    func = _weighted_general_function
    args += (1.0/asarray(sigma),)

Basically what it means is that of sigma is not provided, then the unweighted Levenberg-Marquardt method in minipack will be called. If sigma is provided, then the weighted LM will be called. That implies, if sigma is to be provided, it must be provided as a array of the same length of X or Y.

That means if you want to have asymmetric error residue on Y, you have to come up with some modification to your target function, as @Jaime suggested.

score 1 · Answer 2 · answered Oct 02 '13 at 02:31

Least squares algorithm like curve_fit or scipy.optimize.leastsq will not be able to do this because the loss function is quadratic, and so symmetric for positive and negative error.

I haven't seen any models for this and maybe PAIDA can handle it, as DanHickstein mentioned.

Otherwise, you could use the nonlinear optimizers like optimize.fmin and construct your own asymmetric loss function.

def loss_function(params, ...):
    error = (y - func(x, params))
    error_neg = (error < 0)
    error_squared = error**2 / (error_neg * sigma_low + (1 - error_neg) * sigma_upp))
    return error_squared.sum()

and minimize this with fmin or fmin_bfgs.

(I never tried this.)

score 0 · Answer 3 · answered Oct 01 '13 at 15:34

0

I'm not 100% sure, but it looks like the PAIDA package might do fits with asymmetric errors:

http://paida.sourceforge.net/documentation/fitter/index.html

answered Oct 01 '13 at 15:34

DanHickstein

6,588
13
54
90

score 0 · Answer 4 · answered Jan 30 '19 at 10:23

A solution, which I've used frequently, is to draw realisations (say 100-1000) from a split-normal distribution, and run your fitting algorithm on each of these realisations with the error set to 0.0. You'll then have 100-1000 best-fitting parameters, from which you can simply take the median, along with any error estimate you want to use.

scipy.optimize.curvefit: Asymmetric error in fit

4 Answers4

Linked