import numpy as np
import matplotlib.pyplot as plt
def func(x, a1, a2, a3):
return a1*np.exp(-(x-a2)**2/a3**2)
ndata = 100
xdata = np.linspace(1, ndata, ndata)
y = func(xdata, 20, 60, 15)
ydata = y
plt.plot(xdata, ydata, 'b-', label='data')
popt, pcov = curve_fit(func, xdata, ydata)
plt.plot(xdata, func(xdata, *popt), 'r-', label='fit')
plt.show()
from pprint import pprint
import scipy.linalg as linalg
def dfda1(x,a1,a2,a3):
return np.exp(-(x - a2) ** 2 / a3 ** 2 / 2)
def dfda2(x,a1,a2,a3):
return a1 * (x - a2) / a3 ** 2 * np.exp(-(x - a2) ** 2 / a3 ** 2 / 2)
def dfda3(x,a1,a2,a3):
return a1 * (x - a2) ** 2 / a3 ** 3 * np.exp(-(x - a2) ** 2 / a3 ** 2 / 2)
nparam = 3
guess1 = [10,50,10]
df=np.zeros([ndata])
for i in range(0,ndata):
dy = ydata[i]-func(xdata[i], *guess1)
df[i]=dy
#pprint(df)
Jac=np.zeros([ndata,nparam])
for i in range(0,ndata):
Jac[i,0] = dfda1(xdata[i], *guess1)
Jac[i,1] = dfda2(xdata[i], *guess1)
Jac[i,2] = dfda3(xdata[i], *guess1)
# pprint(Jac)
iJac = linalg.inv(np.dot(np.transpose(Jac),Jac))
# print(iJac)
Jdf = np.dot(np.transpose(Jac),df)
# pprint(Jdf)
guess1 = guess1 + np.dot(iJac, Jdf)
pprint(guess1)
plt.plot(xdata, ydata, 'b-', label='data')
popt, pcov = curve_fit(func, xdata, ydata)
plt.plot(xdata, func(xdata, *guess1), 'r-', label='fit')
plt.show()
array([ 19.866, 59.841, 14.878])