文章目录
- 问题描述
- 等价的矩阵形式
- 产生数据
- 使用MindOpt求解器的API
- 方法1：cell中直接输入代码运行
- 方法2：命令行直接运行.py文件
- 求解结果
- 联系我们
回归分析是一种预测技术，目的是建立 自变量x（向量）和 相关变量y（标量）之间的关系。目前有七种常见的回归分析：Linear Regression线性回归（本篇）、Logistic Regression逻辑回归、Polynomial Regression多项式回归、Stepwise Regression逐步回归、Ridge Regression岭回归、Lasso Regression套索回归、ElasticNet回归。
本篇我们讲述的是**Linear Regression线性回归中的鲁棒线性回归。**鲁棒回归又称为稳健回归，是统计学稳健估计的方法之一，主要思路是对离群点十分敏感的最小二乘回归中的的函数进行修改。
下文我们将讲述一个例子，来展示MindOpt优化鲁棒线性回归问题（含源代码）。
问题描述
学生的行为习惯和测试成绩之间的关系。线性回归使用线性函数 $y = <\boldsymbol a,\boldsymbol x> + b $ 来描述这种关系（更准确地说法是仿射函数）。向量 表示“梯度”，标量
表示“截距”。通过
次观测
，我们可以估计出
和
普通的最小二乘回归（least-squares regression）对“离群观测”很敏感。仅一个离群观测就会影响对 和
假设“离群观测”的数量相对 而言很少，其他正常观测又很准确。则我们应该使用鲁棒线性回归去恢复真实的
和
相比最小二乘回归，鲁棒回归可以更好地抵御（但不也能完全消除）离群观测的影响。
通过变换和引入变量
为了编程方便，我们将变量放在不等式一侧，非变量参数放在另一侧：
等价的矩阵形式
将 x_1, …, x_m 作为行放入矩阵 X，将 y_1, …, y_m 作为元素放在列向量 y 中。令 1 为元素为1的向量。原始鲁棒线性回归问题可以简写为.
通过变换，我们可以把这个问题转化为线性规划：
产生数据
我们通过生成随机数据，验证鲁棒线性回归的有效性。现在随机产生真正的
再使用它们产生
外加
因为观测的次序对鲁棒线性回归没有任何影响，所以不失一般性，我们将离群观测放在最后。
我们在下面的代码里面会实现数据生成的过程。
使用MindOpt求解器的API
直接采用求解器的API，需要查阅API文档来理解API的意思，没有建模语言可读性高。请参阅https://solver.damo.alibaba.com/doc/html/API%20reference/API-python/index.html来查看PythonAPI的使用说明。
关于Python的例子，在文档的5.建模与优化求解章节有Python的示例。这里是个LP的问题，我们可以参考：https://solver.damo.alibaba.com/doc/html/model/lp/linear%20optimization-python.html
下面我们分两种方式描述在本平台环境中的运行方法：
方法1：cell中直接输入代码运行
请运行下面cell中的代码，点击本窗口上面的播放△运行，或者摁shift+enter键运行:
# LP_5_RobustRegression.py
from mindoptpy import *
import numpy as np
import random
def simulateData(dimension=4,m_observations=15,is_write_data = False,output_dir = ""):
# use time to seed the random number generator
random.seed(time.time())
# randomly generate vector a and scalar b
d = dimension
true_a = np.random.randn(d) * np.sqrt([d])
true_b = np.random.randn(1)[0]
# randomly generate observations (x[i],y[i]), i = 1,...,m with 10% outliers
m = m_observations
X = np.zeros((m,d))
y = np.zeros((m,1))
for i in range(m):
X[i,] = np.random.randn(d)
if i < np.floor(0.9*m):
y[i] = np.inner(true_a,X[i,]) + true_b
else:
y[i] = np.random.randn(1)
#print("true_a:\n",true_a)
#print("true_b:\n",true_b)
#print("X:\n",X)
#print("y:\n",y)
if is_write_data:
f_true_a = open(output_dir+"true_a.txt", "w")
L = ["{}\n".format(true_a[i]) for i in range(len(true_a))]
f_true_a.writelines(L)
f_true_b = open(output_dir+"true_b.txt", "w")
L = ["{}\n".format(true_b)]
f_true_b.writelines(L)
f_X = open(output_dir+"X.txt", "w")
L = ["{}\n".format(X[i]) for i in range(len(X))]
f_X.writelines(L)
f_y = open(output_dir+"y.txt", "w")
L = ["{}\n".format(y[i]) for i in range(len(y))]
f_y.writelines(L)
return true_a,true_b,X,y
if __name__ == "__main__":
# --- Create simulation data ---------------
d = 4
m = 15
true_a,true_b,X,y = simulateData(d,m)
#--- MindOpt Step 1. Initialize an empty optimization model ---
model = MdoModel()
try:
#--- MindOpt Step 2. Set model inputs ---
# Change to minimization problem.
model.set_int_attr("MinSense", 1)
# add variables, specify their low bounds, upper bounds, and cost coefficients
# set objective: minimize_{a,b,c} 0'a + 0b + 1'c
INF = MdoModel.get_infinity()
var_a = []
for j in range(d):
var_a.append(model.add_var(-INF, INF, 0.0)) # low bnd, upr bnd, 0 cost coeff
var_b = model.add_var(-INF, INF, 0) # low bnd, upr bnd, 0 cost coeff
var_c = []
for i in range(m):
var_c.append(model.add_var(-INF, INF, 1.0)) # low bnd, upr bnd, cost coeff = 1.0
# Add constraints row by row
for i in range(m):
# compute a'x[i,]
expr = MdoExprLinear()
for j in range(d):
expr += X[i,j]*var_a[j]
# add: y[i] <= a'x[i,] + b + c[i]
model.add_cons(y[i], INF, expr + var_b + var_c[i])
# add: a'x[i,] + b - c[i] <= y[i]
model.add_cons(-INF, y[i], expr + var_b - var_c[i])
#--- MindOpt Step 3. Set solver parameters ---
model.set_int_param("Method", 2) # as an example, set to run the interior-point method
model.set_int_param("NumThreads", 2) # as an example, set to run 2 threads
# ------Step3. Solve the problem and populate the result.-------
model.solve_prob()
model.display_results()
time.sleep(0.01) #for print
status_code, status_msg = model.get_status()
if status_msg == "OPTIMAL":
print("----\n")
print("The solver terminated with an OPTIMAL status (code {0}).".format(status_code))
print("目标函数总收益是: {0}".format(model.get_real_attr("PrimalObjVal")))
print("\n----------------\nCompare var_a to true_a, and var_b to true_b:\n----------------")
# Compare var_a to true_a, and var_b to true_b; Display residuals
# Display solutions a[0],...,a[d-1]
print('\n{0:<5} {1:>9} {2:>9}'.format('Entry','True','Soln'))
for j in range(d):
print('{0:>5} {1:>9f} {2:>9f}'.format('a['+'%s'%j+']', true_a[j],
var_a[j].get_real_attr('PrimalSoln')))
# Display solution b
print('\n{0:<5} {1:>9} {2:>9}'.format(' ','True','Soln'))
print('{0:>5} {1:>9f} {2:>9f}'.format('b', true_b, var_b.get_real_attr('PrimalSoln')))
print("\n----------------\nShow Residual(c value）:\n----------------")
# Display solutions c[0],...,c[m]
print('\n{0:>11} {1:>9}'.format('Observation','Residual'))
for i in range(m):
print('{0:>11d} {1:>9f}'.format(i, var_c[i].get_real_attr('PrimalSoln')))
else:
print("Optimizer terminated with a(n) {0} status (code {1}).".format(status_msg, status_code))
except MdoError as e:
print("Received Mindopt exception.")
print(" - Code : {}".format(e.code))
print(" - Reason : {}".format(e.message))
except Exception as e:
print("Received exception.")
print(" - Reason : {}".format(e))
finally:
# MindOpt Step 5. Free the model.
model.free_mdl()运行后得到如下结果：
Start license validation (current time : 01-MAR-2023 20:59:47).
License validation terminated. Time : 0.002s
Model summary.
- Num. variables : 20
- Num. constraints : 30
- Num. nonzeros : 180
- Bound range : [5.6e-01,9.6e+00]
- Objective range : [1.0e+00,1.0e+00]
- Matrix range : [3.5e-02,2.3e+00]
Presolver started.
Presolver terminated. Time : 0.000s
Interior point method started.
Iter PrimObj DualObj PrimFea DualFea GapFea Mu Time
0 +0.00000000e+00 -3.05311332e-15 2.2e+00 3.7e-02 3.6e+02 1.2e+01 0.01s
1 +3.87012911e+01 +1.10363647e+01 7.4e-02 1.0e-04 2.5e+00 1.2e+00 0.01s
2 +2.30985411e+01 +1.91163774e+01 2.0e-04 5.5e-06 2.1e-01 1.3e-01 0.01s
3 +1.98974978e+01 +1.98777823e+01 5.6e-07 1.5e-08 1.0e-03 6.6e-04 0.01s
4 +1.98802551e+01 +1.98802008e+01 1.5e-09 4.2e-11 2.7e-06 1.8e-06 0.01s
5 +1.98802076e+01 +1.98802075e+01 4.3e-12 1.2e-13 7.5e-09 5.0e-09 0.02s
Terminated.
- Method : Interior point method.
- Primal objective : 1.9880207641660E+01
- Dual objective : 1.9880207492544E+01
- Num. threads : 2
- Num. iterations : 5
- Solver details : Solver terminated with a primal/dual optimal status.
Iteration Objective Dual Inf. Primal Inf. Time
0 1.98802e+01 0.0000e+00 0.0000e+00 0.02s pushes: P(0) D(8)
8 1.98802e+01 0.0000e+00 0.0000e+00 0.02s
Interior point method terminated. Time : 0.027s
Optimizer summary.
- Optimizer used : Interior point method
- Optimizer status : OPTIMAL
- Total time : 0.028s
Solution summary. Primal solution
- Objective : 1.9880207511e+01
----
The solver terminated with an OPTIMAL status (code 1).
目标函数总收益是: 19.880207510885242
----------------
Compare var_a to true_a, and var_b to true_b:
----------------
Entry True Soln
a[0] 3.381614 3.381614
a[1] 2.077166 2.077166
a[2] 1.788548 1.788548
a[3] 1.662762 1.662762
True Soln
b -2.102494 -2.102494
----------------
Show Residual(c value）:
----------------
Observation Residual
0 0.000000
1 0.000000
2 0.000000
3 0.000000
4 0.000000
5 0.000000
6 0.000000
7 0.000000
8 0.000000
9 0.000000
10 0.000000
11 0.000000
12 0.000000
13 11.513372
14 8.366835
OpenBLAS WARNING - could not determine the L2 cache size on this system, assuming 256k方法2：命令行直接运行.py文件
上面是直接在cell中运行所有的脚本，我们也可以建立个新文档，将Python代码存在 src/python_src文件夹的LP_5_RobustRegression.py文件。然后在Launcher中打开Terminal，执行 python xx.py文件来运行。
您也可以下载本.py文件，在自己的电脑上安装MindOpt求解器，然后在自己电脑的环境运行。
Luancher可以点击左上角的+打方块打开，Terminal在最下方，如截图示意。打开的窗口可以拖动调整位置。
然后在Terminal命令行里运行如下指令:
cd src/python_src
python LP_5_RobustRegression.py运行得到的结果，特性类似算法1，由于每次数据随机生成的，所以数据不一样：
求解结果
本案例的数据是随机生成的，因此每次结果会不一样，但是从结果的True原值和鲁棒线性回归的结果Soln几乎一样，拟合度很好。
如下面是一个示例的结果打印。
Optimizer summary.
- Optimizer used : Interior point method
- Optimizer status : OPTIMAL
- Total time : 0.013s
Solution summary. Primal solution
- Objective : 9.7151651023e+00
----
The solver terminated with an OPTIMAL status (code 1).
目标函数总收益是: 9.715165102347509
----------------
Compare var_a to true_a, and var_b to true_b:
----------------
Entry True Soln
a[0] 2.710648 2.710648
a[1] 0.477721 0.477721
a[2] 1.164849 1.164849
a[3] 0.251547 0.251547
True Soln
b 0.878660 0.878660
----------------
Show Residual(c value）:
----------------
Observation Residual
0 0.000000
1 0.000000
2 0.000000
3 0.000000
4 0.000000
5 0.000000
6 0.000000
7 0.000000
8 0.000000
9 0.000000
10 0.000000
11 0.000000
12 0.000000
13 6.173857
14 3.541308
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
