<p>你有两个选择:</p>
<ol>
<li>将系统线性化,并将一条直线拟合到数据日志中。</li>
<li>使用非线性解算器(例如<a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html" rel="noreferrer">^{<cd1>}</a></li>
</ol>
<p>到目前为止,第一个选择是最快、最强大的。但是,它要求你知道y偏移的先验,否则不可能线性化方程。(也就是说,<code>y = A * exp(K * t)</code>可以通过拟合<code>y = log(A * exp(K * t)) = K * t + log(A)</code>来线性化,但是<code>y = A*exp(K*t) + C</code>只能通过拟合<code>y - C = K*t + log(A)</code>来线性化,并且由于<code>y</code>是您的自变量,因此必须事先知道<code>C</code>这是一个线性系统。</p>
<p>如果使用非线性方法,则a)不能保证收敛并得到解,b)速度会慢得多,c)对参数的不确定性给出的估计要差得多,d)通常精度要低得多。然而,与线性反演相比,非线性方法有一个巨大的优势:它可以求解非线性方程组。在您的情况下,这意味着您不必事先知道<code>C</code>。</p>
<p>举个例子,我们用线性和非线性的方法来求解y=A*exp(K*t)和一些噪声数据:</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.optimize
def main():
# Actual parameters
A0, K0, C0 = 2.5, -4.0, 2.0
# Generate some data based on these
tmin, tmax = 0, 0.5
num = 20
t = np.linspace(tmin, tmax, num)
y = model_func(t, A0, K0, C0)
# Add noise
noisy_y = y + 0.5 * (np.random.random(num) - 0.5)
fig = plt.figure()
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
# Non-linear Fit
A, K, C = fit_exp_nonlinear(t, noisy_y)
fit_y = model_func(t, A, K, C)
plot(ax1, t, y, noisy_y, fit_y, (A0, K0, C0), (A, K, C0))
ax1.set_title('Non-linear Fit')
# Linear Fit (Note that we have to provide the y-offset ("C") value!!
A, K = fit_exp_linear(t, y, C0)
fit_y = model_func(t, A, K, C0)
plot(ax2, t, y, noisy_y, fit_y, (A0, K0, C0), (A, K, 0))
ax2.set_title('Linear Fit')
plt.show()
def model_func(t, A, K, C):
return A * np.exp(K * t) + C
def fit_exp_linear(t, y, C=0):
y = y - C
y = np.log(y)
K, A_log = np.polyfit(t, y, 1)
A = np.exp(A_log)
return A, K
def fit_exp_nonlinear(t, y):
opt_parms, parm_cov = sp.optimize.curve_fit(model_func, t, y, maxfev=1000)
A, K, C = opt_parms
return A, K, C
def plot(ax, t, y, noisy_y, fit_y, orig_parms, fit_parms):
A0, K0, C0 = orig_parms
A, K, C = fit_parms
ax.plot(t, y, 'k--',
label='Actual Function:\n $y = %0.2f e^{%0.2f t} + %0.2f$' % (A0, K0, C0))
ax.plot(t, fit_y, 'b-',
label='Fitted Function:\n $y = %0.2f e^{%0.2f t} + %0.2f$' % (A, K, C))
ax.plot(t, noisy_y, 'ro')
ax.legend(bbox_to_anchor=(1.05, 1.1), fancybox=True, shadow=True)
if __name__ == '__main__':
main()
</code></pre>
<p><img src="https://i.stack.imgur.com/S8r8z.png" alt="Fitting exp"/></p>
<p>注意,线性解提供的结果更接近实际值。但是,我们必须提供y偏移值才能使用线性解。非线性解不需要先验知识。</p>