Есть ли способ, где я могу оптимизировать вывод, ограничивая параметры в fmin из scipy.optimize - PullRequest
Новая
       51

Есть ли способ, где я могу оптимизировать вывод, ограничивая параметры в fmin из scipy.optimize

0 голосов
/ 08 апреля 2020

Что я делаю: Я изменил код из системы вторжения zomb ie, чтобы продемонстрировать, как он должен быть написан, и попытался оптимизировать ошибку наименьших квадратов (определенную как функцию оценки) с помощью fmin function. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import integrate
from scipy.optimize import fmin
#=====================================================
#Notice we must import the Model Definition
from zombiewithdata import eq
#=====================================================

#1.Get Data
#====================================================
Td=np.array([0.5,1,1.5,2,2.2,3,3.5,4,4.5,5])#time
Zd=np.array([0,2,2,5,2,10,15,50,250,400])#zombie pop
#====================================================

#2.Set up Info for Model System
#===================================================
# model parameters
#----------------------------------------------------
P = 0       # birth rate
d = 0.0001  # natural death percent (per day)
B = 0.0095  # transmission percent  (per day)
G = 0.0001  # resurect percent (per day)
A = 0.0001  # destroy perecent (per day)
rates=(P,d,B,G,A)

# model initial conditions
#---------------------------------------------------
S0 = 500.               # initial population
Z0 = 0                  # initial zombie population
R0 = 0                  # initial death population
y0 = [S0, Z0, R0]      # initial condition vector

# model steps
#---------------------------------------------------
start_time=0.0
end_time=5.0
intervals=1000
mt=np.linspace(start_time,end_time,intervals)

# model index to compare to data
#----------------------------------------------------
findindex=lambda x:np.where(mt>=x)[0][0]
mindex=map(findindex,Td)
#=======================================================



#3.Score Fit of System
#=========================================================
def score(parms):
    #a.Get Solution to system
    F0,F1,F2,T=eq(parms,y0,start_time,end_time,intervals)
    #b.Pick of Model Points to Compare
    Zm=F1[mindex]
    #c.Score Difference between model and data points
    ss=lambda data,model:((data-model)**2).sum()
    return ss(Zd,Zm)
#========================================================


#4.Optimize Fit
#=======================================================
fit_score=score(rates)
answ=fmin(score,(rates),full_output=1,maxiter=1000000)
bestrates=answ[0]
bestscore=answ[1]
P,d,B,G,A=answ[0]
newrates=(P,d,B,G,A)
#=======================================================

#5.Generate Solution to System
#=======================================================
F0,F1,F2,T=eq(newrates,y0,start_time,end_time,intervals)
Zm=F1[mindex]
Tm=T[mindex]
#======================================================

Теперь в разделе #optimize fit можно ли как-нибудь получить наилучшие возможные значения скоростей, когда я ограничиваю значения «норм», таких как lb <= P, d, B, G, A <= ub, где lb = нижняя граница и ub = верхняя граница, и удалось получить минимум баллов в этой ограниченной области? Это не должно быть наиболее оптимизированным значением. fmin использует алгоритм Нелдера-Мида (симплекс). </p>

Я новичок в этом, поэтому любая помощь в правильном направлении была бы потрясающей. Не стесняйтесь задавать любые сомнения относительно кода, и я отвечу, насколько мне известно. , Спасибо.

1 Ответ

0 голосов
/ 09 апреля 2020

Я не уверен, почему первоначальный автор Приключения в Python: Подгонка системы дифференциальных уравнений к данным перепрыгивает через обручи, чтобы получить выборки, соответствующие заданным точкам данных, процедура может можно значительно упростить, передав eq массив времени вместо его параметров построения 

#=======================================================
def eq(par,initial_cond,t):
     #differential-eq-system----------------------
     def funct(y,t):
        Si, Zi, Ri=y
        P,d,B,G,A=par
        # the model equations (see Munz et al. 2009)
        f0 = P - B*Si*Zi - d*Si
        f1 = B*Si*Zi + G*Ri - A*Si*Zi
        f2 = d*Si + A*Si*Zi - G*Ri
        return [f0, f1, f2]
     #integrate------------------------------------
     ds = odeint(funct,initial_cond,t)
     return ds.T
#=======================================================

Это можно затем назвать 

T = np.linspace(0, 5.0, 1000+1)
S,Z,R=eq(rates,y0,T)

, но также с целью получения только значения, необходимые в функции score 

Tm=np.append([0],Td)
Sm,Zm,Rm=eq(rates,y0,Tm)

Это затем упрощает функцию оценки до 

def score(parms):
    #a.Get Solution to system
    Sm,Zm,Rm=eq(parms,y0,Tm)
    #c.Score Difference between model and data points
    ss=lambda data,model:((data-model)**2).sum()
    return ss(Zd,Zm[1:])

Теперь, если вы хотите, например, категорически отвергать отрицательные параметры, вы можете изменить возврат значение 

    return ss(Zd,Zm[1:]) + 1e6*sum(max(0,-x)**2 for x in parms)

, которое действительно делает все параметры положительными (где ранее в моей записной книжке был отрицательный первый параметр).

...