Python: решатели оптимизации возвращают начальное предположение для задачи нелинейной регрессии - PullRequest
Новая
       14

Python: решатели оптимизации возвращают начальное предположение для задачи нелинейной регрессии

0 голосов
/ 08 февраля 2019

Ниже приведен код для подгонки параметров методом наименьших квадратов для ODE.Были использованы функции «минимизации» и «наименьших квадратов» Python.Были опробованы различные методы и решения / шаги ODE (scipy ode / odeint).Эта проблема легко решается в MATLAB, но Python продолжает возвращать первоначальное предположение.Я надеюсь, что вы найдете ошибку кодирования, или я буду разочарован функциями оптимизации Python.Obj показывает цель (остаточная сумма квадратов), а функция Оде (firstorder) показывает уравнения с неизвестными параметрами.Набор данных прилагается. 

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 доступны онлайн:

https://1drv.ms/t/s!AoVu1vtlAOiLasJxR7rzubDr8YE

1 Ответ

0 голосов
/ 08 февраля 2019

Несмотря на то, что я использовал более новые решатели од в scipy, я склоняюсь к хорошей старой функции odeint, которая немного старше, но все же довольно солидная и во многих случаях дает лучшую производительность по причинам, которые я не делаюполностью понимаю.В любом случае, я в значительной степени реструктурировал ваш код анализа, чтобы использовать как scipy.optimize.least_squares, так и scipy.integrate.odeint.Прогресс достигнут, но я получаю какое-то предупреждение о недопустимых значениях во власти.Вам нужно будет продолжить расследование, но это должно привести вас на правильный путь 

%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])

enter image description here

...