如何使用R或Python在非线性模型中拟合参数值


我有一个函数和一组点。功能是:
s(t) = m*g*t/k-(m*m/(k*k))*(exp(-k*t/m)-1)

其中CCD_ 1和CCD_。我有(t, s(t))的数据

我需要找到最好的mk来将这个等式与我的数据相匹配,而且我几乎没有编程经验。

R:中似乎更容易

df<-read.table("freefall2.txt", header = FALSE, sep = 't', 
strip.white = TRUE, na.strings = "empty")
time<-c(df[1])
position<-c(df[3])
fallData<-data.frame(t=time, position=c(df[3]))
plot(fallData)
m<-60
k<-1
secs=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
ft=c(16,62,138,  242,  366,  504,  652,  808,  971,
1138, 1309, 1483, 1657, 1831, 2005, 2179, 2353, 2527, 2701, 2875)
jitter(secs)
jitter(ft)
fit<-nls(ft~(m*g*secs/k)+(m^2/k^2)*(exp(-k*secs/m)-1),start = list(k=k,m=m))

这是最好的尝试,但它抛出:

Error in nlsModel(formula, mf, start, wts) : 
singular gradient matrix at initial parameter estimates

Python代码:

def func(t,m,k):
g=32.174
return m*g*t/k-(m*m/(k*k))*(np.exp(-k*t/m)-1)
def plotone():
plt.plot(D[0], D[2], 'b-', label='data')
popt, pcov = curve_fit(func, D[0], D[2])
plt.plot(D[0], func(D[0], *popt), 'r-', label='fit')
popt, pcov = curve_fit(func, D[0], D[2], bounds=(0, [120, 120]))
plt.plot(D[0], func(D[0], *popt), 'g--', label='fit-with-bounds')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

它生成了一个线性模型,没有返回所需的参数值。

在Python中重新启动了大量配方采样的尝试:

x = D[0]
y = D[2]
ycount = 0
fitfunc = lambda t, m, k: m*32.174*t/k-(m*m/(k*k))*(np.exp(-k*t/m)-1) # Target function
errfunc = lambda t, m, k: fitfunc(t, m, k) - y[ycount]; ycount+=1# Distance to the target function
p0 = [-15., 0.8, 0., -1.] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, D[2], args=(x, y))
print(type(p1))
print(type(success))
print(p1, success)

这个看起来更有前景,但我可能理解它的1/2。我不理解的部分是我绘制它时看到的是什么,如何为它提供我的模型函数,什么是成功,以及一系列其他问题。

如果有人说如果你做xym0,这会起作用,我会做的,但除此之外,我想避开这个,因为它更复杂,我几乎没有时间。

请帮助我将这些参数与我的模型函数相匹配。我很困惑。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
xdata2 = np.linspace(1,20,20)
xdata2
ydata2=([16,62,138,  242,  366,  504,  652,  808,  971,1138, 1309, 1483, 1657, 1831, 2005, 2179, 2353, 2527, 2701, 2875])
ydata2
def func2(t,m,k):
g=32.174
return m*g*t/k-(m*m/(k*k))*(np.exp(-k*t/m)-1)
plt.plot(xdata2, ydata2, 'b-', label='xdata2')
#plt.show()
#xdata2
#ydata2
popt2, pcov2 = curve_fit(func2, xdata2, ydata2)
plt.plot(xdata2, func2(xdata2, *popt2), 'r-', label='fit')
#plt.show()
#popt3, pcov3 = curve_fit(func2, xdata2, ydata2, bounds=(0, [3., 2.]))
#plt.plot(xdata2, func2(xdata2, *popt3), 'g--', label='fit-with-bounds')

plt.show()
pcov2
popt2
#popt3

你的s方程只取决于m和k的比值。那是

s(t) = r*g*t-(r*r)*(exp(-t/r)-1)

其中

r = m/k

这将给许多求解器带来问题,因为例如,数据无法区分值m、k和10*m、10*k。

你最好单独解决r。

正如@dmuir所指出的,mk是完全相关的。

您可能会发现pythonlmfit模块(http://lmfit.github.io/lmfit-py/)这一点很有用。这比具有更好的参数对象抽象的curve_fit稍高一点。您的示例如下:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
secs = np.linspace(1, 20, 20)
ft = np.array([16, 62, 138, 242, 366, 504, 652, 808, 971, 1138, 1309, 1483,
1657, 1831, 2005, 2179, 2353, 2527, 2701, 2875])
def func(t, m, k, g=32.174):
mk = m / k
return mk*g*t + (mk*mk)*(np.exp(-t/mk)-1)
model = Model(func)
# make parameters with initial values for parameters
params = model.make_params(m=-50, g=32.174, k=-1)
# you can now fix parameters or keep them varying in the fit:
params['g'].vary = False
params['m'].vary = False
# or set bounds on parameters (you might want to avoid divide by 0)
params['k'].max=-1.e-8
# now do the fit:
result = model.fit(ft, params, t=secs, nan_policy='omit')
# print results
print(result.fit_report())
# plot results:
plt.plot(secs, ft, '.', label='data' )
plt.plot(secs, result.best_fit, '--', label='fit')
plt.xlabel('secs')
plt.ylabel('ft')
plt.legend()
plt.show()

在我看来,你的贴合度不太好,但这可能有助于你更好地探索数据和模型功能。

相关内容

最新更新