Python:优化求解器返回非线性回归问题的初始猜测
Python:优化求解器返回非线性回归问题的初始猜测
提问于 2019-02-08 11:37:34
下面是用于常微分方程参数的最小二乘拟合的代码。使用了Python的“最小化”和“最小二乘”函数。已经尝试了不同的方法和ODE解算器/步骤(scipy ode/odeint)。这是一个在MATLAB中很容易解决的问题，但是Python总是返回最初的猜测。我希望你找出一个编码错误，否则我会对Python优化函数感到失望。Obj表示目标(剩余平方和)，ode函数(一阶)表示参数未知的方程。数据集已附加。
import numpy as np
from scipy.integrate import ode
from scipy.optimize import least_squares
from scipy.optimize import minimize
from scipy.optimize import SR1
import matplotlib.pyplot as plt
import math
Minput=np.loadtxt('C:\\Users\\Ladan\\Documents\\Python Scripts\\Python\\moisturesmoothopt.txt')
Minput=Minput.flatten()
time=np.linspace(0,1800,901)
A=np.zeros(3)
XC,RC,alpha=A
#bnds=([0,0,0],[Minput[0],math.inf,math.inf])
bnds=((0,Minput[0]),(0,math.inf),(0,math.inf))
def firstorder(X,time,A):
if X>=XC:
dX=-RC
if X<XC:
dX=-RC*(X/XC)**alpha
return dX
def obj(A):
X0=Minput[0]
# Xpred=odeint(firstorder,X0,time,args=(A,))
Xpred=ode(firstorder).set_integrator('vode', method='bdf',
order=15).set_initial_value(Minput[0],0).set_f_params(A)
#Xpred=ode(firstorder).set_integrator('lsoda').set_initial_value(Minput[0],0).set_f_params(A)
EPR=Xpred
EPR2=EPR.y.flatten()
ERRone=np.sum(np.power((EPR2-Minput),2))
ERR=ERRone/((901-3)) # residual sum of squares deivided by dof
return ERR
XC=1
RC=0.005
alpha=1.5
A0=[XC,RC,alpha]
Parameters=minimize(obj,A0,method='SLSQP',bounds=bnds,options={'ftol':1e-10,
'maxiter': 1000})
print('parameters',Parameters) Minput数组的数据在线共享：
回答 1
Stack Overflow用户
发布于 2019-02-08 12:15:54
虽然我在scipy中使用了较新的ode解算器，但我倾向于使用good ol‘odeint函数，它有一点旧，但仍然相当可靠，在许多情况下，由于我不完全理解的原因，它的性能更好。无论如何，我在很大程度上重构了您的分析代码，以便同时使用scipy.optimize.least_squares和scipy.integrate.odeint。虽然取得了一些进展，但我确实收到了关于power中的无效值的某种警告。您将不得不进一步研究，但这应该会让您走上正确的道路
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import least_squares
Minput=np.loadtxt('Data.txt').T
time=np.linspace(0,1800,901)
bnds=([0,0,0],[Minput[0],np.inf,np.inf])
def firstorder(X,time, XC, RC, alpha):
if X >= XC:
dX = -RC
else:
dX = -RC * (X/XC)**alpha
return dX
XC=1
RC=0.005
alpha=1.5
A0=(XC, RC, alpha)
def residuals(x0):
Mcalc = odeint(firstorder, y0=Minput[0], t=time, args=tuple(x0))[:,0]
return Mcalc - Minput
Mcalc = odeint(firstorder, y0=Minput[0], t=time, args=A0)[:,0]
result = least_squares(residuals, x0=A0, bounds=bnds)
print(result)
Mfit = odeint(firstorder, y0=Minput[0], t=time, args=tuple(result.x))[:,0]
plt.plot(time, Minput, label='data')
plt.plot(time, Mcalc, label='initial values')
plt.plot(time, Mfit, label='fit')
plt.legend()打印输出：
/---/python3.6/site-packages/ipykernel_launcher.py:20: RuntimeWarning: invalid value encountered in power
active_mask: array([0, 0, 0])
cost: 50.571520689865935
fun: array([ 0.00000000e+00, 8.14148814e-03, 1.61829763e-02, 2.44244644e-02,
...])
grad: array([-1.18907831, nan, -7.75389712])
jac: array([[ 0. , 0. , 0. ],
[ 0. , -2. , 0. ],
[ 0. , -3.99999994, 0. ],
...,
[ 0.07146036, nan, 0.1084222 ],
[ 0.07130456, nan, 0.10827456],
[ 0.07114924, nan, 0.1081272 ]])
message: '`xtol` termination condition is satisfied.'
nfev: 30
njev: 12
optimality: nan
status: 3
success: True
x: array([1.0000002 , 0.00717926, 1.50000032])
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/54585676
复制相关文章
相似问题
腾讯云开发者
