1. 问题描述

在回归中，我们往往会碰到有约束的系数问题。
例如，在回归方程中No.006 Python 最优化问题 optimize - 图1中，人为规定：①回归系数全部在[0,1]之间，②系数之和为1。那么为了求解OLS方程，有约束的回归模型求解数学表达式为
No.006 Python 最优化问题 optimize - 图2

那么该问题需要用到scipy.optimize解决。

2. 数据

data.csv

3. 代码准备

  1. import numpy as np
  2. import pandas as pd
  3. from scipy.optimize import minimize
  5. def get_coefficients(init_coeff, X, y, method):
  7.     # 计算sum squared errors
  8.     funcSSE = lambda k, X, y: sum((y - np.dot(X, np.array(k).reshape((-1, 1)))) ** 2)[0]
  10.     # 等式约束条件：eq
  11.     cons = ({'type':'eq', 'fun': lambda k: sum(k) - 1})
  12.     # 如果有不等式约束条件，那么是ineq，且是≥。例如，k0>0 & k0<1。
  13.     # cons = ({'type':'ineq', 'fun': lambda k: 1 - k[0]},
  14.     #         {'type':'ineq', 'fun': lambda k: k[0] - 0})
  16.     # 参数的上下限
  17.     bnds = ((0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1))
  19.     # optimize函数：
  20.     # 输入函数fun，估计参数的初始值init_coeff，估计参数之外的参数args = (X,y)，
  21.     # 求解方法method，上下限bounds，精确度tol，约束条件constraints。
  22.     result = minimize(fun=funcSSE, 
  23.                       x0=init_coeff, 
  24.                       args=(X, y), 
  25.                       method=method, 
  26.                       bounds=bnds, 
  27.                       tol=0.001, 
  28.                       constraints=cons)
  29.     # 返回估计参数值
  30.     return result

这里尝试了全部的求解方法，不仅要得到全局最优解，还要满足约束条件。

  1. # 读取数据
  2. data = pd.read_csv(r'C:\Users\Administrator\Downloads\data.csv').T
  3. data.columns = data.iloc[0]
  4. data = data.iloc[2:]
  6. # 数据格式指定为float，否则后续回归会有问题。
  7. for var in list(data.columns):
  8.     data[var] = data[var].apply(lambda x: float(x))
  10. # 设定X和y，以及初始参数
  11. X = np.array(data.iloc[:, :6]).reshape((-1, 6))
  12. y1 = np.array(data.iloc[:, 7]).reshape((-1, 1))
  13. init_coeff = [0, 0, 0, 0, 0.2, 0.8]
  15. # 设定SSE计算公式。
  16. funcSSE = lambda k, X, y: sum((y - np.dot(X, np.array(k).reshape((-1, 1)))) ** 2)[0]
  18. for m in ['Nelder-Mead', 'Powell', 'CG', 'BFGS', 'Newton-CG', 'L-BFGS-B', 'TNC', 'COBYLA', 
  19.           'SLSQP', 'trust-constr', 'dogleg', 'trust-ncg', 'trust-exact', 'trust-krylov']:
  20.     try:        
  21.         result = get_coefficients(init_coeff, X, y1, m)
  22.         k_hat = result.x
  23.         print('method = %s.'%m)
  24.         print('sum of coefficients is %0.4f.'%sum(k_hat))
  25.         print('estimated coeff is ', k_hat)
  26.         print('SSE of estimated coeff is %0.4f'%funcSSE(k_hat, X, y1))
  27.         print('\n')
  28.     except:
  29.         continue

结果如下

  1. method = Nelder-Mead.
  2. sum of coefficients is 0.0082.
  3. estimated coeff is  [-0.00275149  0.00436154  0.00239604  0.00166577  0.01761348 -0.01503714]
  4. SSE of estimated coeff is 0.9582
  7. method = Powell.
  8. sum of coefficients is 0.0101.
  9. estimated coeff is  [5.81772244e-03 2.11612816e-05 6.16658126e-05 6.61069614e-05
  10.  4.05846939e-03 6.61069614e-05]
  11. SSE of estimated coeff is 0.9780
  14. method = CG.
  15. sum of coefficients is 0.0091.
  16. estimated coeff is  [ 0.00839106  0.00693108 -0.0215021   0.0058082   0.03196229 -0.02249176]
  17. SSE of estimated coeff is 0.9322
  20. method = BFGS.
  21. sum of coefficients is 0.0091.
  22. estimated coeff is  [ 0.00839105  0.00693109 -0.02150201  0.00580812  0.03196221 -0.0224917 ]
  23. SSE of estimated coeff is 0.9322
  26. method = L-BFGS-B.
  27. sum of coefficients is 0.0101.
  28. estimated coeff is  [0.00591175 0.         0.         0.         0.00415395 0.        ]
  29. SSE of estimated coeff is 0.9775
  32. method = TNC.
  33. sum of coefficients is 0.0101.
  34. estimated coeff is  [0.00591182 0.         0.         0.         0.00415385 0.        ]
  35. SSE of estimated coeff is 0.9775
  38. method = SLSQP.
  39. sum of coefficients is 1.0000.
  40. estimated coeff is  [1.68700952e-14 3.56512119e-01 4.81162901e-01 1.88647109e-15
  41.  1.62324980e-01 1.21545881e-14]
  42. SSE of estimated coeff is 160.4392
  45. method = trust-constr.
  46. sum of coefficients is 1.0000.
  47. estimated coeff is  [3.48435748e-06 3.56479768e-01 4.81183407e-01 3.62541898e-05
  48.  1.62291535e-01 5.55114801e-06]
  49. SSE of estimated coeff is 160.4395

根据结果显示，只有SLSQP和trust-constr得到的解能够满足约束条件，而其他方案得到的解根本不满足约束。尽管SLSQP和trust-constr计算得到SSE很大，但是为了满足约束，因此推荐使用SLSQP和trust-constr进行规划求解。

4. 思考

那些指定了约束条件、但是依然没有满足的求解方法，得到的结果是怎么来的？这里，我们尝试了OLS直接回归。

  1. import statsmodels.api as sm
  3. linear_reg = sm.OLS(np.array(data.iloc[:, 6:7]), np.array(data.iloc[:, :6])).fit()
  5. linear_reg.summary()

结果如下。观察到，OLS直接回归的系数结果和BFGS/CG一致。因此，BFGS/CG求解得到的是全局最优解，且置于约束条件不顾。
image.png
尽管SLSQP/trust-constr并非最优解，但是至少满足了约束条件。