Created
March 23, 2015 09:40
-
-
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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