0

I am fitting some data I have with a function. With the fitted function I want to see where it crosses some threshold value, determine this threshold value and get an error on this value as well. When using the curve_fit toolbox with the absolute_sigma = True, I get a large error in the threshold value. When I turn it off, I get an extremely small error in the threshold value, both which seem not very realistic to me. 

#The numbers I use for the fit
absaverage = 1.0979598702453246

Nlist = [0.31974162, 0.52558248, 0.77172549, 1.34829036, 1.91809528, 3.08342098, 5.33183816, 6.60457399, 5.93014992]

averagelist = [0.294913849040662, 0.4648514538342791,0.6496807339899529,1.014641822518085,1.2981207560776595,1.703857575892128,2.0964571453613123,2.1999054339799295,2.151973289007941]


#%%Now fitting the threshold data to obtain <threshold value>
def fittie(x, a, b):
    return a*(1-np.exp(-b*x))

#guess values for the fit
aguess = 2.3
bguess = 0.42


popt,pcov = curve_fit(fittie, Nlist, averagelist, p0=[aguess,bguess], maxfev = 20000)

Nlistforplot = np.arange(0.1, 7, 0.05)


f=plt.figure()
plt.plot(Nlist, averagelist, 'bo', label='data',  markersize = 5)
plt.plot(Nlistforplot, fittie(Nlistforplot, *popt), 'r--', label='fit', linewidth = 1)
plt.axhline(y = absaverage, color = 'black', lw = '1')
plt.axvline(1.54, color = 'black', lw = '1')
plt.fill_between([0,12], absaverage, absaverage+max(averagelist)+100, color = 'blue', alpha = 0.15)
plt.ylabel('-Value',fontsize='xx-large')
plt.xlabel('Nlist',fontsize='xx-large')
plt.yscale('log')
plt.xscale('log')
plt.xlim(0.2,max(Nlist)+1)
plt.ylim(0.15,max(averagelist)+1)
plt.title('Threshold determination')
plt.legend()
plt.show()

afromfit = popt[0]
bfromfit = popt[1]




print "your threshold value is"
thresholdvalue = -(1/bfromfit)*np.log(1-(absaverage/afromfit))
print thresholdvalue

#ERROR IN THRESHOLD PROPAGATION
dfdb = (1/bfromfit**2)*np.log(1-(absaverage/afromfit))
dfda = -(1/bfromfit)*(1/(1-absaverage/afromfit))*(absaverage/(afromfit**2))

sigmax2 = dfda**2*pcov[0,0]+dfdb**2*pcov[1,1]+2*dfda*dfdb*pcov[1,0]
print "sigma in threshold value is"
print sigmax2**0.5

So I obtain the same threshold value of about 1.50. The errors seem completely off though, either too large or too small. Any idea?

Jaco Geuchies
  • 1
  • Add this line of code to supply the missing import statements: "import numpy as np;import matplotlib.pyplot as plt;from scipy.optimize import curve_fit" – James Phillips Sep 10 '19 at 15:56
  • I do not know where the data for `averagelist` comes from, but the residuals for the best fit are of order `1e-9`. So the super small `pcov` seems reasonable. – mikuszefski Sep 11 '19 at 14:05

0 Answers0