有没有一种方法可以优化输出，限制来自 scipy.optimize 的 fmin 中的参数答案

【问题标题】：Is there a way where I can optimize the output restricting the parameters in fmin from scipy.optimize有没有一种方法可以优化输出，限制来自 scipy.optimize 的 fmin 中的参数
【发布时间】：2020-07-20 16:29:15
【问题描述】：

我在做什么：我修改了僵尸入侵系统的代码以演示它应该如何编写，并尝试使用 fmin 函数优化最小二乘误差（定义为分数函数）。

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

我对此很陌生，所以在正确方向上的任何帮助都会很棒。随时询问有关代码的任何疑问，我将尽我所知回答。 .谢谢。

【问题讨论】：

如果您使用minimize 而不是 fmin，您可以指定变量的界限。
您可以将惩罚项K*(max(0,x_lb-x)+max(0,x-x_ub)) 与大K 添加到得分函数。要获得更平滑的函数，请取这些项的平方。还要考虑障碍项，但这更明智，因为算法必须保持在允许值的范围内。
是的，但问题是，这是一个较小的模型，我正在处理冠状病毒数据，并且在定义最小化函数时遇到了一些问题
没有冠状病毒数据，只有关于阳性检测结果的数据。如果专家们在最近的一次采访中像 Streek 一样诚实，他们会告诉你，该测试只检测 RNA 片段，这并不能证明存在病毒。这忽略了测试本身是一个超过 40 次迭代的分形加倍和折叠过程，并且由于测试材料本身的污染和自我出现，每次测试在某种程度上都是阳性的。其他所有病毒检测方案都使用这种 RT-PCR 测试，就好像这些问题不存在一样，从而产生了循环论证。
样本中最大的RNA来源是被破坏的人体细胞，因此感染越严重，样本中污染的细胞碎片就越多。因此，该测试最有可能检测到的只是一小部分传染病，并且根据采样程序的设计，主要是呼吸道感染。这又偏向于更严重的病例，即处于寒冷期的年老体弱的人。不必在世界各地都存在一种特定的病毒，根本不需要病毒，细菌感染也可以，...

标签： python numpy scipy mathematical-optimization ode

【解决方案1】：

我不知道为什么Adventures in Python : Fitting a Differential Equation System to Data 的原作者会跳过箍来获取与给定数据点对应的样本，通过将时间数组而不是其构造传递给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)

这确实使所有参数都为正（以前在我的笔记本中，第一个参数是负的）。

【讨论】：