如何约束耦合单变量线性回归的系数？含m_i求和为0约束

你需要拟合满足m_i * x + c_i = y_i且sum(m_i) = 0的线性模型，以下是基于scipy和statsmodels的两种解决方案，均支持量化抽样变异性与置信区间：

方法一：使用scipy.optimize实现带约束最小二乘拟合 #

scipy.optimize.least_squares支持添加线性约束，适合这类问题。由于scipy不直接输出置信区间，我们用bootstrap自助法估计抽样变异性。

代码示例 #

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

# 生成模拟数据（替换为你的观测时间序列）
np.random.seed(42)
n_groups = 3  # 组数对应m_i的数量
x = np.linspace(0, 10, 50)
true_m = np.array([1, -0.5, -0.5])  # 满足sum(m_i)=0
true_c = np.array([2, 5, 1])
y = []
for m, c in zip(true_m, true_c):
    y.append(m * x + c + np.random.normal(0, 0.5, size=x.shape))
y = np.concatenate(y)
x_full = np.tile(x, n_groups)
group_idx = np.repeat(np.arange(n_groups), len(x))  # 标记数据所属组

# 定义残差函数：参数格式为[m0, m1, m2, c0, c1, c2]
def residuals(params, x, y, group_idx, n_groups):
    m = params[:n_groups]
    c = params[n_groups:]
    y_pred = np.zeros_like(y)
    for i in range(n_groups):
        mask = group_idx == i
        y_pred[mask] = m[i] * x[mask] + c[i]
    return y_pred - y

# 定义斜率和为0的约束
cons = {'type': 'eq',
        'fun': lambda params: np.sum(params[:n_groups])}

# 初始参数猜测
init_params = np.concatenate([np.zeros(n_groups), np.mean(y)*np.ones(n_groups)])

# 拟合模型
result = least_squares(residuals, init_params, args=(x_full, y, group_idx, n_groups), constraints=cons)
estimated_m = result.x[:n_groups]
estimated_c = result.x[n_groups:]

print("估计斜率m_i:", estimated_m)
print("估计截距c_i:", estimated_c)
print("斜率和:", np.sum(estimated_m))

# Bootstrap估计95%置信区间
n_bootstrap = 1000
bootstrap_params = []
for _ in range(n_bootstrap):
    idx = np.random.choice(len(y), len(y), replace=True)
    boot_result = least_squares(residuals, init_params, args=(x_full[idx], y[idx], group_idx[idx], n_groups), constraints=cons)
    bootstrap_params.append(boot_result.x)

bootstrap_params = np.array(bootstrap_params)
ci_m = np.percentile(bootstrap_params[:, :n_groups], [2.5, 97.5], axis=0)
ci_c = np.percentile(bootstrap_params[:, n_groups:], [2.5, 97.5], axis=0)

print("\n斜率m_i的95%置信区间:")
for i in range(n_groups):
    print(f"m{i}: [{ci_m[0,i]:.4f}, {ci_m[1,i]:.4f}]")

print("\n截距c_i的95%置信区间:")
for i in range(n_groups):
    print(f"c{i}: [{ci_c[0,i]:.4f}, {ci_c[1,i]:.4f}]")

# 可视化拟合结果
plt.figure(figsize=(10,6))
for i in range(n_groups):
    mask = group_idx == i
    plt.scatter(x_full[mask], y[mask], label=f"组{i}原始数据")
    plt.plot(x, estimated_m[i]*x + estimated_c[i], label=f"组{i}拟合线", linewidth=2)
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("带斜率和为0约束的线性拟合结果")
plt.show()

方法二：使用statsmodels实现带约束OLS拟合 #

statsmodels的OLS支持通过fit_constrained直接添加线性约束，且内置统计推断功能，可直接输出置信区间。

代码示例 #

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

# 生成模拟数据（替换为你的观测时间序列）
np.random.seed(42)
n_groups = 3
x = np.linspace(0, 10, 50)
true_m = np.array([1, -0.5, -0.5])
true_c = np.array([2, 5, 1])
y = []
for m, c in zip(true_m, true_c):
    y.append(m * x + c + np.random.normal(0, 0.5, size=x.shape))
y = np.concatenate(y)

# 构建设计矩阵：每组对应一个x项和一个截距项
X = np.zeros((len(y), 2*n_groups))
for i in range(n_groups):
    start = i*len(x)
    end = start + len(x)
    X[start:end, i] = x  # 对应斜率m_i
    X[start:end, n_groups+i] = 1  # 对应截距c_i

# 定义约束：m0 + m1 + m2 = 0
constraint = 'x0 + x1 + x2 = 0'

# 拟合带约束的OLS模型
model = sm.OLS(y, X)
result = model.fit_constrained(constraint)

# 输出详细统计结果
print(result.summary())

# 提取参数与置信区间
estimated_m = result.params[:n_groups]
estimated_c = result.params[n_groups:]
ci_m = result.conf_int()[:n_groups]
ci_c = result.conf_int()[n_groups:]

print("\n估计斜率m_i:", estimated_m)
print("估计截距c_i:", estimated_c)
print("斜率和:", np.sum(estimated_m))

print("\n斜率m_i的95%置信区间:")
for i in range(n_groups):
    print(f"m{i}: [{ci_m.iloc[i,0]:.4f}, {ci_m.iloc[i,1]:.4f}]")

print("\n截距c_i的95%置信区间:")
for i in range(n_groups):
    print(f"c{i}: [{ci_c.iloc[i,0]:.4f}, {ci_c.iloc[i,1]:.4f}]")

# 可视化拟合结果
plt.figure(figsize=(10,6))
for i in range(n_groups):
    start = i*len(x)
    end = start + len(x)
    plt.scatter(x, y[start:end], label=f"组{i}原始数据")
    plt.plot(x, estimated_m[i]*x + estimated_c[i], label=f"组{i}拟合线", linewidth=2)
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("statsmodels带约束OLS拟合结果")
plt.show()

方法对比 #

• scipy方案：灵活性高，适合复杂约束场景，但需手动实现bootstrap获取置信区间；
• statsmodels方案：内置统计推断，直接输出置信区间与统计量，代码更简洁，适合需要严谨统计分析的场景。

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