RまたはPythonを使用して非線形モデルにパラメータ値をどのように適合させるか

Question

私は関数と点集合を持っています。機能は次のとおりです。mとkは、問題のパラメータです

s(t) = m*g*t/k-(m*m/(k*k))*(exp(-k*t/m)-1) 

。 は、私は私のデータにこの式を合わせてm最高とkを見つける必要がある、と私は少しプログラミングの経験を持っている(t, s(t))

のデータを持っています。

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つは、より有望に見えるが、私は多分それの1/2を理解しています。私が理解していない部分は、それをプロットするときに私が見ているもの、それに私のモデル関数を与える方法、成功したこと、そして他の多くの問題です。

誰かがx、y、およびzをすれば、これは動作しますと言えば、私はそれを行うだろうが、それはより複雑だと、私は少し時間があるので、それ以外の場合は私はこの1つの明確な舵取りしたいと思います。

私のモデルの機能にこれらのパラメータを適合させるのを助けてください。私はとても混乱しています。

Answer 1

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 

Answer 2

あなたの方程式はmとkの比率によって決まります。すなわち、これは多くのソルバに問題を与える

r = m/k 

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

あるデータがないので、例えば、値M、Kと10×m個、10 * K区別する。

あなたはr単独で解決する方がよいでしょう。

Answer 3

@dmuirが指摘しているように、mとkは完全に相関しています。

この場合、python lmfitモジュール（http://lmfit.github.io/lmfit-py/）が役に立ちます。これは、Parameterオブジェクトのより良い抽象化を伴う、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() 

あなたのフィット感は、私には非常に良い見ていないが、これは、あなたがデータとモデルの機能を探索するより良いことができるように役立つかもしれません。