5

Hi I want to calculate errors in slope and intercept which are calculated by scipy.polyfit function. I have (+/-) uncertainty for ydata so how can I include it for calculating uncertainty into slope and intercept? My code is, 

from scipy import polyfit
import pylab as plt
from numpy import *

data = loadtxt("data.txt")
xdata,ydata = data[:,0],data[:,1]


x_d,y_d = log10(xdata),log10(ydata)
polycoef = polyfit(x_d, y_d, 1)
yfit = 10**( polycoef[0]*x_d+polycoef[1] )


plt.subplot(111)
plt.loglog(xdata,ydata,'.k',xdata,yfit,'-r')
plt.show()

Thanks a lot

bmu
  • 35,119
  • 13
  • 91
  • 108
physics_for_all
  • 2,193
  • 4
  • 19
  • 20

1 Answers1

4

You could use scipy.optimize.curve_fit instead of polyfit. It has a parameter sigma for errors of ydata. If you have your error for every y value in a sequence yerror (so that yerror has the same length as your y_d sequence) you can do:

polycoef, _ = scipy.optimize.curve_fit(lambda x, a, b: a*x+b, x_d, y_d, sigma=yerror)

For an alternative see the paragraph Fitting a power-law to data with errors in the Scipy Cookbook.

halex
  • 16,253
  • 5
  • 58
  • 67
  • Thanks for the reply. Yes I have seen that power law function but how can I combine my +/- error together with ydata? For example my ydata looks like, Y = 5 (+0.1, -0.4), 4.7 (+0.7,-0.4),..etc. – physics_for_all Oct 03 '12 at 12:46
  • @viralparekh You have an asymmetric deviation of your values? Never seen that before :) Could you elaborate a little on why the positive deviation is different from the negative. – halex Oct 03 '12 at 12:51
  • It shows confidence region. so 1st Y value is between (5.09 and 4.59) and so on. It just showing +ve(high) and -ve(low) error. – physics_for_all Oct 03 '12 at 12:59