Skip to content

Instantly share code, notes, and snippets.

@raggleton
Created March 23, 2015 09:40
Show Gist options
  • Save raggleton/3722da82cca3d40aa1a1 to your computer and use it in GitHub Desktop.
Save raggleton/3722da82cca3d40aa1a1 to your computer and use it in GitHub Desktop.
Show fn, fit it, show the fit. Not perfect. For exploring possible shapes of fns.
from collections import namedtuple
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, RadioButtons
from scipy.optimize import curve_fit
from itertools import izip
# 0 to 0.347
xpt = [20.9499,20.7179,21.4734,22.0312,22.9523,24.6784,26.0982,28.7842,31.0813,33.0672,35.4857,38.3296,40.9982,43.7918,46.1765,49.2252,52.6305,55.078,58.345,61.1288,63.9831,67.1981,70.5813,73.5319,76.9119,79.7721,82.9528,86.3631,89.1184,93.0097,96.0561,99.2763,103.673,105.919,109.948,112.495,116.37,119.104,122.172,126.262,129.863,133.374,135.358,140.215,143.906,146.964,151.518,152.608,157.312,161.298,163.645,168.758,171.285,174.861,179.208,182.253,186.099,189.914,193.942]
ypt = [0.903289,1.31444,1.48306,1.57877,1.67689,1.75208,1.81371,1.78577,1.78443,1.77425,1.76701,1.74634,1.72242,1.68689,1.67924,1.66272,1.63702,1.62118,1.59681,1.5808,1.57993,1.55617,1.54106,1.53675,1.52014,1.50499,1.49345,1.48556,1.46981,1.46973,1.46255,1.45699,1.43602,1.42543,1.42378,1.42657,1.4099,1.40933,1.39935,1.38712,1.38268,1.36957,1.37589,1.35065,1.34913,1.34718,1.34197,1.35128,1.33667,1.33335,1.33269,1.32322,1.32131,1.3194,1.32064,1.31893,1.30569,1.30823,1.2943]
errx = [0.196885,0.181703,0.188646,0.168748,0.152036,0.172042,0.159209,0.184003,0.187899,0.174821,0.180278,0.196163,0.211693,0.235332,0.198388,0.222302,0.262741,0.238318,0.260262,0.261012,0.269579,0.26719,0.296326,0.297783,0.315932,0.311536,0.338922,0.367324,0.348821,0.370203,0.355142,0.374982,0.414767,0.408017,0.421579,0.435536,0.477844,0.440446,0.460131,0.477675,0.517247,0.467203,0.50454,0.50201,0.510912,0.539791,0.591163,0.580386,0.628283,0.606503,0.613952,0.688128,0.655641,0.620509,0.714075,0.740542,0.739552,0.787039,0.78965]
erry = [0.00404186,0.00288825,0.00281867,0.00187371,0.00195666,0.0019489,0.00252868,0.00210502,0.00233906,0.00241006,0.00275737,0.00293676,0.0025061,0.00233384,0.00267519,0.00265185,0.00262453,0.00284371,0.00259042,0.00258948,0.00274044,0.00273967,0.00306156,0.00263784,0.00308559,0.0028993,0.00287578,0.00342309,0.00327168,0.00346894,0.00312061,0.00338969,0.00305964,0.00282609,0.00440922,0.0035585,0.00316558,0.00326939,0.00289818,0.00296051,0.00306868,0.00337324,0.00334857,0.00304836,0.00319062,0.00352735,0.00329275,0.00323298,0.00328902,0.00337987,0.00346634,0.00285199,0.00375651,0.00360116,0.00341235,0.00342422,0.00330265,0.0035793,0.00312612]
# 2.172 to 3
# xpt = [22.2085,23.2103,24.4791,26.2364,28.392,31.0028,33.628,36.5432,39.5458,42.456,45.5817,49.0033,52.2821,55.545,58.8706,62.3055,65.6127,69.2169,72.57,76.1445,79.3819,82.6008,86.3982,89.783,93.4862,96.3865,100.25,104.034,107.642,110.441,114.215,118.109,121.636,125.134,128.132,132.763,135.596,140.26,143.363,146.553,150.083,154.892,157.734,161.857,164.186,168.584,172.949,177.078,179.632,183.31,187.095,191.937,193.011,198.364,200.809,204.657,208.971,212.259,216.375]
# ypt = [0.728784,0.87024,0.986051,1.32806,1.36371,1.35797,1.35721,1.35317,1.35284,1.33256,1.33445,1.31841,1.29611,1.3012,1.29829,1.29689,1.2801,1.27344,1.28257,1.27418,1.25238,1.25729,1.25045,1.24194,1.23692,1.22859,1.23409,1.2317,1.21002,1.22504,1.22018,1.21706,1.20166,1.20822,1.21683,1.19282,1.20538,1.19287,1.19148,1.18872,1.18706,1.19123,1.17653,1.17874,1.19442,1.17511,1.17854,1.16684,1.17282,1.17633,1.16966,1.16956,1.17182,1.17214,1.16891,1.1571,1.1598,1.15728,1.17224]
# errx = [.0493529,0.0580339,0.0630995,0.0701758,0.0773509,0.0917409,0.0971733,0.10551,0.121573,0.124725,0.133266,0.145182,0.153974,0.17171,0.17882,0.187128,0.199434,0.205396,0.226151,0.250594,0.233843,0.258088,0.267883,0.274999,0.280288,0.293742,0.30805,0.342408,0.333085,0.336997,0.370125,0.385893,0.399284,0.415534,0.426522,0.423709,0.456023,0.481854,0.494917,0.499839,0.514572,0.555277,0.545515,0.611515,0.57265,0.583315,0.659813,0.629008,0.659342,0.713475,0.680735,0.761322,0.751734,0.82094,0.774906,0.779585,0.842794,0.844081,0.88837]
# erry = [0.00228162,0.00230889,0.00223459,0.00729482,0.00596361,0.0053997,0.00471368,0.00557473,0.00505018,0.00454245,0.00506027,0.00450069,0.00388773,0.00409387,0.00443314,0.00482938,0.0041626,0.00461252,0.00578856,0.00518956,0.0036419,0.00395935,0.00389911,0.00451444,0.00416828,0.00378463,0.00461308,0.00438594,0.00345439,0.00481211,0.00391957,0.00493201,0.00369832,0.0050391,0.00519069,0.00368376,0.00383331,0.00472744,0.00405363,0.00374567,0.00392859,0.00392305,0.00421712,0.00452083,0.00499394,0.00391656,0.00445568,0.00394982,0.00545924,0.00476066,0.00480687,0.00379121,0.0049103,0.00484042,0.00414844,0.004481,0.00381101,0.00422995,0.00727298]
# 3 to 3.5
# xpt = [24.7695,25.6412,27.2112,29.4488,32.0511,34.1513,37.1472,39.591,42.3752,49.8829,63.1688,76.7433,91.1405,104.476,117.279,129.283,140.696,147.795,158.737]
# ypt = [0.749784,0.929346,1.03696,1.10314,1.13075,1.19174,1.22041,1.22072,1.26186,1.26897,1.31677,1.34706,1.3538,1.36199,1.37682,1.38985,1.41726,1.47868,1.55616]
# errx = [0.0952339,0.0979564,0.115596,0.114342,0.131148,0.141657,0.167114,0.178358,0.202079,0.109623,0.154173,0.200008,0.250445,0.321077,0.395747,0.491903,0.592245,0.798579,1.21704]
# erry = [0.019503,0.0196896,0.0157077,0.0162616,0.0134304,0.0151752,0.013843,0.00987511,0.0137001,0.00459906,0.00429588,0.00477587,0.00417857,0.00388191,0.00555698,0.00380698,0.00572628,0.00489103,0.0119748]
# 3.5 to 4
# xpt = [27.124,28.618,30.6127,32.9814,35.2829,37.6925,40.5013,42.9939,45.7114,52.9045,66.4928,79.6829,93.8372,106.698,120.09,129.572,139.078,144.074,157.435]
# ypt = [0.642894,0.763023,0.847863,0.913882,0.96903,1.01952,1.04279,1.08258,1.09688,1.15782,1.21956,1.26237,1.28311,1.32475,1.32287,1.37975,1.38819,1.44214,1.51371]
# errx = [0.0455256,0.0573744,0.0694312,0.0824988,0.0937487,0.106206,0.120707,0.135547,0.145856,0.0895192,0.132787,0.182127,0.262153,0.354762,0.491251,0.723062,1.21425,2.41203,2.03903]
# erry = [0.00900841,0.00976107,0.00825952,0.00952726,0.0105273,0.0101546,0.00812626,0.00950969,0.00740393,0.00405225,0.00334315,0.00349899,0.00329159,0.00565939,0.00617425,0.00425231,0.00837805,0.00815123,0.00884184]
# 4 to 4.5
# xpt = [31.5412,33.149,35.3766,38.0688,40.5904,43.0464,45.704,48.7603,51.292,58.2896,71.7278,84.4046,98.5355,108.665,123.669,128.611,162.25]
# ypt = [0.546626,0.645736,0.744541,0.784184,0.845308,0.888975,0.919132,0.932815,0.965419,1.03987,1.12093,1.1796,1.20463,1.25421,1.23498,1.25091,1.25]
# errx = [0.0618204,0.0791839,0.0979284,0.118723,0.137205,0.155045,0.175677,0.203366,0.223156,0.140967,0.235191,0.391987,0.698247,1.51816,1.73453,11.3738,4.06586]
# erry = [0.01323,0.0120326,0.0181955,0.0126068,0.0190929,0.0132642,0.0114843,0.0107492,0.00930005,0.00562028,0.00727727,0.00993301,0.0124697,0.0155407,0.0740973,0.01089,0.0117851]
# 4.5 to 5
# xpt = [22.623,24.323,26.3178,28.7714,30.8779,33.467,35.9519,38.3701,41.1725,47.4329,61.4801,73.275,90.75]
# ypt = [0.776655,0.902537,0.98518,1.02147,1.08243,1.10017,1.1326,1.16492,1.18513,1.23035,1.2431,1.37815,1.39535]
# errx = [0.0588717,0.0782651,0.0988073,0.12254,0.154739,0.181288,0.210894,0.276969,0.31735,0.252091,0.650735,1.9611,4.77297]
# erry = [0.0129404,0.0135614,0.0157836,0.0117613,0.0126293,0.0106766,0.00905822,0.00923386,0.0129776,0.00953222,0.00989179,1.47283,0.0471405]
print len(xpt)
print len(ypt)
print len(errx)
print len(erry)
# Draw the plot
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=0.36, top=0.95)
plt.errorbar(xpt, ypt, xerr=errx, yerr=erry, fmt='-')
# fitting function
def pf_func(et, p0, p1, p2, p3, p4, p5):
return p0 + (p1/(np.power(np.log10(et), 2)+p2)) + p3 * np.exp(-1.*p4*np.power(np.log10(et)-p5, 2))
# return p0*np.exp(np.divide(p1, et) + p2 + (p3*et) + (p4*np.power(et,2)) + (p5*np.power(et, 3)))
# PF
# central
# p0,p1,p2,p3,p4,p5 = [9.91227E-01, 9.16885, 1.47522, -1.21928E+03, 1.15238E-02, -2.24305E+01]
p0,p1,p2,p3,p4,p5 = [0.013411830610674243, 13.663922866277359, 5.4896362671416199, 10384.0, 2.0, -19.949999999999999]
# flat fwd
# p0,p1,p2,p3,p4,p5 = [1.516, -0.018, -1.722, 1.025E4, 0.84, -29]
# slight upwards fwd
# p0,p1,p2,p3,p4,p5 = [1.38, -0.02, -1.72, 10250, 0.84, 6.01]
# fwd one
# p0,p1,p2,p3,p4,p5 = [1.85, -2.12, 0.67, 10384, 2, -19.95]
# p0,p1,p2,p3,p4,p5 = [1., 1., 1., 1., 1., 1.]
# from fit
# p0,p1,p2,p3,p4,p5 = [ 1.53631577e+00, -4.19291248e-01, -1.29101049e+00, 1.03840000e+04, 2.00000000e+00, -1.99500000e+01]
# p0,p1,p2,p3,p4,p5 = [1.53631606e+00, -4.19292252e-01, -1.29100884e+00, 6.88495735e+05, 1.18894155e+00, 1.53447224e+01]
# et values
etmin = 00
etmax = xpt[-1]+2*errx[-1]
et = np.arange(etmin, etmax, 0.5)
lfit, = plt.plot(et, pf_func(et, p0, p1, p2, p3, p4, p5), lw=2, color='red')
ymin, ymax = 0.5, 2.5
plt.axis([etmin, etmax, ymin, ymax])
plt.xlabel("ET")
plt.ylabel("Correction Factor")
plt.grid(True)
axcolor = 'lightgoldenrodyellow'
ax0 = plt.axes([0.12, 0.24, 0.6, 0.03], axisbg=axcolor)
ax1 = plt.axes([0.12, 0.2, 0.6, 0.03], axisbg=axcolor)
ax2 = plt.axes([0.12, 0.16, 0.6, 0.03], axisbg=axcolor)
ax3 = plt.axes([0.12, 0.12, 0.6, 0.03], axisbg=axcolor)
ax4 = plt.axes([0.12, 0.08, 0.6, 0.03], axisbg=axcolor)
ax5 = plt.axes([0.12, 0.04, 0.6, 0.03], axisbg=axcolor)
# add sliders
s0 = Slider(ax0, 'p0', -5, 5, valinit=p0)
s1 = Slider(ax1, 'p1', -20, 20, valinit=p1)
s2 = Slider(ax2, 'p2', -10, 10, valinit=p2)
s3 = Slider(ax3, 'p3', -20000, 20000, valinit=p3)
s4 = Slider(ax4, 'p4', -5, 5, valinit=p4)
s5 = Slider(ax5, 'p5', -50, 50, valinit=p5)
# behaviour when changed
def update(val):
p0 = s0.val
p1 = s1.val
p2 = s2.val
p3 = s3.val
p4 = s4.val
p5 = s5.val
lfit.set_ydata(pf_func(et, p0, p1, p2, p3, p4, p5))
fig.canvas.draw_idle()
s0.on_changed(update)
s1.on_changed(update)
s2.on_changed(update)
s3.on_changed(update)
s4.on_changed(update)
s5.on_changed(update)
# add reset button
resetax = plt.axes([0.85, 0.025, 0.14, 0.04])
reset_button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')
def reset(event):
s0.reset()
s1.reset()
s2.reset()
s3.reset()
s4.reset()
s5.reset()
fig.canvas.draw_idle()
reset_button.on_clicked(reset)
# add print point button
printax = plt.axes([0.85, 0.09, 0.14, 0.04])
print_button = Button(printax, 'Print', color=axcolor, hovercolor='0.975')
def print_params(event):
for i, p in enumerate([p0, p1, p2, p3, p4, p5]):
print "p%d = %g" % (i,p)
print_button.on_clicked(print_params)
# fitting
minpt = 6
popt, pcov = curve_fit(pf_func, xpt[minpt:], ypt[minpt:], p0=[p0,p1,p2,p3,p4,p5], sigma=erry[minpt:])
print list(popt)
showax = plt.axes([0.85, 0.15, 0.14, 0.08])
show_button = Button(showax, 'Show fitted\n fn', color=axcolor, hovercolor='0.975')
def show_fit(event):
for s,p in izip([s0,s1,s2,s3,s4,s5], list(popt)):
s.set_val(p)
for par,p in izip([p0,p1,p2,p3,p4,p5], list(popt)):
par = p
show_button.on_clicked(show_fit)
#
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment