使用 scipy 拟合微分方程，得到“对象对于所需数组来说太深”

Question

我正在尝试将曲线拟合到微分方程。为了简单起见，我只是在这里做逻辑方程。我写了下面的代码，但我得到下面显示的错误。我不太确定自己做错了什么。

import numpy as np
import pandas as pd
import scipy.optimize as optim
from scipy.integrate import odeint

df_yeast = pd.DataFrame({'cd': [9.6, 18.3, 29., 47.2, 71.1, 119.1, 174.6, 257.3, 350.7, 441., 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9, 659.6], 'td': np.arange(18)})

N0 = 1
parsic = [5, 2]

def logistic_de(t, N, r, K):
    return r*N*(1 - N/K)

def logistic_solution(t, r, K):
    return odeint(logistic_de, N0, t, (r, K))

params, _ = optim.curve_fit(logistic_solution, df_yeast['td'], df_yeast['cd'], p0=parsic);

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
ValueError: object too deep for desired array

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-94-2a5a467cfa43> in <module>
----> 1 params, _ = optim.curve_fit(logistic_solution, df_yeast['td'], df_yeast['cd'], p0=parsic);

~/SageMath/local/lib/python3.9/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    782         # Remove full_output from kwargs, otherwise we're passing it in twice.
    783         return_full = kwargs.pop('full_output', False)
--> 784         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    785         popt, pcov, infodict, errmsg, ier = res
    786         ysize = len(infodict['fvec'])

~/SageMath/local/lib/python3.9/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    420         if maxfev == 0:
    421             maxfev = 200*(n + 1)
--> 422         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
    423                                  gtol, maxfev, epsfcn, factor, diag)
    424     else:

error: Result from function call is not a proper array of floats.

Answer 1

@hpaulj 指出了来自logistic_solution 的返回值形状的问题，并表明修复可以消除您报告的错误。

但是，代码中还有另一个问题。该问题不会产生错误，但会导致您的测试问题（逻辑微分方程）的解不正确。默认情况下，odeint 期望计算微分方程的函数的 t 参数是 second 参数。更改logistic_de 的前两个参数的顺序，或者将参数tfirst=True 添加到您对odeint 的调用中。第二个选项更好一些，因为如果您决定尝试使用该功能而不是odeint，它将允许您使用logistic_de 和scipy.integrate.solve_ivp。

Answer 2

logistic_solution 的示例运行产生 (18,1) 结果：

In [268]: logistic_solution(df_yeast['td'], *parsic)
Out[268]: 
array([[ 1.00000000e+00],
       [ 2.66666671e+00],
       [ 4.33333337e+00],
       [ 1.00000004e+00],
       [-1.23333333e+01],
       [-4.06666666e+01],
       [-8.90000000e+01],
       [-1.62333333e+02],
       [-2.65666667e+02],
       [-4.04000000e+02],
       [-5.82333333e+02],
       [-8.05666667e+02],
       [-1.07900000e+03],
       [-1.40733333e+03],
       [-1.79566667e+03],
       [-2.24900000e+03],
       [-2.77233333e+03],
       [-3.37066667e+03]])
In [269]: _.shape
Out[269]: (18, 1)

但y 的值是

In [281]: df_yeast['cd'].values.shape
Out[281]: (18,)

定义一个返回一维数组的替代函数：

In [282]: def foo(t,r,K):
     ...:     return logistic_solution(t,r,K).ravel()

这行得通：

In [283]: params, _ = optim.curve_fit(foo, df_yeast['td'], df_yeast['cd'], p0=parsic)
In [284]: params
Out[284]: array([16.65599815, 15.52779946])

测试params:

In [287]: logistic_solution(df_yeast['td'], *params)
Out[287]: 
array([[  1.        ],
       [  8.97044688],
       [ 31.45157847],
       [ 66.2980814 ],
       [111.36464226],
       [164.50594767],
       [223.5766842 ],
       [286.43153847],
       [350.92519706],
       [414.91234658],
       [476.24767362],
       [532.78586477],
       [582.38160664],
       [622.88958582],
       [652.1644889 ],
       [668.0610025 ],
       [668.43381319],
       [651.13760758]])
In [288]: df_yeast['cd'].values
Out[288]: 
array([  9.6,  18.3,  29. ,  47.2,  71.1, 119.1, 174.6, 257.3, 350.7,
       441. , 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9, 659.6])

too deep 在此上下文中表示 2d 数组，当它应该是 1d 时，以便与 ydata 进行比较

ydata : array_like
    The dependent data, a length M array - nominally ``f(xdata, ...)``.