@最小二乘法
前言本次知识来自贾俊平《统计学》第七版,初次学习,如有纰漏,请各位多多包含。诸位之意见笔者会一一听取
实现原理通过实验或调查获得测量数据后,可假定-一个函数关系,通过一定的方法确定其各项系数,即求取有关数量之间关系的经验公式从几何上看,就是要选择一条曲线,使之与所获得的数据更好地吻合.因此,求取经验公式的过程也就是曲线拟合的过程.怎样才能正确地获得与实验数据配合的最佳曲线呢?常用的方法有:散点图估计法:最小二乘法.最小二乘法是以严格的统计学理论为基础,是-种科学而可靠的曲线拟合方法.此外,还是方差分析、变量筛选、数字滤波、回归分析的数学基础.其思想如下:
设在实验中获得了自变量x与因变量y;的若干组对应数据(x,y),在使偏差平方和2[v:- f(x,)]取最小值时,找出一个已知类型的函数y= f(x)(即确定关系式中的参数).这种求解f(x)的方法称为最小二乘法。
Visual Studio Code(预装 python、jupyter)
Jupyter Notebook
Python 3.7.0(预装 numpy、pandas、matplotlib)
线性回归是利用数理统计中回归分析,来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法,在线性回归分析中,只包括一个自变量和一个因变量,且二者的关系可用一条直线近似表示,这种回归分析称为一元线性回归分析。如果回归分析中包括两个或两个以上的自变量,且因变量和自变量之间是线性关系,则称为多元线性回归分析。
最小二乘参数估计对于第i个x值,一元线性回归方程可表示为:
y=ax+b
这跟初中所学一元二次方程一致,然后其中x,y是观测值,描述直线的总数,现在的问题是,在这么多直线之中,我们应该选择哪一条作为我们的线性回归方程才会使得实验误差达到最小值,德国数学家Karl Gauss提出了最小化图的方式,也就也是使用其中的垂直方向离差平方和来估计参数a,b这一方法就称为最小二乘法,其是线性回归求方程的重要方法。
如图所示:
最小二乘化简过程详见《统计学》第七版 [贾俊平] p245
由过程化简获得方程为:
即为:
由于数据过多过大,所以使用python来使用最小二乘法对线性回归的参数进行估计。
首先要安装所需的东西(分割数据集和检测数据使用)
安装步骤为:
(1) 打开cmd
(2) 输入命令 pip install scikit-learn(此处牢记是scikit-learn而非sklearn)
根据上面参数编写代码估计线性回归参数:
分割数据集:
x_train,x_test,y_train,y_test=train_test_split(x,y,test_size=0.2,random_state=1)
求均值:
#求均值
y_avg=y_test.mean()
#求均值
x_avg=np.mean(x_test)
#分数线上部分
high=0.0
#分数线下部分
low=0.0
for i in range(len(x_train)):
#求上部分
high += (x_train[i][0]-x_avg)*(y_train[i][0]-y_avg)
#求下部分
low += (x_train[i][0]-x_avg)**2
self.b=high/low
#求出a
self.a=y_avg-self.b*x_avg
return [self.a,self.b]
最小二乘法的线性回归检验
写完之后使用python封装一下,然后开始检验线性回归。
# 引入高性能运算库
import numpy as np
#引入划分测试集和数据集
from sklearn.model_selection import train_test_split
#引入画图库
import matplotlib.pyplot as plt
#引入机器学习线性回归用作检验最小二乘法的线性回归方程
from sklearn.linear_model import LinearRegression
#引入自己设置的最小二乘法线性回归
import Linear
#导入pandas为后期
做数据处理做准备
import pandas as pd
#引入jupyter notebook魔法命令,此魔法命令可以使用鼠标悬停在图中查看线性回归的两个参数并且可以使多行图片出现在一个图片之中。
%matplotlib notebook
此处应注意的是,引用自己定义的python文件需要是使用面向对象封装的,即如上图,且与使用的文件位于同一文件夹下(与jupyter notebook(.ipynb后缀文件)处于同一文件夹,无法找到按win键然后打开jupyter文件位置,修改空白后面为文件位置(注意加上空格))
导入自己建立的参数估计方法与各种第三方库,开始回归预测。
自己自定义了一个满足线性回归的参数来检测预测结果,如上图所示,转换类型为ndarray,然后传入自己的数据作为测试,使用自定义数据并且使用自己建立的线性回归与机器学习分别进行预测,然后画图看其结果:
plt.scatter(x,y)
plt.scatter(x,y)
plt.plot(x,y_ppp,color="r")
plt.plot(x,y_hat)
图片如下:
可以看出,我们自己使用最小二乘法参数估计实现的线性回归跟机器学习所提供的是完全重合的(虽然数据集太小)。
下面使用机器学习经典案例之波士顿房价来进行线性回归预测与检验。
首先导入数据查看:
由于本次分析案例的重点不在于机器学习,所以对于CHAS这一种需要数据二值化的数据选择直接删掉。
根据数据画出图表:
由上图容易看出,其中RM和LSTAT与价格是比较明显的线性关系,所以选取其中一个进行线性回归预测都可以,在此选择RM。
首先画出RM的散点图(其实上面图表之中也有,放大有助于观察)
然后与上面步骤一致,进行预测(此处代码便不再贴出,后附完全源码),使用预测的模型来画出图,与原图进行对比。查看所拟合的线性回归方程是否与原图点的趋势相差太大。(此处因为考虑到并不是实现机器学习模型,故未剪枝)所以会存在一定误差影响(被图片中脱离线性方程周围过远的点所影响)。
使用最小二乘法参数估计所实现的线性回归与机器学习所提供的线性回归的回归方程对比:
其中红色的是使用机器学习所预测的线性回归方程。可以看出,与使用最小二乘法参数估计所求的基本一致。
对回归方程进行评价
残差:残差是因变量的观测值yi与根据线性回归方程所预测的预测值之间的差,第i个残差可以表示为:
ei=yi-y(pri)i
由图可见,由于未剪枝,些许数据还是偏离x=0较远,但是大部分数据还是处于x=0周围。因此本次数据还是一个比较满意的模型。
欧几里得距离(欧式距离)
欧几里得距离与其它几种距离(汉明距离、曼哈顿距离)是回归方程的重要评判方式。其公式为:
此处判定距离为预测值和真实值的距离,x坐标都是一致,即为
判定系数:判定系数是对估计的回归方程的拟合优度的度量。因变量y的取值不同,y的这种波动叫做变差。变差产生来自于x取值不同和误差测量等(如本次未剪枝等),其值等于实际观测值与均值之差的平方和,即
判定系数R^2 测度来了回归直线对观测数据的拟合程度。R^2越接近于1,表明回归平方占总平方和比例越大,回归直线拟合程度越好。
最小二乘法拟合:array([0.48496756])
机器学习相关拟合:array([0.48518179])
不难看出,本次回归方程与机器学习相关预测基本一致。
附:本文代码以及说明import numpy as np
#导入分割训练数据集
from sklearn.model_selection import train_test_split
class LinearRegression:
def __init__(self) -> None:
self.a=0
self.b=0
#训练函数
def fit(self,x,y):
#分割数据集
high=0
low=0
x_train,x_test,y_train,y_test=train_test_split(x,y,test_size=0.2,random_state=1)
x_avg=x.mean()
y_avg=y.mean()
for i in range(len(x_train)):
high+=(x_train[i]-x_avg)*(y_train[i]-y_avg)
low+=(x_train[i]-x_avg)**2
self.b=high/low
self.a=y_avg-self.b*x_avg
return [self.a,self.b]
def predict(self,x):
ans=[]
for i in x:
ans.append(self.a+self.b*i)
y_ppp=np.array(ans)
return y_ppp
import numpy as np
#引入划分测试集和数据集
from sklearn.model_selection import train_test_split
#引入画图库
import matplotlib.pyplot as plt
#引入机器学习线性回归用作检验最小二乘法的线性回归方程
from sklearn.linear_model import LinearRegression
#引入自己设置的最小二乘法线性回归
import Linear
#导入pandas为后期做数据处理做准备
import pandas as pd
#引入jupyter notebook魔法命令,此魔法命令可以使用鼠标悬停在图中查看线性回归的两个参数
%matplotlib notebook
#自定义一个数据集来测试
x = [170,150,168,185,189,190]
y = [62,45,60,70,72,80]
#使用numpy转换
x = np.reshape(x,newshape=(6,1))
y = np.reshape(y,newshape=(6,1))
lr=Linear.LinearRegression()
#训练线性回归模型
lr.fit(x,y)
y_ppp=lr.predict(x)
#使用SKlearn学习的相关库预测一遍
lrr = LinearRegression()
#训练数据集
lrr.fit(x,y)
# 计算y_hat
y_hat = lrr.predict(x)
#画出预测后的图片与预测后的比较,不难看出,最小二乘和SKlearn预测的数据是重合的
plt.scatter(x,y)
plt.scatter(x,y)
plt.plot(x,y_ppp,color="r")
plt.plot(x,y_hat)
#读取波士顿房价数据集
house=pd.read_csv("D://QQ//house_data.csv")
#查看前五项
house.head()
#删除不必要的列
house.drop(['RAD','TAX','CHAS'],axis=1,inplace=True)
#使用数据填充
house['CRIM'].fillna(3.55,inplace=True)
house['NOX'].fillna(0.5,inplace=True)
house['RM'].fillna(6.29,inplace=True)
house['DIS'].fillna(3.86,inplace=True)
house['PTRATIO'].fillna(18.26,inplace=True)
house['B'].fillna(356.66,inplace=True)
house['LSTAT'].fillna(12.39,inplace=True)
#为画图做准备
names=['CRIM','ZN','INDUS','NOX','RM','AGE','DIS','PTRATIO','B','LSTAT']
house.head()
#开始画图,y轴为价格
y_data=house[['MEDV']]
for i in range(len(names)):
#找到画图的位置
plt.subplot(4,4,i+1)
#散点图的数据
plt.scatter(house[names[i]],y_data,s=5)
#x轴标签
plt.xlabel(names[i])
#y轴标签
plt.ylabel('MEDV')
plt.tight_layout()
#把图表画出来
plt.show()
#将pandas的dataframe转换为ndarray作为后期预测的原始数据
x_boston=house[['RM']]
y_boston=house[['MEDV']]
x_boston=np.array(x_boston)
y_boston=np.array(y_boston)
# 再转换两个作为求相关系数的数据集
x_bos=house['RM']
y_bos=house['MEDV']
my_rho = np.corrcoef(x_bos,y_bos)
#计算皮尔逊相关系数
my_rho
#认真比对上面的十一个小图表,发现RM可以很容易看出是线性的关系,而且相关系数已经从侧面印证了这一理论,所以把它单独拿出来做一个图表
plt.scatter(x_boston,y_boston,s=5)
plt.xlabel("RM")
plt.ylabel('MEDV')
plt.tight_layout()
plt.show()
lr.fit(x_boston,y_boston)
x_new=lr.predict(x_boston)
#将预测出的回归图画到散点图中
plt.scatter(x_boston,y_boston,s=5,color = 'c')#颜色为青色
plt.xlabel("RM")
plt.ylabel('MEDV')
plt.tight_layout()
plt.plot(x_boston,x_new,"c:")#颜色为青色
plt.show()
e_boston=[y_boston[i]-x_new[i] for i in range(len(y_boston))]
e_boston=np.array(e_boston)#转化为ndarray
e_boston.shape
# 画出图表,不难看出数据多集中于y=0上下,因此最小二乘法是一个比较满意的模型
plt.scatter(x_boston,e_boston,s=5,color = 'c')
#在零区间画一条直线来观察
plt.plot(range(4,10), [x*0 for x in range(4,10)], 'r:')
lis
plt.scatter(y_boston,lis,s=3,color="c")
plt.plot(range(50), [20 for x in range(50)],'r:')
y_avg_boston=y_boston.mean()
SST=[pow(y_boston[i]-y_avg_boston,2) for i in range(len(y_boston))]
SSE=[pow((y_boston[i]-x_new[i]),2) for i in range(len(y_boston))]
R1=sum(SST)
R2=sum(SSE)
1-(R2/R1)
#使用SKlearn学习的相关库预测一遍
lrr = LinearRegression()
#训练数据集
lrr.fit(x_boston,y_boston)
# 计算y_hat
y_hat = lrr.predict(x_boston)
SSE1=[pow((y_boston[i]-y_hat[i]),2) for i in range(len(y_boston))]
R21=sum(SSE1)
1-(R21/R1)
#将预测出的回归图画到散点图中
plt.scatter(x_boston,y_boston,s=5,color = 'c')#颜色为青色
plt.xlabel("RM")
plt.ylabel('MEDV')
plt.tight_layout()
plt.plot(x_boston,x_new,"c:")#颜色为青色
plt.show()
plt.plot(x_boston,y_hat,color="r")
注:
由于%matplotlib notebook魔法命令的存在,因此个别代码会出现上下文不太连贯问题,因此如需使用代码,请务必在jupyter notebook中加上此魔法命令。并且注意结束图片,让图片不会全部出现在同一个框。
结束一个框的图片显示:
点击图片中红圈圈中的地方即可结束本次。奈何笔者才疏学浅,无法予以深度理解,如有过错之处,还请大家多多包涵。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)