如何拟合具有共享和非共享参数组合的多个数据集
How to fit multiple datasets which have a combination of shared and non-shared parameters
我正在尝试拟合多个数据集,这些数据集应该有一些变量,这些变量在数据集之间共享,而其他变量则没有。但是我不确定我需要采取的步骤来做到这一点。下面我展示了我尝试使用的方法(来自 'Issues begin here' 不起作用,它仅用于说明目的)。
In this answer 有人能够共享跨数据集的参数是否有某种方法可以对此进行调整,以便我也可以拥有一些非共享参数?
有没有人知道我如何实现这一目标,或者有人可以建议一种更好的方法来实现相同的结果?谢谢
import numpy as np
from scipy.stats import gamma
import matplotlib.pyplot as plt
import pandas as pd
from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit, Model
# Create datasets to fit
a = 1.99
start = gamma.ppf(0.001, a)
stop = gamma.ppf(.99, a)
xvals = np.linspace(start, stop, 100)
yvals = gamma.pdf(xvals, a)
data_dict = {}
for dataset in range(4):
name = 'dataset_' + str(dataset)
rand_offset = np.random.uniform(-.1, .1)
noise = np.random.uniform(-.05, .05,len(yvals)) + rand_offset
data_dict[name] = yvals + noise
df = pd.DataFrame(data_dict)
# Create some non-shared parameters
non_shared_param1 = np.random.uniform(0.009, .21, 4)
non_shared_param2 = np.random.uniform(0.01, .51, 4)
# Create the independent variable
ind_var = np.linspace(.001,100,100)
# Create a model
def model_func(time, Qi, at, vw, R, rhob_cb, al, NSP1, NSP2):
Dt = at * vw
Dl = al * vw
t = time
first_bot = 8 * np.pi * t * rhob_cb
sec_bot = np.sqrt((np.pi * (Dl * R) * t))
exp_top = R * np.power((NSP1 - ((t * vw)/R)), 2)
exp_bot = 4 * Dl * t
exp_top2 = R * np.power(NSP2, 2)
exp_bot2 = 4 * Dt * t
return (Qi / first_bot * sec_bot) * np.exp(- (exp_top / exp_bot) - (exp_top2 / exp_bot2))
model = Model(model_func)
### Issues begin here ###
all_results = {}
index = 0
for col in df:
# This block assigns the correct non-shared parameter for the particular fit
nsp1 = non_shared_param1[index]
nsp2 = non_shared_param2[index]
index += 1
params = Parameters()
at = 0.1
al = 0.15
vw = 10**-4
Dt = at * vw
Dl = al * vw
# Non-shared parameters
model.set_param_hint('NSP1', value = nsp1)
model.set_param_hint('NSP2', value = nsp2)
# Shared and varying parameters
model.set_param_hint('vw', value =10**-4, min=10**-10)
model.set_param_hint('at', value =0.1)
model.set_param_hint('al', value =0.15)
# Shared and fixed parameters
model.set_param_hint('Qi', value = 1000, vary = True)
model.set_param_hint('R', value = 1.7, vary = True)
model.set_param_hint('rhob_cb', value =2895, vary = True)
# One set of parameters should be returned
result = model.fit(df[col], time = ind_var)
all_results[index] = result
与 lmfit 的拟合总是使用 Parameters 对象的单个实例,它不采用多个 Parameters 对象。
为了同时用相似的模型拟合多个数据集(也许是相同的数学模型,但期望每个模型的参数值不同),您需要一个 objective 函数来连接来自不同的组件模型。这些模型中的每一个都必须具有从 Parameters() 的单个实例中获取的参数,每个参数都有一个唯一的名称。
因此,为了拟合具有相同函数的 2 个数据集(让我们使用高斯,参数为 "center"、"amplitude" 和 "sigma"),您可以将参数定义为
params = Parameters()
params.add('center_1', 5., vary=True)
params.add('amplitude_1', 10., vary=True)
params.add('sigma_1', 1.0, vary=True
params.add('center_2', 8., vary=True)
params.add('amplitude_2', 3., vary=True)
params.add('sigma_2', 2.0, vary=True)
然后使用'center_1'、'amplitude_1'和'sigma_1'计算第一个数据集的模型,'center_2'等计算第二个数据集的模型, 也许是
def residual(params, x, datasets):
model1 = params['amplitude_1'] * gaussian(x, params['center_1'], params['sigma_1'])
model2 = params['amplitude_2'] * gaussian(x, params['center_2'], params['sigma_2']
resid1 = datasets[0] - model1
resid2 = datasets[1] - model2
return np.concatenate((resid1, resid2))
fit = lmfit.minimize(residual, params, fcn_args=(x, datasets))
从这里您可能会看到,默认情况下参数值是独立的。为了共享要在不同数据集中使用的参数值,您必须明确地执行此操作(如您提供的链接答案所示)。
比如你要要求sigma
的值相同,不改变残差函数,只需将上面的Parameter定义为:
params.add('sigma_2', expr='sigma_1')
您可能需要将两个振幅相加到某个值:
params.add('amplitude_2', expr='10 - amplitude_1')
或者您可能希望确保 'center_2' 大于 'center_1',但要确定合适的数量:
params.add('center_offset', value=0.5, min=0)
params.add('center_2', expr='center_1 + center_offset')
这些都是绑定参数值的方法。默认情况下,它们是独立的。当然,你也可以有一些在所有模型中使用的参数(比如,只需调用参数 'sigma' 并在所有模型中使用它)。
我正在尝试拟合多个数据集,这些数据集应该有一些变量,这些变量在数据集之间共享,而其他变量则没有。但是我不确定我需要采取的步骤来做到这一点。下面我展示了我尝试使用的方法(来自 'Issues begin here' 不起作用,它仅用于说明目的)。
In this answer 有人能够共享跨数据集的参数是否有某种方法可以对此进行调整,以便我也可以拥有一些非共享参数?
有没有人知道我如何实现这一目标,或者有人可以建议一种更好的方法来实现相同的结果?谢谢
import numpy as np
from scipy.stats import gamma
import matplotlib.pyplot as plt
import pandas as pd
from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit, Model
# Create datasets to fit
a = 1.99
start = gamma.ppf(0.001, a)
stop = gamma.ppf(.99, a)
xvals = np.linspace(start, stop, 100)
yvals = gamma.pdf(xvals, a)
data_dict = {}
for dataset in range(4):
name = 'dataset_' + str(dataset)
rand_offset = np.random.uniform(-.1, .1)
noise = np.random.uniform(-.05, .05,len(yvals)) + rand_offset
data_dict[name] = yvals + noise
df = pd.DataFrame(data_dict)
# Create some non-shared parameters
non_shared_param1 = np.random.uniform(0.009, .21, 4)
non_shared_param2 = np.random.uniform(0.01, .51, 4)
# Create the independent variable
ind_var = np.linspace(.001,100,100)
# Create a model
def model_func(time, Qi, at, vw, R, rhob_cb, al, NSP1, NSP2):
Dt = at * vw
Dl = al * vw
t = time
first_bot = 8 * np.pi * t * rhob_cb
sec_bot = np.sqrt((np.pi * (Dl * R) * t))
exp_top = R * np.power((NSP1 - ((t * vw)/R)), 2)
exp_bot = 4 * Dl * t
exp_top2 = R * np.power(NSP2, 2)
exp_bot2 = 4 * Dt * t
return (Qi / first_bot * sec_bot) * np.exp(- (exp_top / exp_bot) - (exp_top2 / exp_bot2))
model = Model(model_func)
### Issues begin here ###
all_results = {}
index = 0
for col in df:
# This block assigns the correct non-shared parameter for the particular fit
nsp1 = non_shared_param1[index]
nsp2 = non_shared_param2[index]
index += 1
params = Parameters()
at = 0.1
al = 0.15
vw = 10**-4
Dt = at * vw
Dl = al * vw
# Non-shared parameters
model.set_param_hint('NSP1', value = nsp1)
model.set_param_hint('NSP2', value = nsp2)
# Shared and varying parameters
model.set_param_hint('vw', value =10**-4, min=10**-10)
model.set_param_hint('at', value =0.1)
model.set_param_hint('al', value =0.15)
# Shared and fixed parameters
model.set_param_hint('Qi', value = 1000, vary = True)
model.set_param_hint('R', value = 1.7, vary = True)
model.set_param_hint('rhob_cb', value =2895, vary = True)
# One set of parameters should be returned
result = model.fit(df[col], time = ind_var)
all_results[index] = result
与 lmfit 的拟合总是使用 Parameters 对象的单个实例,它不采用多个 Parameters 对象。
为了同时用相似的模型拟合多个数据集(也许是相同的数学模型,但期望每个模型的参数值不同),您需要一个 objective 函数来连接来自不同的组件模型。这些模型中的每一个都必须具有从 Parameters() 的单个实例中获取的参数,每个参数都有一个唯一的名称。
因此,为了拟合具有相同函数的 2 个数据集(让我们使用高斯,参数为 "center"、"amplitude" 和 "sigma"),您可以将参数定义为
params = Parameters()
params.add('center_1', 5., vary=True)
params.add('amplitude_1', 10., vary=True)
params.add('sigma_1', 1.0, vary=True
params.add('center_2', 8., vary=True)
params.add('amplitude_2', 3., vary=True)
params.add('sigma_2', 2.0, vary=True)
然后使用'center_1'、'amplitude_1'和'sigma_1'计算第一个数据集的模型,'center_2'等计算第二个数据集的模型, 也许是
def residual(params, x, datasets):
model1 = params['amplitude_1'] * gaussian(x, params['center_1'], params['sigma_1'])
model2 = params['amplitude_2'] * gaussian(x, params['center_2'], params['sigma_2']
resid1 = datasets[0] - model1
resid2 = datasets[1] - model2
return np.concatenate((resid1, resid2))
fit = lmfit.minimize(residual, params, fcn_args=(x, datasets))
从这里您可能会看到,默认情况下参数值是独立的。为了共享要在不同数据集中使用的参数值,您必须明确地执行此操作(如您提供的链接答案所示)。
比如你要要求sigma
的值相同,不改变残差函数,只需将上面的Parameter定义为:
params.add('sigma_2', expr='sigma_1')
您可能需要将两个振幅相加到某个值:
params.add('amplitude_2', expr='10 - amplitude_1')
或者您可能希望确保 'center_2' 大于 'center_1',但要确定合适的数量:
params.add('center_offset', value=0.5, min=0)
params.add('center_2', expr='center_1 + center_offset')
这些都是绑定参数值的方法。默认情况下,它们是独立的。当然,你也可以有一些在所有模型中使用的参数(比如,只需调用参数 'sigma' 并在所有模型中使用它)。