Revision history [back]
 | 1 | initial version |
I'm not sure what you mean by "with considering the errorbars" (which isn't quite English.) Do you mean merely that you want the covariance of the fit (which leastsq returns), or that your data has errors in y (the usual model) and you want that to affect the fitting process, or that your data has errors in both x and y and you want that to affect the fit?
Anyway, here are two relatively straightforward ways to fit and get covariances.
sage: import numpy as np
sage: import scipy.optimize
sage: import scipy.odr
sage: import random
sage: random.seed(3)
sage:
sage: #
sage: # generate the data
sage: #
sage:
sage: x = np.linspace(0, 10, 100)
sage: y = np.sin(0.7*x)+2.4*x # underlying function
sage: y += np.random.randn(len(x))/10 # add some noise
sage: x_err = [0.1]*len(x) # estimated errors
sage: y_err = [0.3]*len(x) # estimated errors
sage:
sage: # first, we can use curve_fit
sage:
sage: def fit1(x, a, b):
....: return np.sin(x*a)+b*x
....:
sage: # note: can't handle errors in x
sage: popt, pcov = scipy.optimize.curve_fit(fit1, x, y, sigma=y_err)
sage: popt
array([ 0.70124584, 2.39768117])
sage: pcov
array([[ 6.97634001e-06, -2.58153187e-06],
[ -2.58153187e-06, 4.22745005e-06]])
sage:
sage: # or use the more powerful odr approach, which unfortunately has
sage: # some different conventions
sage:
sage: def fit2(B, x):
....: return np.sin(x*B[0])+B[1]*x
....:
sage: mymodel = scipy.odr.Model(fit2)
sage: mydata = scipy.odr.RealData(x, y, sx=x_err, sy=y_err)
sage: myodr = scipy.odr.ODR(mydata, mymodel, beta0=[1., 1.], maxit=1000)
sage: myoutput = myodr.run()
sage: myoutput.pprint()
Beta: [ 0.70139136 2.39769177]
Beta Std Error: [ 0.00258362 0.00198363]
Beta Covariance: [[ 9.08706871e-05 -2.35028948e-05]
[ -2.35028948e-05 5.35658971e-05]]
Residual Variance: 0.0734568446584
Inverse Condition #: 0.679929358499
Reason(s) for Halting:
Numerical error detected
| | 2 | No.2 Revision |
I'm not sure what you mean by "with considering the errorbars" (which isn't quite English.) errorbars". Do you mean merely that you want the covariance of the fit (which leastsq returns), or that your data has errors in y (the usual model) and you want that to affect the fitting process, or that your data has errors in both x and y and you want that to affect the fit?
Anyway, here are two relatively straightforward ways to fit and get covariances.
sage: import numpy as np
sage: import scipy.optimize
sage: import scipy.odr
sage: import random
sage: random.seed(3)
sage:
sage: #
sage: # generate the data
sage: #
sage:
sage: x = np.linspace(0, 10, 100)
sage: y = np.sin(0.7*x)+2.4*x # underlying function
sage: y += np.random.randn(len(x))/10 # add some noise
sage: x_err = [0.1]*len(x) # estimated errors
sage: y_err = [0.3]*len(x) # estimated errors
sage:
sage: # first, we can use curve_fit
sage:
sage: def fit1(x, a, b):
....: return np.sin(x*a)+b*x
....:
sage: # note: can't handle errors in x
sage: popt, pcov = scipy.optimize.curve_fit(fit1, x, y, sigma=y_err)
sage: popt
array([ 0.70124584, 2.39768117])
sage: pcov
array([[ 6.97634001e-06, -2.58153187e-06],
[ -2.58153187e-06, 4.22745005e-06]])
sage:
sage: # or use the more powerful odr approach, which unfortunately has
sage: # some different conventions
sage:
sage: def fit2(B, x):
....: return np.sin(x*B[0])+B[1]*x
....:
sage: mymodel = scipy.odr.Model(fit2)
sage: mydata = scipy.odr.RealData(x, y, sx=x_err, sy=y_err)
sage: myodr = scipy.odr.ODR(mydata, mymodel, beta0=[1., 1.], maxit=1000)
sage: myoutput = myodr.run()
sage: myoutput.pprint()
Beta: [ 0.70139136 2.39769177]
Beta Std Error: [ 0.00258362 0.00198363]
Beta Covariance: [[ 9.08706871e-05 -2.35028948e-05]
[ -2.35028948e-05 5.35658971e-05]]
Residual Variance: 0.0734568446584
Inverse Condition #: 0.679929358499
Reason(s) for Halting:
Numerical error detected
