1

如何允许函数的一个参数在组之间变化而其他参数适合所有组？

我正在使用 lmfit 来拟合疾病传播的模型。我希望函数的指数适合所有数据点，但是比例因子需要在不同组之间有所不同（作为不同年份繁殖率不同的疾病的代表）。

请参阅下面的代码：

#### create parameters ####
params = Parameters()
params.add('tau_1', value=1.,min=0.01)
params.add('tau_2', value=1.,min=0.01)
params.add('rho_1', value=1.,min=0.01)
for i in range(0, len(sorted(variable_data.flu_season.unique()))):
    season = str(sorted(variable_data.flu_season.unique())[i])
    params.add('theta_' + season, value=1., min = 0.01)

#### model ####
def pop_dist_inverse_grav(params, dist, host_pop, targ_pop, flu_sea, data):

    parvals = params.valuesdict()
    tau_1 = parvals['tau_1']
    tau_2 = parvals['tau_2']
    rho_1 = parvals['rho_1']

    theta_1 = parvals['theta_' + str(season)]

    grav_model = theta_1 * (np.power(dist, rho_1)) / ((np.power(host_pop, tau_1)) * (np.power(targ_pop, tau_2)))

    return grav_model - data

按照 M Newville 的建议，我得到了合适的工作。下面的代码：

import pandas as pd
from lmfit import minimize, Parameters, report_fit
import numpy as np

variable_data = pd.read_csv("../Data/variables_for_model_fit.csv", sep=",", header='infer')

season_list = sorted(variable_data.flu_season.unique())


#### create parameters ####
params1 = Parameters()
params1.add('tau_host', value=0.24,min=0, max = 3)
params1.add('tau_targ', value=0.14,min=0, max = 3)
params1.add('rho', value=0.29,min=0, max = 3)
# "global" parameters

for i, season in enumerate(season_list):
    params1.add('theta_%d' % season_list[i], value=1000., min=1, max=1e6)
    # creates theta parameters for each season

#### define model ####
def grav_dist_over_pop(dist, hist_pop, targ_pop, theta, rho, tau_host, tau_targ):
    return theta**(-1) * dist**rho * host_pop**(-tau_host) * targ_pop**(-tau_targ)

#### objective function ####
def objective_1(params, dist, host_pop, targ_pop, flu_sea, data):
    parvals1 = params1.valuesdict()
    resid = np.zeros((len(season_data),len(season_data)))

    for i, data in enumerate(data):

        theta = parvals1['theta_%d' % flu_sea[i]]
        rho = parvals1['rho']#_%d' % flu_sea[i]]
        tau_host = parvals1['tau_host']#_%d' % flu_sea[i]]
        tau_targ = parvals1['tau_targ']#_%d' % flu_sea[i]]
        model = grav_dist_over_pop(dist, host_pop, targ_pop, theta, rho, tau_host, tau_targ)
        resid[i, :] = model - data
    return resid.flatten()

#### fit global variables ####
season_data = variable_data.sample(500)
# my dataset is huge so lmfit takes an age when fitting all of the
# theta values
host_pop = np.asarray(season_data.host_city_pop.values.tolist())
targ_pop = np.asarray(season_data.target_city_pop.values.tolist())
dist = np.asarray(season_data.distance.values.tolist())
data = np.asarray(season_data.time_to_spread.values.tolist())
flu_sea = np.asarray(season_data.flu_season.values.tolist())

result = minimize(objective_1, params1, args=(dist, host_pop, targ_pop, flu_sea, data))

report_fit(result.params)
4

1 回答 1

0

与lmfit（或scipy.optimize我知道的所有其他工具）的拟合始终是“全局”拟合，为一组用于最小化一维残差数组的参数找到最佳值。需要明确的是，您可以优化参数以适应多个数据集（或组或季节），但您必须将问题简化为计算一维残差数组的单个参数集。

对于您的问题，我建议从定义为单个季节或数据集建模数据的“重力模型”开始。我对这类领域一无所知（但很高兴 lmfit 可能有用！），但从您的示例来看，这可能看起来像这样（如果不是，请更正）

def grav_season(dist, host_pop, targ_pop, theta, rho, tau_host, tau_targ):
    return theta * dist**rho * host_pop**(-tau_host) * targ_pop**(-tau_targ)

在我看来，它有 3 个自变量（dist, host_pop, targ_pop）和 4 个可能是变量的参数：theta, rho, tau_host, tau_targ. 同样，如果需要，请更正，但出于此处的目的，这些细节并不那么重要。

要适合单个季节/数据集，您可以这样做

from lmfit import Parameters, minimize, report_fit

def objective(params, data, dist, host_pop, targ_pop):
     pars = params.valuesdict()
     model = grav_season(dist, host_pop, targ_pop, pars['theta'],
                         pars['rho'], pars['tau_host'], pars['tau_targ'])
     return (model - data)

params = Parameters()
param.add('theta', value=1., min = 0.01)
param.add('rho', value=1., min = 0.01)
param.add('tau_host', value=1., min = 0.01)
param.add('tau_targ, value=1., min = 0.01)

result = minimize(objective, params, args=(data, dist, host_pop, targ_pop))
report_fit(result.params)

现在，对于多个季节，您可以简单地为每个季节制作一个theta, rho, tau_host,参数。tau_targ但是，将所有这些数据一起使用并没有多大意义。如果我理解正确，您希望tau_host,tau_targ和rho在所有季节都具有相同的值。

为此，请为全局应用的指数创建参数：

params = Parameters()
param.add('rho', value=1., min = 0.01)
param.add('tau_host', value=1., min = 0.01)
param.add('tau_targ', value=1., min = 0.01)

和每个季节的参数：

for i, season in enumerate(sorted(variable_data.flu_season.unique())):
     param.add('theta_%d' % i, value=1., min = 0.01)
     param.add('rho_%d' % i,  expr='rho')
     param.add('tau_host_%d' % i, expr='tau_host')
     param.add('tau_targ_%d' % i, expr='tau_targ')

请注意，theta_i将独立变化，但是rho_i, etc 将被限制为变量rho等的值。这会为每个季节提供一整套参数，但符合您的约束条件的参数。如果看起来可能需要更多测试，这种方法允许您轻松更改其中的一个或多个以分别改变。为此，您可以说：

    params['rho_7'].set(value=2.0, vary=True, expr='')  # vary independently

要使用这组多季节参数，您还需要多季节数据。我不确定和是否会随季节变化，或者dist仅。我会假设只会随着季节而变化（但如果不是，它应该很容易改变）。为每个季节构建一个包含数据的列表，然后修改您的目标函数，如下所示：host_poptarg_popdatadata

def objective(params, season_data, dist, host_pop, targ_pop):
    pars = params.valuesdict()
    resid = np.zeros((len(season_data), len(season_data[0]))
    for i, data in enumerate(season_data):
        theta = pars['theta_%d' % i]
        rho = pars['rho_%d' % i]
        tau_host = pars['tau_host_%d' % i]
        tau_targ = pars['tau_targ_%d' % i]
        model = grav_season(dist, host_pop, targ_pop,
                            theta, rho, tau_host, tau_targ)
        resid[i, :] = model - data
    return resid.flatten()

result = minimize(objective, params, args=(season_data, dist, host_pop, targ_pop))
report_fit(result.params)

希望这有助于您入门。同样，主要建议是：

  1. 为您正在建模的现象创建一个每个季节/每个数据集的模型功能。
  2. 使用 'expr' 来约束每个季节的参数以采用“全局值”。尽管您希望所有季节都具有相同的值，但可能某些参数可能会随季节线性变化，或者需要添加到某个值或其他意味着它们并不是真正独立的值。
于 2018-02-14T04:05:29.787 回答