【数据分析】Python数据分析之回归分析

回归分析

---

回归分析的核心在于确定变量之间存在着的函数关系

各个变量之间可以分为确定关系和非确定关系(相对关系)，我们要做的就是对这种关系进行建模和解释。

其主要流程可以解释如下：

• 收集一组包含因变量和自变量的数据
• 根据因变量和自变量之间的关系，初步设定回归模型
• 求解合理的回归系数
• 进行相关性检验，确定相关系数
• 利用模型对因变量做出预测或解释，并计算预测值的置信区间。

---

一、一元线性回归分析

表达式如下：
$$y=f(x,\theta )+\epsilon ={\beta }_0+{\beta }_{1}x+\epsilon$$
其中，$\epsilon$表示误差项，其期望$E(\epsilon )=0$，方差等于$D(\epsilon )={\sigma }^2$,${\beta }_0$为常数项，也称为截距，${\beta }_1$为斜率。

求解参数的主流方法有最小二乘法、最大似然法、矩方法。下面介绍最小二乘法。

最小二乘法(Least Squares Estimation,LSE)通过最小化误差的平方和来寻找数据的最佳匹配。

我们定义残差平方和(Residual Sum of Squares,RSS),$△y=(y-\hat{y})$表示残差：
$$Q({\beta }_0,{\beta }_1)=\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2=\sum _{i=1}^n(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i{)}^2$$
根据微积分知识，这玩意的极值点应该在导数为0的时候取得，我们对$Q$求偏导，得到：
$$\{\begin{array}{l}\frac{\sigma Q}{\sigma {\hat{\beta }}_0}=-2{\sum }_{i=1}^n(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i)\\ \frac{\sigma Q}{\sigma {\hat{\beta }}_1}=-2{\sum }_{i=1}^n(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i)x_i\end{array}$$
求解方程组得到：
$$\{\begin{array}{l}{\hat{\beta }}_0=\bar{y}-{\hat{\beta }}_1\bar{x}\\ {\hat{\beta }}_1=\frac{{\sum }_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{{\sum }_{i=1}^n(x_i-\bar{x}{)}^2}\end{array}$$
将其带入方程，即可得到最佳拟合曲线。

误差估计

SSE：
$$\sum _{i=1}^n{e}_i^2=\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$
MSE是对SSE的无偏估计量
$${\sigma }^2=\frac{SSE}{n-2}=MSE$$

Python实现一元线性回归

Step 1️⃣ 数据准备

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.linear_model import LinearRegression

# 导入数据集
iris=load_iris()
data=pd.DataFrame(iris.data)
data.columns=['sepal-length','sepal-width','petal-length','petal-width']
print(data.head())

   sepal-length  sepal-width  petal-length  petal-width
0           5.1          3.5           1.4          0.2
1           4.9          3.0           1.4          0.2
2           4.7          3.2           1.3          0.2
3           4.6          3.1           1.5          0.2
4           5.0          3.6           1.4          0.2

# 使用scikit-learn完成回归
x=data['petal-length'].values
y=data['petal-width'].values
x=x.reshape(len(x),1)
y=y.reshape(len(y),1)
clf=LinearRegression()
clf.fit(x,y)
pre=clf.predict(x)

# 绘制图形
plt.scatter(x,y,s=50)
plt.plot(x,pre,'r-',linewidth=2)
plt.xlabel("petal-length")
plt.ylabel("petal-wdith")
for idx,m in enumerate(x):
    # 绘制长条
    # 从(m,y[idx])到(m,pre[idx])
    plt.plot([m,m],[y[idx],pre[idx]],'g-')
plt.show()

step 2️⃣ 显示回归参数

print("斜率",clf.coef_)
print("截距",clf.intercept_)
print("MSE",np.mean((y-pre)**2))

斜率 [[0.41575542]]
截距 [-0.36307552]
MSE 0.04206730919499318

step 3️⃣ 进行预测

print(clf.predict([[3.9]]))

[[1.2583706]]

---

二、多元线性回归

也就是有多个参数啦。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression

# 导入数据集
d=load_boston()
data=pd.DataFrame(d.data)
data['price']=d.target
print(data.sample(5))

            0    1      2    3      4  ...      9    10      11     12  price
373  11.10810  0.0  18.10  0.0  0.668  ...  666.0  20.2  396.90  34.77   13.8
491   0.10574  0.0  27.74  0.0  0.609  ...  711.0  20.1  390.11  18.07   13.6
91    0.03932  0.0   3.41  0.0  0.489  ...  270.0  17.8  393.55   8.20   22.0
363   4.22239  0.0  18.10  1.0  0.770  ...  666.0  20.2  353.04  14.64   16.8
322   0.35114  0.0   7.38  0.0  0.493  ...  287.0  19.6  396.90   7.70   20.4

[5 rows x 14 columns]

多元线性回归

y=d.target
x=d.data
clf=LinearRegression()
from sklearn.model_selection import train_test_split
# 分割训练集
x_train,x_test,y_train,y_test=train_test_split(x,y)
clf.fit(x_train,y_train)
print("多元线性回归模型参数",clf.coef_)
print("多元线性回归模型常数项",clf.intercept_)
print("预测",clf.predict([x_test[0]]))

多元线性回归模型参数 [-1.11747256e-01  4.05201935e-02 -6.69439553e-04  3.34919157e+00
 -1.83818082e+01  3.95199967e+00 -9.12733246e-03 -1.31523502e+00
  2.44628300e-01 -1.08309725e-02 -1.00522555e+00  7.56771086e-03
 -4.23492114e-01]
多元线性回归模型常数项 36.524193135861886
预测 [15.9054318]

模型分析

y_predict=clf.predict(x_test)
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

print("预测值均方误差",mean_squared_error(y_test,y_predict))
print("R2得分",r2_score(y_test,y_predict))
print("回归得分",clf.score(x_test,y_test))

print("各个特征间的系数矩阵",clf.coef_)
print("影响房价的特征排序",np.argsort(clf.coef_))
print("影响房价的特征排序",d.feature_names[np.argsort(clf.coef_)])

预测值均方误差 26.66065801123315
R2得分 0.7170059315243467
回归得分 0.7170059315243467
各个特征间的系数矩阵 [-1.24752898e-01  4.23381228e-02  7.89030069e-03  2.76191464e+00
 -1.86055326e+01  3.76015663e+00 -3.25002550e-03 -1.49233753e+00
  3.12843628e-01 -1.40160600e-02 -8.47213267e-01  7.64996205e-03
 -5.32883469e-01]
影响房价的特征排序 [ 4  7 10 12  0  9  6 11  2  1  8  3  5]
影响房价的特征排序 ['NOX' 'DIS' 'PTRATIO' 'LSTAT' 'CRIM' 'TAX' 'AGE' 'B' 'INDUS' 'ZN' 'RAD'
 'CHAS' 'RM']

---

三、逻辑回归

如果说线性回归偏向数学，那么逻辑回归就是机器学习从统计领域借鉴来的技术

逻辑回归用来分析二分类或有序的因变量与解释变量之间的关系，算是广义上的线性回归分析方法。他在线性回归的基础上利用Sigmoid函数对事件发生的概率进行预测。

线性回归可以得到一个预测值，然后通过S函数封装后，就能得到一个概率值，再通过概率值进行分类。(上清下浊)

Sigmoid的函数能够将任意值转化为[0,1]范围内，其定义如下：
$$g(z)=\frac{1}{1+e^{-z}}$$
我们来看下函数曲线

宏观尺度

微观尺度

def sigmoid(x):
	return 1./(1.+np.exp(-x))

Python 实现

x=load_iris().data
y=load_iris().target

# 归一化
from sklearn.model_selection import train_test_split
x_train,x_test,y_train,y_test=train_test_split(x,y,test_size=0.25,random_state=0)
from sklearn.preprocessing import StandardScaler

sc=StandardScaler()
x_train=sc.fit_transform(x_train)
x_test=sc.transform(x_test)

# 进行逻辑回归
from sklearn.linear_model import LogisticRegression
classifier=LogisticRegression(random_state=0)
classifier.fit(x_train,y_train)
y_pred=classifier.predict(x_test)
# 测试准确性
print("Accuracy of LR %.3f"%classifier.score(x_test,y_test))

Accuracy of LR 0.974

我们可以看一下经过逻辑回归后的数据：

[2 1 0 2 0 2 0 1 1 1 2 1 1 1 1 0 1 1 0 0 2 1 0 0 2 0 0 1 1 0 2 1 0 2 2 1 0 2]

这是个多分类的回归。

---

四、多项式回归

适用于非线性关系。

# 线性回归
lin_reg=LinearRegression()
lin_reg.fit(X,y)
y_pre=lin_reg.predict(X)
plt.rcParams['font.family']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
ax[1].scatter(x,y)
ax[1].plot(x,y_pre,color="r")

# 多项式回归
from sklearn.preprocessing import PolynomialFeatures
POLY=PolynomialFeatures(degree=2) # 设置最多几次幂
POLY.fit(X)
x2=POLY.transform(X)
# 这个多项式回归是对x进行处理后，让其成为非线性关系
# 譬如：
print(x.shape)
print(x2.shape)
# (100,)
# (100, 3)
# 之后的操作与LR完全相同，所以Polynomial并没有作为独立的API
# 而是放在preprocessing
lin_reg2=LinearRegression()
lin_reg2.fit(x2,y)
y_pre2=lin_reg2.predict(x2)

ax[2].scatter(x,y)
# 此时的关系并不再是原先的对应了
ax[2].plot(np.sort(x),y_pre2[np.argsort(x)],color="b")

# degree调成10 后
POLY=PolynomialFeatures(degree=10) # 设置最多几次幂
POLY.fit(X)
x2=POLY.transform(X)
# 这个多项式回归是对x进行处理后，让其成为非线性关系
# 譬如：
print(x.shape)
print(x2.shape)
# (100,)
# (100, 11)
# 之后的操作与LR完全相同，所以Polynomial并没有作为独立的API
# 而是放在preprocessing
lin_reg2=LinearRegression()
lin_reg2.fit(x2,y)
y_pre2=lin_reg2.predict(x2)
ax[3].scatter(x,y)
ax[3].plot(np.sort(x),y_pre2[np.argsort(x)],color="g")
plt.show()

---

五、岭回归

岭回归(Ridge Regression)是一种专用于共线性数据分析的有偏估计回归方法，实质上是改良的最小二乘法。通过放弃无偏性，降低部分信息为代价，使回归系数更加可靠，对病态数据的耐受性远高于最小二乘法。

通常的岭回归是在顺势函数中加入L2正则项：
$$L(\theta )=\frac{1}{N}\sum _{i=1}^N(f(x_i;\theta )-y_i{)}^2+\frac{\lambda }{2}||\theta |{|}^2$$
其中，$||\theta ||$表示向量$\theta$的$L2$范数。岭回归的$R^2$往往会小于线性回归，但其具有更强的泛化能力，也能解决线性回归汇总的不可逆问题。

from sklearn.linear_model import Ridge,RidgeCV
# Ridge CV是广义交叉验证的岭回归
X,y=load_iris(return_X_y=True)
x=X[:,1].reshape(len(X),-1)
y=X[:,0].reshape(len(X),-1)

model=Ridge(alpha=0.5)
model1=RidgeCV(alphas=[0.1,1.0,10.0]) # cross validation

model.fit(x,y)
model1.fit(x,y)

print("系数矩阵",model.coef_)
print("线性回归模型",model)
print("CV最优alpha值",model1.alpha_)

# 模型预测
pre=model.predict(x)
plt.scatter(x,y)
plt.plot(x,pre)
plt.show()

系数矩阵 [[-0.2194842]]
线性回归模型 Ridge(alpha=0.5)
CV最优alpha值 10.0

---

六、Lasso回归

$L2$正则只能削弱影响，而不能剔除变量。Lasso(Least Absolute Shrinkage and Selection Operator)模型将惩罚项换为了$L1$正则，从而达到剔除变量的作用。

from sklearn.metrics import r2_score
# 产生一些稀疏数据
np.random.seed(42)
n_samples,n_features=50,100

X=np.random.randn(n_samples,n_features)
coef=3*np.random.randn(n_features) # 每个特征对应一个系数

inds=np.arange(n_features)
np.random.shuffle(inds)
coef[inds[10:]]=0 # 随机将向量中的10个变为0 稀疏化
y=np.dot(X,coef)

# 添加高斯噪声
y+=0.01*np.random.normal(size=n_samples)

# 划分数据集
n_samples=X.shape[0]
X_train,y_train=X[:n_samples//2],y[:n_samples//2]
X_test,y_test=X[n_samples//2:],y[n_samples//2:]

# 训练Lasson模型
from sklearn.linear_model import Lasso
alpha=0.1
lass=Lasso(alpha=alpha)
y_pre=lass.fit(X_train,y_train).predict(X_test)
r2_score_lasso=r2_score(y_test,y_pre)

print("R^2 socre",r2_score_lasso)
plt.plot(lass.coef_,color='gold')
plt.title("Lasso R^2 %s"%r2_score_lasso)
plt.show()

---

回归案例–波士顿数据处理

import numpy as np
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression,SGDRegressor,Ridge,Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error,accuracy_score
from sklearn.model_selection import train_test_split

def linearModel():

    # 1. load dataset
    X,y=load_boston(return_X_y=True)
    X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.25)

    # 2. Standard
    std=StandardScaler()
    X_train=std.fit_transform(X_train)
    X_test=std.transform(X_test)

    # 添加维度
    # y=y[:,np.newaxis]
    y_test=y_test[:,np.newaxis]
    y_train=y_train[:,np.newaxis]

    std_y=StandardScaler()
    y_train=std_y.fit_transform(y_train)
    y_test=std_y.transform(y_test)

    # 3. training

    lr=LinearRegression()
    sgd=SGDRegressor()
    rid=Ridge(alpha=0.5)
    las=Lasso(alpha=0.5)

    lr.fit(X_train,y_train)
    y_lr_pre=lr.predict(X_test)
    # 还原真实值
    # 这是因为我们做处理用的都是标准化
    # 所以最后跟原始数据比较需要反标准化
    y_lr_pre=std_y.inverse_transform(y_lr_pre)

    sgd.fit(X_train,y_train)
    y_sgd_pre=sgd.predict(X_test)
    y_sgd_pre=std_y.inverse_transform(y_sgd_pre)

    rid.fit(X_train, y_train)
    y_rid_pre = rid.predict(X_test)
    y_rid_pre = std_y.inverse_transform(y_rid_pre)

    las.fit(X_train, y_train)
    y_las_pre = las.predict(X_test)
    y_las_pre = std_y.inverse_transform(y_las_pre)

    # 4. evaluation
    y_test=std_y.inverse_transform(y_test)
    print("线性回归的均方误差为: ",mean_squared_error(y_test,y_lr_pre))

    print("梯度下降的均方误差为: ",mean_squared_error(y_test,y_sgd_pre))

    print("岭回归的均方误差为: ",mean_squared_error(y_test,y_rid_pre))

    print("Lasson均方误差为: ",mean_squared_error(y_test,y_las_pre))

linearModel()

线性回归的均方误差为:  24.474589501317524
梯度下降的均方误差为:  23.82902862765662
岭回归的均方误差为:  24.433234356526174
Lasson均方误差为:  61.95307932912878