Nic*_*ick 5 python curve-fitting scipy
我希望y=x ** pw使用一个常数拟合指数函数pw来拟合两个数据点.该scipy curve_fit功能应优化adj1和adj2.我已尝试使用下面的代码,但无法使其工作.曲线不会通过数据点.我该如何解决?
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, adj1,adj2):
return np.round(((x+adj1) ** pw) * adj2, 2)
x = [0.5,0.85] # two given datapoints to which the exponential function with power pw should fit
y = [0.02,4]
pw=15
popt, pcov = curve_fit(func, x, y)
xf=np.linspace(0,1,50)
plt.figure()
plt.plot(x, y, 'ko', label="Original Data")
plt.plot(xf, func(xf, *popt), 'r-', label="Fitted Curve")
plt.show()
Run Code Online (Sandbox Code Playgroud)
如果您只想从两个数据点找到目标函数中的两个参数,这不一定是最小二乘拟合的问题:只需求解联立方程 y1 = b(x1+a)^p 和 y2 = b (x2+a)^p 为参数a和b:
import numpy as np
import matplotlib.pyplot as plt
def func(x, adj1,adj2):
return ((x+adj1) ** pw) * adj2
x = [0.5,0.85] # two given datapoints to which the exponential function with power pw should fit
y = [0.02,4]
pw = 15
A = np.exp(np.log(y[0]/y[1])/pw)
a = (x[0] - x[1]*A)/(A-1)
b = y[0]/(x[0]+a)**pw
xf=np.linspace(0,1,50)
plt.figure()
plt.plot(x, y, 'ko', label="Original Data")
plt.plot(xf, func(xf, a, b), 'r-', label="Fitted Curve")
plt.show()
Run Code Online (Sandbox Code Playgroud)
当然,结果函数将准确地通过这两个点。
这是解决方案。我认为对于曲线拟合 lmfit 是 scipy 的一个很好的替代品。
from lmfit import minimize, Parameters, Parameter, report_fit
import numpy as np
# create data to be fitted
xf = [0.5,0.85] # two given datapoints to which the exponential function with power pw should fit
yf = [0.02,4]
# define objective function: returns the array to be minimized
def fcn2min(params, x, data):
pw = params['pw'].value
adj1 = params['adj1'].value
adj2 = params['adj2'].value
model = adj1 * np.power(x + adj2, pw)
return model - data
pw=2
# create a set of Parameters
params = Parameters()
params.add('pw', value= pw, vary=False)
params.add('adj1', value= 1)
params.add('adj2', value= 1)
# do fit, here with leastsq model
result = minimize(fcn2min, params, args=(xf, yf))
# calculate final result
final = yf + result.residual
# write error report
report_fit(result.params)
adj1=result.params['adj1']
adj2=result.params['adj2']
# try to plot results
x = np.linspace(0, 1, 100)
y = adj1 * np.power(x + adj2, pw)
import pylab
pylab.plot(xf, yf, 'ko')
pylab.plot(x, y, 'r')
pylab.show()
Run Code Online (Sandbox Code Playgroud)
| 归档时间: |
|
| 查看次数: |
1708 次 |
| 最近记录: |