107 lines
7.6 KiB
Python
Executable File
107 lines
7.6 KiB
Python
Executable File
###################################################################################################
|
|
# Fitting with offset using least square optimization
|
|
###################################################################################################
|
|
|
|
from mathutils import *
|
|
from plotutils import *
|
|
|
|
###################################################################################################
|
|
#Fitting the gaussian function with offset f(x) = a + b * exp(-(pow((x - c), 2) / (2 * pow(d, 2))))
|
|
###################################################################################################
|
|
|
|
xdata=[-200.30429237268825, -200.2650700434188, -200.22115208318002, -199.9457671375377, -199.86345548879072, -199.85213073174933, -199.35687977133284, -199.13811861090275, -197.97304970346386, -197.2952215624348, -195.09076092936948, -192.92276048970703, -191.96871876227698, -189.49577852322938, -187.9652790409825, -183.63756456925222, -180.04899765472996, -178.43839623242422, -174.07311671294445, -172.0410133577918, -165.90824309893102, -160.99771795989466, -159.30176653939253, -154.27688897558514, -152.0854103810786, -145.75652847587313, -140.80843828908465, -139.23982133191495, -134.27073891256106, -132.12649284133064, -125.95947209775511, -121.00309550337462, -119.26736932643232, -114.2706655484383, -112.07393889578914, -105.72295990367157, -100.8088439880125, -99.2034906238494, -94.30042325164636, -92.15010048151461, -85.92203653534293, -81.03913275494665, -79.27412793784428, -74.33487658582118, -72.06274362408762, -65.76562628131825, -60.91255356825276, -59.20334389560392, -54.33286972659312, -52.19387171350535, -45.94978737932291, -41.03014719193582, -39.301602568238906, -34.35572209014114, -32.04464301272608, -25.8221033382824, -20.922074315528747, -19.21590299233186, -14.31090212502093, -12.217203140101386, -5.9283722049240435, -0.9863587170369246, 0.7408048387279834, 5.71126832601389, 7.972628957879352, 14.204559894256546, 19.11839959633025, 20.8218087836657, 25.678748486941828, 27.822718344586864, 34.062659474970715, 38.9745656819391, 40.77409719734158, 45.72080631619803, 47.974156754056835, 54.23453768983539, 59.12020360609568, 60.77306570712026, 65.70734521458867, 67.8344660434617, 74.03187028154134, 78.96532114824849, 80.76070945985495, 85.74802197591286, 87.9140889204674, 94.18082276873524, 99.25790470037091, 100.68454787413205, 105.7213026221542, 107.79483801526698, 113.99555681638138, 119.0707052529143, 120.72715813056156, 125.77551384921307, 127.91257836719551, 134.2011330887875, 139.23043006997628, 140.71673537840158, 145.76288138835983, 147.80216629676042, 154.06420451405637, 159.0846626604798, 160.76183155710717, 165.73699067536242, 167.9265357747636, 173.96705069576544, 178.2522282751915, 179.9042617354548, 183.54586165856657, 185.23269803071796, 189.41678143751972, 191.87149157986588, 192.8741468985015, 195.0241934550453, 195.966634211846, 197.9821647518146, 198.99006812859284, 199.33202054855676, 199.91897441965887, 200.11536227958896, 200.22280936469997, 200.25181179127208]
|
|
ydata=[11.0, 6.0, 8.0, 5.0, 11.0, 7.0, 18.0, 11.0, 12.0, 10.0, 8.0, 6.0, 16.0, 4.0, 12.0, 9.0, 15.0, 14.0, 8.0, 20.0, 15.0, 8.0, 9.0, 11.0, 13.0, 12.0, 13.0, 15.0, 13.0, 20.0, 10.0, 7.0, 17.0, 11.0, 20.0, 13.0, 13.0, 23.0, 14.0, 10.0, 17.0, 15.0, 20.0, 16.0, 14.0, 13.0, 18.0, 22.0, 9.0, 20.0, 12.0, 14.0, 17.0, 19.0, 14.0, 14.0, 23.0, 19.0, 15.0, 20.0, 20.0, 21.0, 20.0, 23.0, 22.0, 15.0, 10.0, 17.0, 21.0, 15.0, 23.0, 23.0, 25.0, 18.0, 16.0, 21.0, 22.0, 16.0, 16.0, 14.0, 19.0, 20.0, 18.0, 20.0, 23.0, 13.0, 16.0, 20.0, 25.0, 15.0, 15.0, 17.0, 22.0, 26.0, 19.0, 30.0, 25.0, 17.0, 17.0, 23.0, 16.0, 27.0, 21.0, 21.0, 26.0, 27.0, 21.0, 17.0, 20.0, 20.0, 21.0, 19.0, 25.0, 19.0, 13.0, 23.0, 20.0, 20.0, 18.0, 20.0, 19.0, 25.0]
|
|
#ydata = [t-8.0 for t in ydata]
|
|
#xdata = [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0]
|
|
#ydata = [24.0, 36.0, 66.0, 121.0, 474.0, 989.0, 357.0, 175.0, 50.0, 40.0, 30.0, 22.0]
|
|
|
|
#xdata = [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
|
|
#ydata = [36.0, 66.0, 121.0, 183.0, 263.0, 365.0, 473.0, 603.0, 753.0, 917.0]
|
|
|
|
weigths = [ 1.0] * len(xdata)
|
|
|
|
p = plot(ydata, "Data" , xdata)[0]
|
|
p.setLegendVisible(True)
|
|
|
|
guess = None
|
|
def fit_gaussian_offset(y, x, start_point = None, weights = None):
|
|
"""Fits data into a gaussian with offset.
|
|
f(x) = a + b * exp(-(pow((x - c), 2) / (2 * pow(d, 2))))
|
|
|
|
Args:
|
|
x(float array or list): observed points x
|
|
y(float array or list): observed points y
|
|
start_point(optional tuple of float): initial parameters (normalization, mean, sigma)
|
|
weights(optional float array or list): weight for each observed point
|
|
Returns:
|
|
Tuples of gaussian parameters: (offset, normalization, mean, sigma)
|
|
"""
|
|
|
|
# For normalised gauss curve sigma=1/(amp*sqrt(2*pi))
|
|
if start_point is None:
|
|
off = min(y) # good enough starting point for offset
|
|
com = x[y.index(max(y))]
|
|
amp = max(y) - off
|
|
sigma = trapz([v-off for v in y], x) / (amp*math.sqrt(2*math.pi))
|
|
start_point = [off, amp, com , sigma]
|
|
print "Start point", start_point
|
|
|
|
class Model(MultivariateJacobianFunction):
|
|
def value(self, variables):
|
|
value = ArrayRealVector(len(x))
|
|
jacobian = Array2DRowRealMatrix(len(x), 4)
|
|
for i in range(len(x)):
|
|
(a,b,c,d) = (variables.getEntry(0), variables.getEntry(1), variables.getEntry(2), variables.getEntry(3))
|
|
v = math.exp(-(math.pow((x[i] - c), 2) / (2 * math.pow(d, 2))))
|
|
model = a + b * v
|
|
value.setEntry(i, model)
|
|
jacobian.setEntry(i, 0, 1) # derivative with respect to p0 = a
|
|
jacobian.setEntry(i, 1, v) # derivative with respect to p1 = b
|
|
v2 = b*v*((x[i] - c)/math.pow(d, 2))
|
|
jacobian.setEntry(i, 2, v2) # derivative with respect to p2 = c
|
|
jacobian.setEntry(i, 3, v2*(x[i] - c)/d ) # derivative with respect to p3 = d
|
|
return Pair(value, jacobian)
|
|
|
|
model = Model()
|
|
target = [v for v in y] #the target is to have all points at the positios
|
|
(parameters, residuals, rms, evals, iters) = optimize_least_squares(model, target, start_point, weigths)
|
|
print evals, iters
|
|
return parameters
|
|
|
|
(off, amp, com, sigma) = fit_gaussian_offset(ydata, xdata, None, weigths)
|
|
print "\nFit with Offset:\noff: ", off , " amp: ", amp, " com: ", com, " sigma: ", sigma
|
|
f = Gaussian(amp, com, sigma)
|
|
gauss = [f.value(i)+off for i in xdata]
|
|
s = LinePlotSeries("Fit with offset")
|
|
p.addSeries(s)
|
|
s.setData(xdata, gauss)
|
|
e=0
|
|
for i in range(len(ydata)) : e = e + abs(ydata[i]-gauss[i])
|
|
print "E = ", e
|
|
|
|
(amp, com, sigma) = fit_gaussian(ydata, xdata, None, weigths)
|
|
print "\nNormal Fit:\namp: ", amp, " com: ", com, " sigma: ", sigma
|
|
f = Gaussian(amp, com, sigma)
|
|
gauss = [f.value(i) for i in xdata]
|
|
s = LinePlotSeries("Normal Fit")
|
|
p.addSeries(s)
|
|
s.setData(xdata, gauss)
|
|
e=0
|
|
for i in range(len(ydata)) : e = e + abs(ydata[i]-gauss[i])
|
|
print "E = ", e
|
|
|
|
|
|
|
|
run("CPython/GaussFit_wrapper")
|
|
[off, amp, com, sigma] = gfitoff(xdata, ydata, off=None, amp=None, com=None, sigma=None)
|
|
print "\nCPython Fit:\n", off , " amp: ", amp, " com: ", com, " sigma: ", sigma
|
|
from mathutils import Gaussian
|
|
g = Gaussian(amp, com, sigma)
|
|
gauss = [g.value(i)+off for i in xdata]
|
|
s = LinePlotSeries("CPython Fit")
|
|
p.addSeries(s)
|
|
s.setData(xdata, gauss)
|
|
e=0
|
|
for i in range(len(ydata)) : e = e + abs(ydata[i]-gauss[i])
|
|
print "E = ", e
|