多元线性模型的分位数回归
 
一、为什么要使用分位数回归？
二、分位数回归基本模型
三、分位数回归估计--线性规划
3.1损失函数
3.2目标函数
3.3线性规划
3.4回归算法
  
四、实际案例分析与python编程计算
4.1引入数据集
4.2计算$\hat{\beta }$
  
五、参考文献

 
一、为什么要使用分位数回归？
 
       众所周知，对于线性模型$Y=X\beta +\epsilon$，人们往往习惯于使用均值回归。但是均值回归往往更关注的只是均值，对于数据的“其他部位”往往照顾不够。有时不能较为客观地反映一组数据的各个层次的实际情况。比如我是一名还在读大三的本科生，我和当前的世界首富埃隆▪马斯克人均财富千亿美元级别，显然这个均值对于反映我的财富水平来说是毫无意义的，甚至是有信息误导的作用。但是如果试想使用分位数回归，那我就可以避免被马斯克平均了。
        我们知道在均值回归的普通最小二乘法中，我们是通过求残差的平方和最小来估计参数的。而在分位数回归中我们通常求残差的绝对值的加权求和最小来估计参数。就我目前了解，这样做的目的之一是减小离群点的大误差对整体回归估计的影响。例如，通常离群点处的残差是远远大于1的，在这种情况下残差的绝对值就要远小于残差的平方值，所以使用残差的绝对值可以减少离群点处的误差对整体回归的影响，而这也算是我们进行分位数回归的初衷之一。
        分位数回归通常是采用最小一乘法，而最小一乘法对误差项$\epsilon$并没有要求其服从正态分布，所以相比于均值回归通常采用的最小二乘法对误差$\epsilon$要求服从正态分布来说，最小一乘法的误差$\epsilon$更具有普遍性。
        总的来说相比于均值回归，（1）分位数回归对数据分布的情况掌握的更全面客观。（2）使用分位数回归，离群点对于数据整体的影响要比较使用均值回归小的多。所以我们也可以说分位数回归更加稳健。（3）分位数回归对于误差项更具有普适性。
 
二、分位数回归基本模型
 
$$Y^{(\theta )}=X{\beta }^{(\theta )}+{\epsilon }^{(\theta )}\text{\hspace{1em}(1)}$$
        其中$\theta$为分位数，对于实随机变量$Y$，其右连续分布函数为$F(y)=P(Y\le y)$,$Y$的$\theta$分位数函数为$F^{-1}(\theta )=inf(y:F(y)\ge \theta )$，也即第$100\theta \%$的$y$。
        其中，$$\begin{gathered} Y^{(\theta )}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right],X=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{1p}\\ 1& x_{21}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& \cdots & x_{np}\end{array}\right] \\ {\beta }^{(\theta )}=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right],{\epsilon }^{(\theta )}=\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_n\end{array}\right] \end{gathered}$$
 
三、分位数回归估计–线性规划
 
3.1损失函数
 
       定义损失函数为$${\rho }_{\theta }(\epsilon )=\epsilon (\theta -I(\epsilon ))\text{\hspace{1em}(2)}$$
        其中，
 $$I(\epsilon )=\{\begin{array}{l}0,\epsilon \ge 0\\ 1,\epsilon <0\end{array}\text{\hspace{1em}(3)}$$
 
3.2目标函数
 
       一般地，根据$Koenker(1978)$，对于$\theta$分位数回归，我们的目标函数通常采用如下函数：
 $$\underset{\beta \in R^m}{\min }\sum _{i=1}^n{\rho }_{\theta }(y_i-x_i^{\prime }\beta )\text{\hspace{1em}(4)}$$
 
3.3线性规划
 
       对于式$(4)$，我们可以用另一种等价形式表达：$$\underset{\beta \in R^m}{\min }\sum _{i=1}^n[\theta u_i+(1-\theta )v_i]\text{\hspace{1em}(5)}$$
        其中$u_i,v_i$分别是$(y_i-x_i^{\prime }\beta )$的正部和负部：
 $$u_i=\{\begin{array}{l}{y}_i-x_i^{\prime }\beta ,y_i-x_i^{\prime }\beta \ge 0\\ 0,y_i-x_i^{\prime }\beta <0\end{array}\text{\hspace{1em}(6)}$$
 $$v_i=\{\begin{array}{l}0,y_i-x_i^{\prime }\beta \ge 0\\ y_i-x_i^{\prime }\beta ,y_i-x_i^{\prime }\beta <0\end{array}\text{\hspace{1em}(7)}$$
        我们知道，一个函数等于它的正部$-$负部，一个函数的绝对值等于它的正部$+$负部。所以，我们可以得到下面的线性规划：
 $$\{\begin{array}{l}\min {\sum }_{i=1}^n[\theta u_i+(1-\theta )v_i]\\ x_i^{\prime }\beta +u_i-v_i=y_i\\ \beta \in R^m,u_i\ge 0,v_i\ge 0,i=1,2,\cdots ,n\end{array}\text{\hspace{1em}(8)}$$
        对于该线性规划，
 $$\{\begin{array}{l}A=\left[\begin{array}{ccc}X& I_n& -I_n\end{array}\right],B=Y\\ C=[0_{1\times (p+1)},{\theta }_{1\times n},(1-\theta {)}_{1\times n}]\end{array}\text{\hspace{1em}(9)}$$
 
3.4回归算法
 
$step(1)$：输入回归数据$X$和$Y$
 $step(2)$：根据$(9)$计算出线性规划的相关矩阵$A,B,C$
 $step(3)$：将$A,B,C$带入线性规划算法模块，求得最优解$(\beta ,u,v{)}^T$
 
注：关于相关的线性规划算法本文不再赘述。
 
四、实际案例分析与python编程计算
 
4.1引入数据集
 
       我们以著名的$Hald$数据为例。
 
 

 图4.1.1
 
4.2计算$\hat{\beta }$
 
       下面给出计算多元线性模型分位数回归方程完整$python$源代码：
 
import pandas as pd
import numpy as np
from scipy import optimize
#----------------------多元线性模型的分位数回归---------------------------------------------------------

#导入数据
dataset1=pd.read_excel('Hald.xlsx')
dataset2=pd.read_excel('Hald.xlsx')

#输入分位数theta
theta=[0.1,0.25,0.5,0.75,0.9]
#计算X,Y
Y=dataset1['Y'].values
dataset2['Y']=1
X=dataset2.values

#计算A,B
n=len(X)
B=Y
In=np.eye(n)
A=np.hstack((X,In,-In))
p=len(X[0])-1
l=np.ones(n,np.int)
o=np.zeros(p+1,np.int)

#给出线性规划的自变量取值范围
b=[]
for i in range(p+1):
    b.append(((None,None)))
for i in range(2*n):
    b.append((0,None))

#使用optimize包的linprog函数求解线性规划
for i in theta:
    C = np.hstack((o, i * l, (1 - i) * l))
    r = optimize.linprog(C,A_eq=A,b_eq=B,bounds=b)
    x=r.x
    #计算β估计
    beta=[]
    print('多元线性模型的{}分位数回归方程为：\ny='.format(i), end='')
    for i in range(p+1):
        beta.append(x[i])
    print(beta[0], end='')
    for i in range(1, p + 1):
        if beta[i] > 0:
            print('+{}x{}'.format(beta[i], i), end='')
        else:
            print('{}x{}'.format(beta[i], i), end='')
    print('\n',end='')
 
       下面给出程序运行结果：
 

 图4.2.1

 
五、参考文献
 
 
 
[1]分位数回归理论及其在金融风险测量中的应用/王新宇著.——北京：科学出版社,2010.6
 [2]吕书龙. 最小一乘估计快速算法的研究[D].福州大学,2003.
 

实验目的：了解广义多元线性模型的性质以及R语言的实现
 实验要求：对给出的两组数据建立线性回归模型。
 实验数据：分析习题3.1、习题3.2并给解释
 实验素材：在我的资源中
 
实验代码：
 1、
 
 
 
 2、
 
 
 实验结果分析：
 1、第一个数据分析得到x对应的P值比较大，即X不显著，采用回归分析，从回归分析中我们看到回归模型的系数在0.1%水平上显著。于是得回归模型为：预测出当年龄为30岁，该影片持肯定观点的人可能为49%。
 2、第二个数据分析得到x对应的P值比较小，即X显著，不用采用回归分析，预测出当年龄为30岁的女性，对该项服务产品持认可观点的人可能为79%。

 
文章目录
 
一、多元线性回归模型预测房价
1.基础包导入
2.变量探索
3.多元线性回归建模
4.模型优化
二、用Excel重做上面的多元线性回归，求解回归方程
三、用机器学习库Sklearn库重做上面的多元线性归
四、参考文献

 
一、多元线性回归模型预测房价
 
1.基础包导入
 
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.read_csv('house_prices.csv')
df.info(); df.head()
 
 
2.变量探索
 
# 异常值处理
# ================ 异常值检验函数：iqr & z分数 两种方法 =========================
def outlier_test(data, column, method=None, z=2):
	""" 以某列为依据，使用 上下截断点法 检测异常值(索引) """
	"""
	full_data: 完整数据
	column: full_data 中的指定行，格式 'x' 带引号
	return 可选; outlier: 异常值数据框
	upper: 上截断点; lower: 下截断点
	method：检验异常值的方法（可选, 默认的 None 为上下截断点法），
	选 Z 方法时，Z 默认为 2
	"""
	# ================== 上下截断点法检验异常值 ==============================
	if method == None:
		print(f'以 {column} 列为依据，使用 上下截断点法(iqr) 检测异常值...')
		print('=' * 70)
		# 四分位点；这里调用函数会存在异常
		column_iqr = np.quantile(data[column], 0.75) - np.quantile(data[column], 0.25)
		# 1，3 分位数
		(q1, q3) = np.quantile(data[column], 0.25), np.quantile(data[column], 0.75)
		# 计算上下截断点
		upper, lower = (q3 + 1.5 * column_iqr), (q1 - 1.5 * column_iqr)
		# 检测异常值
		outlier = data[(data[column] <= lower) | (data[column] >= upper)]
		print(f'第一分位数: {q1}, 第三分位数：{q3}, 四分位极差：{column_iqr}')
		print(f"上截断点：{upper}, 下截断点：{lower}")
		return outlier, upper, lower
	# ===================== Z 分数检验异常值 ==========================
	if method == 'z':
		""" 以某列为依据，传入数据与希望分段的 z 分数点，返回异常值索引与所在数据框 """
		"""
		params
		data: 完整数据
		column: 指定的检测列
		z: Z分位数, 默认为2，根据 z分数-正态曲线表，可知取左右两端的 2%，
		根据您 z 分数的正负设置。也可以任意更改，知道任意顶端百分比的数据集合
		"""
		print(f'以 {column} 列为依据，使用 Z 分数法，z 分位数取 {z} 来检测异常值...')
		print('=' * 70)
		# 计算两个 Z 分数的数值点
		mean, std = np.mean(data[column]), np.std(data[column])
		upper, lower = (mean + z * std), (mean - z * std)
		print(f"取 {z} 个 Z分数：大于 {upper} 或小于 {lower} 的即可被视为异常值。")
		print('=' * 70)
		# 检测异常值
		outlier = data[(data[column] <= lower) | (data[column] >= upper)]
		return outlier, upper, lower
 
outlier, upper, lower = outlier_test(data=df, column='price', method='z')
outlier.info(); outlier.sample(5)
 
 
# 这里简单的丢弃即可
df.drop(index=outlier.index, inplace=True)
 
# 类别变量，又称为名义变量，nominal variables
nominal_vars = ['neighborhood', 'style']
for each in nominal_vars:
	print(each, ':')
	print(df[each].agg(['value_counts']).T)
	# 直接 .value_counts().T 无法实现下面的效果
	## 必须得 agg，而且里面的中括号 [] 也不能少
	print('='*35)
	# 发现各类别的数量也都还可以，为下面的方差分析做准备
 
 
# 热力图
def heatmap(data, method='pearson', camp='RdYlGn', figsize=(10 ,8)):
    """
    data: 整份数据
    method：默认为 pearson 系数
    camp：默认为：RdYlGn-红黄蓝；YlGnBu-黄绿蓝；Blues/Greens 也是不错的选择
    figsize: 默认为 10，8
    """
    ## 消除斜对角颜色重复的色块
    # mask = np.zeros_like(df2.corr())
    # mask[np.tril_indices_from(mask)] = True
    plt.figure(figsize=figsize, dpi= 80)
    sns.heatmap(data.corr(method=method), \
        xticklabels=data.corr(method=method).columns, \
        yticklabels=data.corr(method=method).columns, cmap=camp, \
        center=0, annot=True)
    # 要想实现只是留下对角线一半的效果，括号内的参数可以加上 mask=mask
 
# 通过热力图可以看出 area，bedrooms，bathrooms 等变量与房屋价格 price 的关系都还比较强
## 所以值得放入模型，但分类变量 style 与 neighborhood 两者与 price 的关系未知
heatmap(data=df, figsize=(6,5))
 
 
# 刚才的探索我们发现，style 与 neighborhood 的类别都是三类，
 ## 如果只是两类的话我们可以进行卡方检验，所以这里我们使用方差分析
    
## 利用回归模型中的方差分析
## 只有 statsmodels 有方差分析库
## 从线性回归结果中提取方差分析结果
import statsmodels.api as sm
from statsmodels.formula.api import ols # ols 为建立线性回归模型的统计学库
from statsmodels.stats.anova import anova_lm
 
插播一条样本量和置信水平 α_level 的注意点（置信水平 α 的选择经验）
 
样本量 α-level
≤ 100 10%
100 ＜ n ≤ 500 5%
500 ＜ n ≤ 1000 1%
n ＞ 2000 千分之一
 
样本量过大，α-level 就没什么意义了。
 数据量很大时，p 值就没用了，样本量通常不超过 5000，
 为了证明两变量间的关系是稳定的，样本量要控制好。
 
# 数据集样本数量：6028，这里随机选择 600 条，如果希望分层抽样，可参考文章：
df = df.copy().sample(600)
# C 表示告诉 Python 这是分类变量，否则 Python 会当成连续变量使用
## 这里直接使用方差分析对所有分类变量进行检验
## 下面几行代码便是使用统计学库进行方差分析的标准姿势
lm = ols('price ~ C(neighborhood) + C(style)', data=df).fit()
anova_lm(lm)
# Residual 行表示模型不能解释的组内的，其他的是能解释的组间的
# df: 自由度（n-1）- 分类变量中的类别个数减1
# sum_sq: 总平方和（SSM），residual行的 sum_eq: SSE
# mean_sq: msm, residual行的 mean_sq: mse
# F：F 统计量，查看卡方分布表即可
# PR(>F): P 值
# 反复刷新几次，发现都很显著，所以这两个变量也挺值得放入模型中
 
 
3.多元线性回归建模
 
from statsmodels.formula.api import ols
lm = ols('price ~ area + bedrooms + bathrooms', data=df).fit()
lm.summary()
 

 
 
4.模型优化
 
发现精度还不够高，这里通过添加虚拟变量与使用方差膨胀因子检测多元共线性的方式来提升模型精度
 
# 设置虚拟变量
# 以名义变量 neighborhood 街区为例
nominal_data = df['neighborhood']
# 设置虚拟变量
dummies = pd.get_dummies(nominal_data)
dummies.sample() # pandas 会自动帮你命名
# 每个名义变量生成的虚拟变量中，需要各丢弃一个，这里以丢弃C为例
dummies.drop(columns=['C'], inplace=True)
dummies.sample()
 
 
# 将结果与原数据集拼接
results = pd.concat(objs=[df, dummies], axis='columns') # 按照列来合并
results.sample(3)
# 对名义变量 style 的处理可自行尝试
 
 
# 再次建模
lm = ols('price ~ area + bedrooms + bathrooms + A + B', data=results).fit()
lm.summary()
 
 
# 自定义方差膨胀因子的检测公式
def vif(df, col_i):
    """
    df: 整份数据
    col_i：被检测的列名
    """
    cols = list(df.columns)
    cols.remove(col_i)
    cols_noti = cols
    formula = col_i + '~' + '+'.join(cols_noti)
    r2 = ols(formula, df).fit().rsquared
    return 1. / (1. - r2)
 
test_data = results[['area', 'bedrooms', 'bathrooms', 'A', 'B']]
for i in test_data.columns:
    print(i, '\t', vif(df=test_data, col_i=i))
# 发现 bedrooms 和 bathrooms 存在强相关性，可能这两个变量是解释同一个问题
 
 
# 再次进行多元共线性检测
test_data = df[['area', 'bathrooms']]
for i in test_data.columns:
    print(i, '\t', vif(df=test_data, col_i=i))

 
 
二、用Excel重做上面的多元线性回归，求解回归方程
 
基于多元线性回归的房价预测
 摘要
 市场房价的走向受到多种因素的影响，通过对影响市场房价的多种因素进行分析，有助于对未来房价的走势进行较为准确的评估。
 多元线性回归适用于对受到多因素影响的数据进行分析的场景。由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。本文基于数学模型，对过去一段时间某一地区的房屋出售价格等相关数据进行整理，利用多元线性回归的方法对数据进行分析，预测该地区未来的房价走势。
 关键词：多元线性回归；房价预测；数据分析；
 
引言
 对未来房价进行预测，在一定程度上影响着社会经济的发展。广义上讲，精确的房价预测有助于国家对市场房价走势的宏观调控，小范围来讲，未来房价预测是企业战略规划的一部分，对于消费者而言，房价预测为个人经济的合理规划起到了积极作用。由于房屋售价与多因素有关，并且房屋价格与影响房价的一些因素存在线性关系，所以选取多元线性回归模型研究该问题较为合适。
 本次课题研究通过对某段时间某地区的已售房价数据进行线性回归分析，探索影响房价高低的主要因素，并对这些影响因素的影响程度进行分析，利用分析得到的数据，对未来房价的趋势和走向进行预测。
 
线性回归理论基础
 一元线性回归是分析只有一个自变量（自变量x和因变量y）线性相关关系的方法。一元线性回归分析的数学模型为：y = a+bx+ε。
 使用偏差平方和分别对参数a和参数b求偏导，可以得到线性模型的未知参数a、b的最小二乘估计值，其中，偏差平方和定义为∑(yi-a-bXi)2，a和b的唯一解如图所示。
 
 
** 参数的最小二乘估计**
 为了方便回归效果显著性检验，根据b的估计，引入LXX、LYY、LXY三个数学符号，这三个数学符号定义如图所示。
 
 
图 LXX、LYY、LXY的数学定义
 在现实问题研究中，因变量的变化往往受几个重要因素的影响，此时就需要用两个或两个以上的影响因素作为自变量来解释因变量的变化，这就是多元回归。也就是说，当多个自变量与因变量之间是线性关系时，所进行的回归分析就是多元性回归。多元线性回归的数学模型为：y=β0+β1X1+β2X2+…++βpXp+ε。使用残差平方和分别对参数βi（i=0,1,…,p）求偏导，可以得到线性模型的未知参数βi（i=0,1,…,p）的估计值，β矩阵的估计值如图所示。
 
 回归效果的显著性检验
 对平面上杂乱无章的点，利用最小二乘法求解出的线性回归方程是毫无意义的，线性回归反映出的趋势描述是否合理，需要一个数量指标来度量。
 数据总的波动可以用离差平方和LYY来描述。它表示y的各次离差yi-y ̅的平方和。LYY数值越大，说明yi数值波动越大，也就是越分散。离差平方和LYY可以分解为回归直线上y的离差平方和U以及yi与回归直线上的y间的差的平方和Q。其中，U是由于x对y的线性相关关系引起的y的分散性，Q是由随机误差引起的分散性。yi-y ̅分解如图2-4所示。在总和中，U所占比重越大，说明随机误差所占的比重越小，回归效果越显著。故此，可以使用决定系数R2来度量线性回归效果是否显著，R2作为拟合优度，表示用直线来拟合数据的好坏，R2等于U/Lyy。
 
 
R2开方后的结果为皮尔逊相关系数，皮尔逊（Pearson）相关系数可以用来衡量两个数据集合是否在一条线上面，从而衡量定距变量间的线性关系。相关系数的绝对值越大，相关性越强；相关系数越接近于1或-1，相关度越强，相关系数越接近于0，相关度越弱。当|r|>=0.8时，x和y强相关，当|r|<0.3时，x和y弱相关。皮尔逊相关系数定义为如图所示。
 
 对于一元线性回归模型，线性回归模型效果的显著性可以通过假设检验问题H0：b=0；H1：b≠0进行判断，检验方法包括F检验法和t检验法。F检验属于回归方程显著性检验，是检验x与y是否相关的检验方法。t检验是回归系数显著性检验，是检验变量x是否有用的方法。H0成立时，两种检验方法定义如图2-6、图2-7所示。H0不成立时，对于给定的显著性水平α，当F>F1-α(1,n-2)时，回归效果显著。当|t|>t1-α/2(n-2)时，认为回归系数影响显著，否则回归系数的效果不显著。
 
 一元线性回归F检验法
 
 一元线性回归t检验法
 对于多元线性回归模型，回归效果的显著性可以使用F检验法通过假设检验问题H0：β0=β1=β2=…=βp=0；H1：βi（i=0,1,…,p）不全为0进行判断，H0成立时，F检验方法定义如图2-8所示。H0不成立时，对于给定的显著性水平α，当F>F1-α(p,n-p-1)时，回归效果显著。
 
 回归系数的显著性检验可以使用t检验法通过假设检验问题H0：βi=0；H1：βi≠0进行判断。H0成立时，t检验方法定义如图2-9所示。H0不成立时，对于给定的显著性水平α，当|t|>t1-α/2(n-p-1)时，认为回归系数影响显著，否则回归系数的效果不显著。
 
 基于机器学习的线性回归与梯度下降
 机器学习横跨计算机科学、工程技术和统计学等多个学科，渗透到了人们生产和生活中的各个领域，被广泛应用于各行各业之中，在当今世界激烈的竞争中，掌握和理解机器学习的基础模型和基本方法是非常有必要的。
 机器学习中的线性回归模型以数理统计的线性回归模型为基础，它用一条直线对数据点进行拟合，在机器学习中，回归问题的求解过程就是寻找最佳拟合参数集的过程，也就是寻找满足使得估计值与实际值的方差最小时的参数解，这个过程用到了损失函数，损失函数定义如图2-10所示。利用损失函数，可以求解最佳拟合参数集。利用损失函数进行求解可以采用梯度下降法。
 
 梯度下降法的计算过程就是沿梯度下降的方向求解极小值或沿梯度上升方向求解极大值。一般情况下，梯度向量为0的话说明是到了一个极值点，此时梯度的幅值也为0。采用梯度下降算法进行最优化求解时，算法迭代的终止条件是梯度向量的幅值接近0或接近一个非常小的常数阈值。梯度下降的过程如图所示。
 
 使用EXCEL进行数据分析
 

 
 
 
 选择x,y值域
 
 

 字段Multiple R代表复相关系数R，也就是R2的平方根，又称相关系数，用来衡量自变量x与y之间的相关程度的大小。
 R Square是复测定系数，也就是相关系数R的平方。
 Adjusted R Square是调整后的复测定系数R2。
 标准误差用来衡量拟合程度的大小，也用于计算与回归相关的其它统计量，此值越小，说明拟合程度越好。
 观察值是用于估计回归方程的数据的观察值个数，本次数据集共有20条数据，所以观察值为20。
 Coefficients为常数项
 
设因变量房屋售价price为y，自变量面积area为x1，bedrooms为x2，bathrooms为x3
 
所以得到方程：
 y=10072.11+345.911x1-2925.81x2+7345.392x3
 
三、用机器学习库Sklearn库重做上面的多元线性归
 
1. 不处理直接求解
 
import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt #画图
from sklearn import linear_model #线性模型
data = pd.read_csv('house_prices.csv') 
data.head() #数据展示
 
 
new_data=data.iloc[:,1:]#除掉house_id这一列
new_data.head()

 
 
new_data.corr()#相关系数矩阵，只统计数值列
 
 
"""取are、bedrooms和bathroom作为X,price为Y求线性回归。"""
x_data = new_data.iloc[:, 1:4] #are、bedrooms、bathroom对应列
y_data = new_data.iloc[:, -1] #price对应列
print(x_data, y_data, len(x_data))

 
 
# 应用模型
model = linear_model.LinearRegression()
model.fit(x_data, y_data)
print("回归系数：", model.coef_)
print("截距：", model.intercept_)
print('回归方程: price=',model.coef_[0],'*area +',model.coef_[1],'*bedrooms +',model.coef_[2],'*bathromms +',model.intercept_)
 
 
2 数据进行清洗后再求解
 
new_data_Z=new_data.iloc[:,0:]
new_data_IQR=new_data.iloc[:,0:]
def outlier_test(data, column, method=None, z=2):
    
    if method == None:
        print(f'以 {column} 列为依据，使用 上下截断点法(iqr) 检测异常值...')
        print('=' * 70)
        column_iqr = np.quantile(data[column], 0.75) - np.quantile(data[column], 0.25)
        (q1, q3) = np.quantile(data[column], 0.25), np.quantile(data[column], 0.75)
        upper, lower = (q3 + 1.5 * column_iqr), (q1 - 1.5 * column_iqr)
        outlier = data[(data[column] <= lower) | (data[column] >= upper)]
        print(f'第一分位数: {q1}, 第三分位数：{q3}, 四分位极差：{column_iqr}')
        print(f"上截断点：{upper}, 下截断点：{lower}")
        return outlier, upper, lower
    
    if method == 'z':
        
        print(f'以 {column} 列为依据，使用 Z 分数法，z 分位数取 {z} 来检测异常值...')
        print('=' * 70)    
        mean, std = np.mean(data[column]), np.std(data[column])
        upper, lower = (mean + z * std), (mean - z * std)
        print(f"取 {z} 个 Z分数：大于 {upper} 或小于 {lower} 的即可被视为异常值。")
        print('=' * 70)
        outlier = data[(data[column] <= lower) | (data[column] >= upper)]
        return outlier, upper, lower

outlier, upper, lower = outlier_test(data=new_data_Z, column='price', method='z')
outlier.info(); outlier.sample(5)

#这里简单丢弃即可
new_data_Z.drop(index=outlier.index, inplace=True)
 
 
outlier, upper, lower = outlier_test(data=new_data_IQR, column='price')
outlier.info(); outlier.sample(5)

# 这里简单的丢弃即可
new_data_IQR.drop(index=outlier.index, inplace=True)

 
 
print("原数据相关性矩阵")
new_data.corr()
 
 
print("z方法处理的相关性矩阵")
new_data_Z.corr()
 
 
print("IQR方法处理的数据相关性矩阵")
new_data_IQR.corr()
 
 
x_data = new_data_Z.iloc[:, 1:4]
y_data = new_data_Z.iloc[:, -1]
# 应用模型
model = linear_model.LinearRegression()
model.fit(x_data, y_data)
print("回归系数：", model.coef_)
print("截距：", model.intercept_)
print('回归方程: price=',model.coef_[0],'*area +',model.coef_[1],'*bedrooms +',model.coef_[2],'*bathromms +',model.intercept_)
 
 
x_data = new_data_IQR.iloc[:, 1:4]
y_data = new_data_IQR.iloc[:, -1]
# 应用模型
model = linear_model.LinearRegression()
model.fit(x_data, y_data)
print("回归系数：", model.coef_)
print("截距：", model.intercept_)
print('回归方程: price=',model.coef_[0],'*area +',model.coef_[1],'*bedrooms +',model.coef_[2],'*bathromms +',model.intercept_)
 
 
3.3 对比
 
不做任何数据处理： price= 345.911018840024 *area + 
-2925.806324666705 *bedrooms + 7345.391713693825 *bathromms    + 10072.107046726742

Z方式清理数据：price= 226.4211697383351 *area + 49931.50311720713 *bedrooms + -12224.71724496588 *bathromms + 64356.04135007458

IQR方式清理数据：price= 242.6111551782956 *area + 41547.43068790577 *bedrooms + -6415.78250090158 *bathromms + 58018.13845504692
 
四、参考文献
 
https://blog.csdn.net/qq_55691662/article/details/120960932
 
https://blog.csdn.net/weixin_43196118/article/details/108462140
 
https://blog.csdn.net/m0_51120713/article/details/120969812
 
https://blog.csdn.net/weixin_43196118/article/details/108462140

 
多元线性模型最小一乘回归
 
一、算法目的
二、算法推导
2.1目标函数的表示形式
2.2构造一个线性规划
2.3一个证明：揭示(2)和(6)之间的关系
2.4回归算法
  
三、实际案例与python编程计算
3.1引入数据集
3.2计算$\hat{\beta }$
  
四、参考文献

 
一、算法目的
 
       多元线性模型的矩阵形式：
 $$Y=X\beta +\epsilon \text{\hspace{1em}(1)}$$
        其中，$$\begin{gathered} Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right],X=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{1p}\\ 1& x_{21}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& \cdots & x_{np}\end{array}\right] \\ \beta =\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right],\epsilon =\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_n\end{array}\right] \end{gathered}$$
        对于以上多元线性模型我们要通过最小一乘法，也就是要求出使得残差的绝对值的和最小的参数$\beta$的估计$\hat{\beta }$。
 
二、算法推导
 
2.1目标函数的表示形式
 
       我们的目标函数为：
 $$\min Q(\beta )=\sum _{i=1}^n|y_i-x_i^{\prime }\beta |\text{\hspace{1em}(2)}$$
        其中，$x_i^{\prime }=[1,x_{i1},\cdots ,x_{ip}]$。
 
2.2构造一个线性规划
 
       我们知道，任何一个向量都可以表示为两个非负向量的差，为了符合线性规划中对未知变量保持非负的习惯与一般做法。我们令$$\beta =d_1-d_2,d_1\ge 0,d_2\ge 0\text{\hspace{1em}(3)}$$
        下面我们给出线性规划约束条件的系数矩阵和系数矩阵右端项：$$A=\left[\begin{array}{ccc}{I}_n& X& -X\\ I_n& -X& X\end{array}\right],B=\left[\begin{array}{c}Y\\ -Y\end{array}\right]\text{\hspace{1em}(4)}$$
        线性规划的变量为$(r,d_1,d_2{)}^T$，线性规划目标函数的变量系数矩阵为$(l^{\prime },0,0)$。其中，$$r=\left[\begin{array}{c}{r}_1\\ r_2\\ ⋮\\ r_n\end{array}\right],l={\left[\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right]}_{n\times 1}\text{\hspace{1em}(5)}$$
        综上，线性规划表示为：
 $$\{\begin{array}{l}\min (l^{\prime },0,0)(r,d_1,d_2{)}^T={\sum }_{i=1}^n{r}_i\\ A\left[\begin{array}{c}r\\ d_1\\ d_2\end{array}\right]\ge B,\left[\begin{array}{c}r\\ d_1\\ d_2\end{array}\right]\ge 0\end{array}\text{\hspace{1em}(6)}$$
 
2.3一个证明：揭示(2)和(6)之间的关系
 
       下面我们来证明，当线性规划$(6)$求得最优解$(r,d_1,d_2{)}^T$时，$\beta =d_1-d_2$就是目标函数(2)的解。并且线性规划$(6)$的最小值就是目标函数$(2)$的最小值，也即$$\min Q=\sum _{i=1}^n{r}_i\text{\hspace{1em}(7)}$$
 证明：
        线性规划(6)的第一个约束条件可以写成$r+X{d}_1-X{d}_2\ge Y,r-X{d}_1+X{d}_2\ge -Y$。也即$r\ge Y-X\beta ,r\ge -(Y-X\beta )$。也就是$$r_i\ge |y_i-x_i^{\prime }\beta |,i=1,2,\cdots ,n\text{\hspace{1em}(8)}$$
 
       假设当线性规划$(6)$求得最优解$(r,d_1,d_2{)}^T,\beta =d_1-d_2$时，存在$\tilde{\beta }=\widetilde{d_1}-\widetilde{d_2}$，使得$Q(\tilde{\beta })<Q(\beta )$，其中$\widetilde{d_1}\ge 0,\widetilde{d_2}\ge 0$。令$\tilde{r}=(\widetilde{r_1},\widetilde{r_2},\cdots ,\widetilde{r_n})$，不妨令其中$\widetilde{r_i}=|y_i-x_i^{\prime }\tilde{\beta }|,i=1,2,\cdots ,n$.那么
 $$A\left[\begin{array}{c}\tilde{r}\\ \widetilde{d_1}\\ \widetilde{d_2}\end{array}\right]\ge B,\left[\begin{array}{c}\tilde{r}\\ \widetilde{d_1}\\ \widetilde{d_2}\end{array}\right]\ge 0\text{\hspace{1em}(9)}$$
        所以$(\tilde{r},\widetilde{d_1},\widetilde{d_2}{)}^T$也满足线性规划$(6)$的约束条件，并且$(l^{\prime },0,0)(\tilde{r},\widetilde{d_1},\widetilde{d_2}{)}^T={\sum }_{i=1}^n\widetilde{r_i}={\sum }_{i=1}^n|y_i-x_i^{\prime }\tilde{\beta }|=Q(\tilde{\beta })<Q(\beta )$，由$(8)$知，$Q(\beta )={\sum }_{i=1}^n|y_i-x_i^{\prime }\beta |\le {\sum }_{i=1}^n{r}_i$，所以${\sum }_{i=1}^n\widetilde{r_i}<{\sum }_{i=1}^n{r}_i$，这与$(r,d_1,d_2{)}^T$是线性规划的最优解矛盾，所以第一部分得证。
         从上述的推导我们可以得知，在$d_1$和$d_2$已经确定的情况下，$|y_i-x_i^{\prime }\beta |,i=1,2,\cdots ,n$均为定值，由$(8)$可知，为了使得线性规划$(6)$的目标函数最小，必须取$r_i=|y_i-x_i^{\prime }\beta |$。所以$(7)$成立。
 
2.4回归算法
 
$step(1)$：输入回归数据$X$和$Y$
 $step(2)$：根据$2.2$中的$(4),(5),(6)$计算出线性规划的相关矩阵$A,B,C$其中$C=(l^{\prime },0_{(p+1)\times 1},0_{(p+1)\times 1})$。
 $step(3)$：将$A,B,C$带入单纯形法模块，求得最优解$(r,d_1,d_2{)}^T$
 $step(4)$：得到$\hat{\beta }=d_1-d_2$
 
注：关于单纯形法原理本文不再赘述。
 
三、实际案例与python编程计算
 
3.1引入数据集
 
       我们以著名的$Hald$数据为例。
 
 

 图3.1.1
 
3.2计算$\hat{\beta }$
 
       下面给出计算最小一乘回归方程完整$python$源代码：
 
import pandas as pd
import numpy as np
from scipy import optimize
#多元线性模型的最小一乘回归(不等式约束条件)
#导入数据
dataset1=pd.read_excel('Hald.xlsx')
dataset2=pd.read_excel('Hald.xlsx')

#计算X,Y
Y=dataset1['Y'].values
dataset2['Y']=1
X=dataset2.values

#计算A,B,C
n=len(X)
B=np.hstack((Y,-Y))
In=np.eye(n)
A1=np.hstack((In,X,-X))
A2=np.hstack((In,-X,X))
A=np.vstack((A1,A2))
p=len(X[0])-1
l=np.ones(n,np.int)
o=np.zeros(p+1,np.int)
C=np.hstack((l,o,o))

#使用optimize包的linprog函数求解线性规划
r = optimize.linprog(C,A_ub=-A,b_ub=-B,bounds=(0,None))
#得到最优解x=(r,d1,d2)'
x=(r.x)
#计算β估计
beta=[]
for i in range(n,n+p+1):
    beta.append(x[i]-x[i+p+1])
    
#输出结果
print('多元线性模型的最小一乘法回归方程为：\ny=',end='')
print(beta[0],end='')
for i in range(1,p+1):
    if beta[i]>0:
        print('+{}x{}'.format(beta[i],i),end='')
    else:
        print('{}x{}'.format(beta[i], i),end='')
 
       下面给出程序运行结果：
 

 图3.2.1

 
四、参考文献
 
 
 
[1]陈希孺.最小一乘线性回归(下)[J].数理统计与管理,1989(06):48-56.
 [2]吕书龙. 最小一乘估计快速算法的研究[D].福州大学,2003.