用scipy.optimize.least_squares拟合多组Sigmoid及参数模型对比咨询

当然可以用scipy.optimize.least_squares来实现这两种模型的拟合和AIC对比！你遇到的参数扁平化混乱问题，核心是要把参数的逻辑结构和优化用的一维数组做清晰的映射，下面我给你梳理一套清晰易维护的实现思路：

一、先回顾你的模拟数据（保留原逻辑） #

首先把你生成模拟数据的代码保留，确保我们有一致的实验数据基础：

import numpy as np; import matplotlib.pyplot as plt;
# lifted, scaled, stretched and shifted sigmoid
def lssig ( t, base, height, stretch, t50 ):
    return base + height / ( 1 + np.exp ( -1/stretch * ( t - t50 ) ) );
# number of curves per experiment
numcurves = 10;
# number of experiments (time pts)
num_experiments = 20;
# timepoints (with simulated variability in timing)
time_jitter = 2;
tstart = np.linspace ( -10, 30, num = num_experiments );
tstart += time_jitter * np.random.rand ( num_experiments );
# measurements + added noise (this is added noise, not the variability in the estimates)
measurements = np.ndarray ( shape = ( num_experiments, numcurves ) );
measure_noise = np.random.normal ( 0, .2, ( num_experiments, numcurves ) );
# Variability of model parameters between curves:
# par_spreading [0] -> variability in base (start value)
# par_spreading [1] -> variability in height (end - start)
# par_spreading [2] -> variability in stretch (duration of step)
# par_spreading [3] -> variability in t50 (middle of step)
par_spreading = np.ones ( 4 );
# parameters of the model:
# every curve can have its own parameters, variability between curves
# is given by the par_spreading value (see above) for that parameter
height = 8 * np.ones ( numcurves ) + par_spreading [0] * np.random.rand ( numcurves);
base = 2 * np.ones ( numcurves ) + par_spreading [1] * np.random.rand ( numcurves);
stretch = 4 * np.ones ( numcurves ) + par_spreading [2] * np.random.rand ( numcurves);
t50 = 9 * np.ones ( numcurves ) + par_spreading [2] * np.random.rand ( numcurves);
# fill the measurement array
for t in range(num_experiments):
    for c in range (numcurves):
        measurements[ t, c ] = lssig ( tstart [t], base [c], height [c], stretch [c], t50 [c], );
measurements += measure_noise;

二、参数结构的清晰映射 #

为了避免扁平化参数的逻辑混乱，我们用打包/解包函数把逻辑结构和优化用的一维数组隔离开：

1. 全局共享参数模型 #

所有曲线共用一套[base, height, stretch, t50]参数，一维数组直接对应4个元素，逻辑非常直观。

2. 独立参数模型 #

每组曲线有自己的4个参数，一维数组长度为numcurves * 4，我们写一个解包函数来拆分参数：

def unpack_independent_params(theta, numcurves):
    # 将一维theta拆分为(numcurves, 4)的参数矩阵，每行对应一组曲线的4个参数
    return theta.reshape((numcurves, 4))

这样在目标函数里，我们可以清晰地调用对应组的参数，不用再纠结数组索引的含义。

三、针对两种模型的目标函数 #

我们分别为两种模型编写目标函数，逻辑清晰且易于维护：

1. 全局参数模型的目标函数 #

def fun_global(theta, tstart, measurements):
    # theta：[base, height, stretch, t50]，所有曲线共用这套参数
    residuals = []
    for curve_idx in range(numcurves):
        # 计算当前曲线的模型预测值
        model_vals = lssig(tstart, theta[0], theta[1], theta[2], theta[3])
        # 收集残差（模型值-测量值）
        residuals.extend(model_vals - measurements[:, curve_idx])
    return np.array(residuals)

2. 独立参数模型的目标函数 #

def fun_independent(theta, tstart, measurements, numcurves):
    # 先解包参数
    params = unpack_independent_params(theta, numcurves)
    residuals = []
    for curve_idx in range(numcurves):
        # 取出当前曲线的独立参数
        base_c, height_c, stretch_c, t50_c = params[curve_idx]
        # 计算当前曲线的模型预测值
        model_vals = lssig(tstart, base_c, height_c, stretch_c, t50_c)
        # 收集残差
        residuals.extend(model_vals - measurements[:, curve_idx])
    return np.array(residuals)

四、拟合模型并计算AIC对比 #

拟合完成后，我们用**赤池信息准则（AIC）**对比两种模型，AIC公式为：
AIC = 2*k + n*ln(RSS/n)
其中：

• k：模型的参数数量
• n：总数据点数量
• RSS：残差平方和

完整拟合与对比代码 #

from scipy.optimize import least_squares

# ----------------------
# 1. 拟合全局参数模型
# ----------------------
# 初始值用模拟数据的均值，加快收敛
theta0_global = [np.mean(base), np.mean(height), np.mean(stretch), np.mean(t50)]
res_global = least_squares(fun_global, theta0_global, args=(tstart, measurements))

# 计算全局模型的AIC
n_total = measurements.size  # 总数据点数量
k_global = 4  # 全局模型参数数
rss_global = np.sum(res_global.fun ** 2)
aic_global = 2 * k_global + n_total * np.log(rss_global / n_total)

# ----------------------
# 2. 拟合独立参数模型
# ----------------------
# 初始值：给每组曲线复制全局均值作为初始值
theta0_independent = np.tile(theta0_global, numcurves)
res_independent = least_squares(fun_independent, theta0_independent, args=(tstart, measurements, numcurves))

# 计算独立模型的AIC
k_independent = 4 * numcurves  # 独立模型参数数
rss_independent = np.sum(res_independent.fun ** 2)
aic_independent = 2 * k_independent + n_total * np.log(rss_independent / n_total)

# ----------------------
# 3. 模型对比结果
# ----------------------
print(f"全局共享参数模型AIC: {aic_global:.2f}")
print(f"独立参数模型AIC: {aic_independent:.2f}")
print(f"AIC差值（全局-独立）: {aic_global - aic_independent:.2f}")

注意：AIC越小的模型，在拟合度和复杂度之间的平衡越好，更适合你的数据

五、代码维护的小技巧 #

• 把参数解包、目标函数都封装成独立函数，后续修改模型逻辑（比如部分参数共享）时，只需要改动对应函数即可
• 初始值尽量贴近数据的统计特征（比如均值），可以大幅提升优化的收敛速度
• 如果需要可视化拟合结果，可以在拟合后调用lssig函数生成模型曲线，和原始测量值对比

内容的提问来源于stack exchange，提问作者alle_meije