In [3]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
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()
In [4]:
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)
In [16]:
nparam = 3
guess1 = [10,50,10]
In [17]:
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([12.65890851, 65.30826287, 20.88231001])
これを4回ぐらい繰り返す．
guess1 = [10,50,10]
array([16.69058876, 55.85037666, 14.09363923])
array([19.03931932, 61.10281187, 15.0516337 ])
array([19.86610657, 59.84145965, 14.87841873])
array([20.0354474 , 60.01550602, 14.92040764])
下のプロットとはずれているかも．．．．
In [ ]: