回归分析实例

回归分析研究自变量（也称为预测变量）与因变量（也称为响应变量）之间的关系。它旨在通过建立数学模型来描述和预测因变量的变化，基于自变量的取值。

在回归分析中，我们通常假设自变量和因变量之间存在一种函数关系。这个函数可以是线性的，如线性回归模型，也可以是非线性的，如多项式回归模型、指数回归模型等。通过收集一组自变量和对应的因变量的观测值，我们可以使用回归分析来找到最佳拟合该函数关系的模型参数。

在建立回归模型时，我们常常使用最小二乘法来估计模型参数。最小二乘法通过最小化预测值与观测值之间的差异来找到最佳的参数，使得模型与实际数据的拟合度最好。

回归分析概述

相关分析是研究两个或两个以上的变量之间相关程度及大小的一种统计方法

回归分析是寻找存在相关关系的变量间的数学表达式，并进行统计推断的一种统计方法

在对回归分析进行分类时，主要有两种分类方式：

• 根据变量的数目，可以分类一元回归、多元回归

• 根据自变量与因变量的表现形式，分为线性与非线性

所以，回归分析包括四个方向：一元线性回归分析、多元线性回归分析、一元非线性回归分析、多元非线性回归分析。

回归分析一般步骤

回归方程定义

一元线性回归分析

• 因变量(dependent variable)：被预测或被解释的变量，用y表示

• 自变量(independent variable)：预测或解释因变量的一个或多个变量，用x表示

• 对于具有线性关系的两个变量，可以用一个方程来表示它们之间的线性关系

• 描述因变量y如何依赖于自变量x和误差项ε的方程称为回归模型。对于只涉及一个自变量的一元线性回归模型可表示为：

• y叫做因变量或被解释变量
• x叫做自变量或解释变量
• β0 表示截距
• β1 表示斜率
• ε表示误差项，反映除x和y之间的线性关系之外的随机因素对y的影响

一元回归例子

• 人均收入是否会显著影响人均食品消费支出

• 贷款余额是否会影响到不良贷款

• 航班正点率是否对顾客投诉次数有显著影响

回归方程

描述因变量y的期望值如何依赖于自变量x的方程称为回归方程。根据对一元线性回归模型的假设，可以得到它的回归方程为：

误差项定义

最小二乘法推导与求解

对于回归直线，关键在于求解参数，常用高斯提出的最小二乘法，它是使因变量的观察值y与估计值之间的离差平方和达到最小来求解。

回归方程求解小例子

70年代世界制造业总产量与世界制成品总出口量的变化关系如表

利用回归直线进行估计与预测

• 点估计：利用估计的回归方程，对于x的某一个特定的值 ，求出y的一个估计值 就是点估计

• 区间估计：利用估计的回归方程，对于x的一个特定值 ，求出y的一个估计值的区间就是区间估计

估计标准误差的计算

为了度量回归方程的可靠性，通常计算估计标准误差。它度量观察值回绕着回归直线的变化程度或分散程度。

估计平均误差：

• 公式中根号内的分母是n-2,而不是n，因而自由度为n-2。

• 估计标准误差越大，则数据点围绕回归直线的分散程度就越大，回归方程的代表性越小。

• 估计标准误差越小，则数据点围绕回归直线的分散程度越小，回归方程的代表愈大，其可靠性越高。

置信区间估计

在1－α置信水平下预测区间为

某企业从有关资料中发现广告投入和产品销售有较密切的关系。近年该企业广告费和销售额资料见表10-3，若2003年广告费为120万元，请用一元线性回归求2003年产品销售额的置信区间与预测区间（α＝0.05）

影响区间宽度的因素

• 置信水平 (1 - α)，区间宽度随置信水平的增大而增大

• 数据的离散程度Se，区间宽度随离程度的增大而增大

• 样本容量，区间宽度随样本容量的增大而减小

• X0与X均值之间的差异，随着差异程度的增大而增大

回归直线拟合优度

回归直线与各观测点的接近程度称为回归直线对数据的拟合优度

总平方和 SST

总平方和可以分解为回归平方和、残差平方和两部分：SST＝SSR+SSE

反映因变量的 n 个观察值与其均值的总离差

回归平方和 SSR

回归平方和SSR反映了y的总变差中，由于x与y之间的线性关系引起的y的变化部分

残差平方和 SSE

残差平方和SSE反映了除了x对y的线性影响之外的其他因素对y变差的作用，是不能由回归直线来解释的y的变差部分

判定系数

判定系数是衡量回归方称好坏的重要因素,回归平方和占总平方和的比例，用R^2表示,其值在0到1之间。

• R^2 == 0：说明y的变化与x无关，x完全无助于解释y的变差

• R^2 == 1：说明残差平方和为0，拟合是完全的，y的变化只与x有关

（可解释的部分占总体的多少，理想情况下R^2接近于0）

调整的多重判定系数

用样本容量n和自变量的个数k去修正R^2得到，应避免增加自变量而高估R^2

显著性检验

显著性检验的主要目的是根据所建立的估计方程用自变量x来估计或预测因变量y的取值。当建立了估计方程后，还不能马上进行估计或预测，因为该估计方程是根据样本数据得到的，它是否真实的反映了变量x和y之间的关系，则需要通过检验后才能证实。

根据样本数据拟合回归方程时，实际上就已经假定变量x与y之间存在着线性关系，并假定误差项是一个服从正态分布的随机变量，且具有相同的方差。但这些假设是否成立需要检验

显著性检验包括两方面：

• 线性关系检验

线性关系检验是检验自变量x和因变量y之间的线性关系是否显著，或者说，它们之间能否用一个线性模型来表示。

将均方回归 (MSR)同均方残差 (MSE)加以比较，应用F检验来分析二者之间的差别是否显著。

• 均方回归：回归平方和SSR除以相应的自由度(自变量的个数K)

• 均方残差：残差平方和SSE除以相应的自由度(n-k-1)

H0：β1=0 所有回归系数与零无显著差异，y与全体x的线性关系不显著

计算检验统计量F：

• 回归系数检验

回归系数显著性检验的目的是通过检验回归系数β的值与0是否有显著性差异，来判断Y与X之间是否有显著的线性关系。若β=0，则总体回归方程中不含X项(即Y不随X变动而变动)，因此，变量Y与X之间并不存在线性关系；若β≠0，说明变量Y与X之间存在显著的线性关系。

计算检验的统计量：

线性关系检验与回归系数检验的区别

线性关系的检验是检验自变量与因变量是否可以用线性来表达，而回归系数的检验是对样本数据计算的回归系数检验总体中回归系数是否为0

• 在一元线性回归中，自变量只有一个，线性关系检验与回归系数检验是等价的

• 多元回归分析中，这两种检验的意义是不同的。线性关系检验只能用来检验总体回归关系的显著性，而回归系数检验可以对各个回归系数分别进行检验

多元与曲线回归分析

经常会遇到某一现象的发展和变化取决于几个影响因素的情况，也就是一个因变量和几个自变量有依存关系的情况，这时需用多元线性回归分析。

• 多元线性回归分析预测法，是指通过对两上或两个以上的自变量与一个因变量的相关分析，建立预测模型进行预测和控制的方法

• 多元线性回归预测模型一般式为：

使用statsmodels模块进行回归分析

一元回归

# 使用statsmodels模块进行回归分析
import numpy as np 
import matplotlib.pyplot as plt
import statsmodels.api as sm

nsample = 20
# 0~10之间20个样本数据
x = np.linspace(0,10,nsample)
# 一元线性回归
# 添加一列'1'  
X = sm.add_constant(x)
# 第一列全为1 ，第二列为真实值方便计算
X
#β0,β1分别设置成2,5
beta = np.array([2,5])
beta
# 设置误差项(服从高斯分布)
e = np.random.normal(size=nsample)
# 使数据上下抖动
e
# 实际值 y
y = np.dot(X, beta)+e
y
# 最小二乘法
model = sm.OLS(y,X)
# 拟合数据
res = model.fit()
# 回归系数
# 原始指定的β0,β1分别、2,5
# 真实计算得到β0,β1分别1.83214489, 5.08118177
res.params
# 显示最终数据的全部结果
res.summary()

# 拟合的估计值
y_ = res.fittedvalues
# 实际的估计值
y_
fig,ax = plt.subplots(figsize=(8,6))
# 原始数据
ax.plot(x,y,'o',label='data')
# 拟合数据
ax.plot(x,y_,'r--',label='test')
ax.legend(loc='best')
plt.show()

高阶回归

# 高阶回归  Y=5+2⋅X+3⋅X^2
nsamples = 50
x = np.linspace(0,10,nsamples)
print('第一个自变量x:',x)
X = np.column_stack((x,x**2))
print('第二个自变量X:',X)
X = sm.add_constant(X)
print('整合后X:',X)
results.summary()

y_fitted = results.fittedvalues
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label='data')
ax.plot(x, y_fitted, 'r--.',label='OLS')
ax.legend(loc='best')
plt.show()

分类变量

假设分类变量有3个取值（a,b,c）,比如考试成绩有3个等级。那么a就是（1,0,0），b（0,1,0），c(0,0,1),这个时候就需要3个系数β0,β1,β2，也就是β0x0+β1x1+β2x2

nsample = 50
groups = np.zeros(nsample,int)
groups[20:40] = 1
groups[40:] = 2
# 得到指定数据格式
dummy = sm.categorical(groups, drop=True)
dummy

#Y=5+2X+3Z1+6⋅Z2+9⋅Z3.

x = np.linspace(0, 20, nsample)
X = np.column_stack((x, dummy))
X = sm.add_constant(X)
beta = [5, 2, 3, 6, 9]
e = np.random.normal(size=nsample)
y = np.dot(X, beta) + e
result = sm.OLS(y,X).fit()
result.summary()
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, result.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best')
plt.show()

回归分析实例

!pip install -i https://mirrors.aliyun.com/pypi/simple plotly
import pandas as pd
import statsmodels.api as sm
# 动态可视化包
from plotly.offline import init_notebook_mode, iplot
import plotly.graph_objs as go

init_notebook_mode(connected=True)
# 加载数据
data = pd.read_csv('./data/weight.csv')
data = data[data.height > 120]
data

# 将数据以散点图形式动态显示
trace = go.Scatter(x=data.height, y=data.weight, mode='markers')

layout = go.Layout(
    width=900,
    height=500,
    xaxis=dict(title='Height',titlefont=dict(family='Consolas, monospace',size=15)),
    yaxis=dict(title='Weight',titlefont=dict(family='Consolas, monospace',size=15))
)

data2 = [trace]
fig = go.Figure(data=data2, layout=layout)

iplot(fig)

X = data.height
y = data.weight
X = sm.add_constant(X)

model = sm.OLS(y, X).fit()

model.summary()

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(data.height, y, 'o', label="data")
ax.plot(data.height, model.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best')
plt.show()

scikit-learn进行回归分析

汽车价格预测

数据集简介

主要包括3类指标:

• 汽车的各种特性
• 保险风险评级：(-3, -2, -1, 0, 1, 2, 3)
• 每辆保险车辆年平均相对损失支付

类别属性

• make: 汽车的商标（奥迪，宝马。。。）
• fuel-type: 汽油还是天然气
• aspiration: 涡轮
• num-of-doors: 两门还是四门
• body-style: 硬顶车、轿车、掀背车、敞篷车
• drive-wheels: 驱动轮
• engine-location: 发动机位置
• engine-type: 发动机类型
• num-of-cylinders: 几个气缸
• fuel-system: 燃油系统

连续指标

• bore: continuous from 2.54 to 3.94.
• stroke: continuous from 2.07 to 4.17.
• compression-ratio: continuous from 7 to 23.
• horsepower: continuous from 48 to 288.
• peak-rpm: continuous from 4150 to 6600.
• city-mpg: continuous from 13 to 49.
• highway-mpg: continuous from 16 to 54.
• price: continuous from 5118 to 45400.

# loading packages
import numpy as np
import pandas as pd
from pandas import datetime

# data visualization and missing values
import matplotlib.pyplot as plt
import seaborn as sns # advanced vizs
# 对缺失值进行可视化展示
import missingno as msno # missing values
%matplotlib inline

# stats
from statsmodels.distributions.empirical_distribution import ECDF
from sklearn.metrics import mean_squared_error, r2_score

# machine learning
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Lasso, LassoCV 
from sklearn.model_selection import train_test_split, cross_val_score 
from sklearn.ensemble import RandomForestRegressor
# 设置随机种子（控制随机因素）
seed = 123

# importing data ( ? = missing values)
data = pd.read_csv("Auto-Data.csv", na_values = '?')
data.columns
data.dtypes

# first glance at the data itself
print("In total: ",data.shape)
data.head(5)

一个比较小的数据集:只有 205 条数据, 26个特征

缺失值填充

# 1.数据量大的情况下可以直接将缺失值去掉
# 2.使用平均值、众数、中位数填充缺失值
# 3.构建回归模型预测缺失值
sns.set(style = "ticks")

msno.matrix(data)
# 白色代表缺失值
#https://github.com/ResidentMario/missingno

# missing values in normalied-losses
# "pd.isnull()"是用于检查值是否为缺失值的函数。在这里，它被用于判断"data['normalized-losses']"列中的每个元素是否为缺失值，并返回前五条数据
data[pd.isnull(data['normalized-losses'])].head()

sns.set(style = "ticks")
plt.figure(figsize = (12, 5)) 
c = '#366DE8'

# ECDF
"""
ECDF 可以用来可视化和分析数据的累积分布情况。具体来说，它用于估计给定数据集中每个值的累积概率。对于一个给定的观测值 x，ECDF 提供的结果是小于或等于 x 的观测值占总观测值数量的比例。

计算 ECDF 的步骤非常简单：

    对数据进行排序，从小到大排列。
    对于每个数据点，计算该数据点在排序后数据中的位置（排名），并除以总观测值数量。
    绘制每个数据点的值在 x 轴上，以它的 ECDF 值在 y 轴上，形成折线图。

通过观察 ECDF 的图形，你可以看到数据的分布情况，比如数据集中的集中区域、分布的形状等。ECDF 还可以用于比较不同数据集之间的分布差异。
"""
plt.subplot(121)
cdf = ECDF(data['normalized-losses'])
plt.plot(cdf.x, cdf.y, label = "statmodels", color = c);
plt.xlabel('normalized losses'); plt.ylabel('ECDF');

# 总体分布
plt.subplot(122)
plt.hist(data['normalized-losses'].dropna(), 
         bins = int(np.sqrt(len(data['normalized-losses']))),
         color = c);

可以发现 80% 的 normalized losses 是低于200 并且绝大多数低于125.

一个基本的想法就是用中位数来进行填充，但是我们得来想一想，这个特征跟哪些因素可能有关呢？应该是保险的情况，所以我们可以分组来进行填充这样会更精确一些。

首先来看一下对于不同保险情况的统计指标:

# 使用groupby函数来对"data"数据集中的"symboling"列进行分组，并对每个分组中的"normalized-losses"列进行描述性统计。

"""

    groupby('symboling')：首先，根据"symboling"列的不同取值，将数据集分成了多个组。每个组都包含着具有相同"symboling"取值的行。

    ['normalized-losses']：然后，在每个组内，选择了"normalized-losses"列作为要进行描述性统计的列。

    describe()：最后，对每个组内的"normalized-losses"列进行了描述性统计分析。描述性统计主要包括计数（count）、均值（mean）、标准差（std）、最小值（min）、25% 分位数（25%）、中位数（50%）、75% 分位数（75%）和最大值（max）等。

"""
# 目的是为了了解不同"symboling"值下的"normalized-losses"列的分布情况
# 例如上述第一条缺失值的normalized-losses填充根据 symboling=3的所有normalized-losses列的平均值进行填充
# 上述第三条缺失值则使用symboling=1的所有normalized-losses列的平均值进行填充
data.groupby('symboling')['normalized-losses'].describe()

# 缺失值替换
# 从"data"数据集中删除含有缺失值的行，其中包括'price', 'bore', 'stroke', 'peak-rpm', 'horsepower'和'num-of-doors'这些列。
data = data.dropna(subset = ['price', 'bore', 'stroke', 'peak-rpm', 'horsepower', 'num-of-doors'])
# 定位到缺失值比较严重的这一列，对这一列使用当前组的均值进行改变
"""
首先，代码使用了groupby函数，按照'symboling'列对数据进行分组。这将数据集中的每个'symboling'值作为分组依据。

然后，我们选择了'normalized-losses'列，并使用transform函数应用一个匿名函数（lambda函数）来处理该列。lambda函数的功能是将每个分组中的缺失值（NaN）替换为该分组对应的平均值。

最后，我们将填充后的结果赋值给'data['normalized-losses']'列，使得原始数据集中的缺失值被修复。
"""
data['normalized-losses'] = data.groupby('symboling')['normalized-losses'].transform(lambda x: x.fillna(x.mean()))

print('In total:', data.shape)
data.head()

特征相关性

#计算数据集的相关系数矩阵 计算数据集的相关系数矩阵
# 相关系数矩阵描述了数据集中各个变量之间的线性关系强度和方向。
# 相关系数的取值范围在 -1 到 1 之间，其中 -1 表示强反相关，0 表示无相关性，1 表示强相关。
# data.corr()调用了pandas库的corr()函数，该函数作用是计算给定数据集的相关系数。
# 该函数会计算每对变量之间的相关系数，并返回相关系数矩阵。
cormatrix = data.corr()
cormatrix

"""
使用NumPy库来将相关系数矩阵的上三角区域（不包括对角线）设置为零。

代码中的np.tri(*cormatrix.values.shape, k=-1)是一个函数调用，它创建一个与相关系数矩阵具有相同形状的三角矩阵，
上三角部分（不包括对角线）的元素为1，下三角和对角线部分的元素为0。参数k=-1表示上三角部分相对于对角线的偏移量。

然后，使用*=运算符将相关系数矩阵与上述生成的三角矩阵进行逐元素相乘。这会将相关系数矩阵的上三角部分（不包括对角线）的值设置为零，
而保持其他部分不变。

这样做的目的可能是为了在分析和可视化相关系数矩阵时，只关注变量之间的单向相关性而不考虑重复的信息。
对于大型相关系数矩阵，这可以减少冗余和复杂性，从而更清晰地呈现相关关系。
"""

cormatrix *= np.tri(*cormatrix.values.shape, k=-1).T  #返回函数的上三角矩阵，把对角线上的置0，让他们不是最高的。
# 将对角线元素设置为0只取上三角区域元素
cormatrix

# 创建一个叠加的相关性矩阵
"""
使用Pandas库中的stack()函数。stack()函数的作用是将数据按照指定的轴（默认是最内层的轴）进行堆叠，将多层次的索引转换成一个新的二维索引。

在这里，cormatrix是一个相关性矩阵，可能呈现为一个二维的数据结构，其中包含了各个变量之间的相关性信息。
通过调用stack()函数，我们将相关性矩阵中的每个元素都转化成一个二维索引，使得这些相关性数据能够更方便地进行处理和分析。
"""
cormatrix = cormatrix.stack()
cormatrix

# 根据相关性的绝对值对相关性矩阵进行重新排序。
"""
    abs(): 这个函数用来计算相关性矩阵中每个元素的绝对值。

    sort_values(ascending=False): 这个方法对绝对值进行排序，ascending=False表示按降序排列，也就是将绝对值最大的相关性排在前面。

    index: 这个属性返回排序后的索引值，即相关性矩阵元素的位置顺序。

    reindex(): 这个方法通过指定新的索引值，对相关性矩阵进行重新排序。在这里，使用了前面得到的排序后的索引值来重新排列相关性矩阵。

    reset_index(): 这个方法将原始的索引重置为默认的整数索引，同时创建一个新的索引列。

"""
cormatrix = cormatrix.reindex(cormatrix.abs().sort_values(ascending=False).index).reset_index()
cormatrix

# 重命名列名
cormatrix.columns = ["FirstVariable", "SecondVariable", "Correlation"]
cormatrix.head(10)

Correlation值越高表明两个因素相关性越相近删除其中一个指标，对于因素相近的指标应进行合并（例如'width', 'length', 'height'可以用width*length*height代替）

data['volume'] = data.length * data.width * data.height

data.drop(['width', 'length', 'height', 
           'curb-weight', 'city-mpg'], 
          axis = 1, # 1 for columns
          inplace = True) 

# new variables
data.columns

# 计算相关矩阵
corr_all = data.corr()

# 为上部三角形生成热力图
mask = np.zeros_like(corr_all, dtype = np.bool)
mask[np.triu_indices_from(mask)] = True

# 设置matplotlib图形
f, ax = plt.subplots(figsize = (11, 9))

sns.heatmap(corr_all, mask = mask,
            square = True, linewidths = .5, ax = ax, cmap = "BuPu")      
plt.show()

price 跟这几个的相关程度比较大 wheel-base,enginine-size, bore,horsepower

"""
这段代码使用了Seaborn库中的pairplot函数来创建一个散点图矩阵，用于可视化数据集data中的多个变量之间的关系。参数hue指定了一个分类变量，这里是"fuel-type"，它将通过不同颜色来区分不同类别的数据点。参数palette指定了使用的颜色调色板，这里是"plasma"。

通过这个散点图矩阵，你可以同时查看数据集中每对变量之间的散点图和变量的分布情况。这对于发现变量之间的关联性、异常值和数据分布特征非常有用。

如果你想深入了解每个变量对之间的关系，可以观察散点图矩阵中每个小格子中的点分布情况。不同颜色的点表示了不同的燃料类型，在分析数据时可以根据这个分类信息进行更准确的解读。
"""
sns.pairplot(data, hue = 'fuel-type', palette = 'plasma')

# 查看价格和马力变量之间的关系
sns.lmplot('price', 'horsepower', data, 
           hue = 'fuel-type', col = 'fuel-type',  row = 'num-of-doors', 
           palette = 'plasma', 
           fit_reg = True);

预处理问题

如果一个特性的方差比其他的要大得多，那么它可能支配目标函数，使估计者不能像预期的那样正确地从其他特性中学习。这就是为什么我们需要首先对数据进行缩放。

对连续值进行标准化

# target为数据集data中的"price"列的引用或副本
target = data.price
# data数据集中除了'price'以外的所有列，并将它们作为回归器（自变量）添加到regressors列表中
regressors = [x for x in data.columns if x not in ['price']]
features = data.loc[:, regressors]
# data.loc[:, regressors]使用loc函数从data中选择所有行（冒号表示选择所有行）和在regressors列表中列出的列。
# 返回一个新的数据框，其中包含data中的所有行，但只包括regressors列表中指定的列作为特征。
num = ['symboling', 'normalized-losses', 'volume', 'horsepower', 'wheel-base',
       'bore', 'stroke','compression-ratio', 'peak-rpm']

# 缩放数据
"""代码standard_scaler = StandardScaler()创建了一个名为standard_scaler的StandardScaler对象。
StandardScaler是一个在机器学习中常用的数据预处理技术，用于将特征数据进行标准化处理。
标准化过程将特征数据按照均值为0、标准差为1的方式进行缩放，使得数据具有相同的尺度，有助于提高模型的性能。
"""
standard_scaler = StandardScaler()
"""
使用standard_scaler.fit_transform()方法对它们进行标准化操作，其中num是一个索引列表，表示我们要对哪些列进行标准化处理。
通过调用fit_transform()方法，我们首先对数据进行拟合（计算均值和标准差），然后对该列的每个值进行标准化转换，使其符合标准正态分布。
最后，我们将标准化后的值存储回特征集合的对应列features[num]中，以便后续的机器学习模型训练和分析
"""
features[num] = standard_scaler.fit_transform(features[num])

# glimpse
features.head()

对分类属性就行one-hot编码

# categorical vars
classes = ['make', 'fuel-type', 'aspiration', 'num-of-doors', 
           'body-style', 'drive-wheels', 'engine-location',
           'engine-type', 'num-of-cylinders', 'fuel-system']

# 参照实例 abc分类 转换为[[1,0,0],[0,1,0],[0,0,1]]
dummies = pd.get_dummies(features[classes])
features = features.join(dummies).drop(classes, 
                                       axis = 1)

# new dataset
print('In total:', features.shape)
features.head()

划分数据集

# split the data into train/test set
X_train, X_test, y_train, y_test = train_test_split(features, target, 
                                                    test_size = 0.3,
                                                    random_state = seed)
print("Train", X_train.shape, "and test", X_test.shape)

Lasso回归

一般线性回归模型的目标是最小化残差平方和，即通过拟合一个线性方程来预测目标变量。然而，在实际问题中，可能存在大量的自变量，其中一些可能对目标变量的预测能力较弱或冗余。此时，Lasso回归通过引入L1正则化 (即Lasso惩罚项)，可以将系数向量中小的权重变为0，从而实现特征选择和模型稀疏性。(多加了一个绝对值项来惩罚过大的系数，alphas=0那就是最小二乘了)

公式

代码实现

# 生成一个指数序列，作为一组用于正则化参数(alpha)的候选值
"""
首先，np.arange(2, 12)生成一个从2到11的整数序列，不包括12。这个序列表示了指数的幂次，即2的幂次。
例如，生成的序列为[2, 3, 4, 5, 6, 7, 8, 9, 10, 11]。
接下来，2. ** np.arange(2, 12)对上述生成的幂次序列进行指数运算，计算得到每个幂次对应的2的幂次方值。
使用**运算符表示指数运算，即2的x次方。
"""
alphas = 2. ** np.arange(2, 12)
"""
首先，np.empty_like(alphas)使用alphas数组的形状创建了一个空数组。
np.empty_like()函数返回一个新的、大小与给定参数相同的空数组，但不会对数组进行初始化，即数组中的元素值是不确定的。
将这个空数组赋值给变量scores，用于存储后续操作中的分数或者评估结果。
"""
scores = np.empty_like(alphas)

# 使用Lasso回归模型进行模型训练和评估
"""
使用Lasso回归模型进行模型训练和评估。
1. `for i, a in enumerate(alphas):`：遍历正则化参数数组`alphas`的每个元素，并用索引`i`和元素值`a`进行迭代。
2. `lasso = Lasso(random_state = seed)`: 创建一个Lasso回归模型对象`lasso`，并指定随机种子`random_state`来保证结果的可重复性。
3. `lasso.set_params(alpha = a)`: 设置当前迭代中Lasso模型的正则化参数为`a`。
4. `lasso.fit(X_train, y_train)`: 使用训练数据`X_train`和目标变量`y_train`对Lasso模型进行训练。
5. `scores[i] = lasso.score(X_test, y_test)`: 使用测试数据`X_test`和目标变量`y_test`对训练好的Lasso模型进行评估，
得到模型的得分，并将得分保存在`scores`数组的第`i`个位置。

通过遍历不同的正则化参数值，训练多个Lasso回归模型，并计算每个模型在测试数据上的得分。
最终，得分将保存在`scores`数组中，以便后续分析和比较模型性能。
"""
for i, a in enumerate(alphas):
    lasso = Lasso(random_state = seed)
    lasso.set_params(alpha = a)
    lasso.fit(X_train, y_train)
    scores[i] = lasso.score(X_test, y_test)
    
    
# 方法2.交叉验证 使用LassoCV（Lasso交叉验证）来选择最佳的正则化参数，并训练模型
"""
使用LassoCV模型自动选择效果最好的正则化参数，并用该参数训练模型。
然后计算模型在训练数据上的得分，并获取最佳正则化参数的值。这些信息可用于评估模型的性能以及模型参数的选择。
"""
# 创建一个LassoCV对象lassocv，其中cv参数指定了交叉验证的折数（这里是10折），
# random_state参数用于设置随机种子以保证结果的可重复性。
lassocv = LassoCV(cv = 10, random_state = seed)
# 使用特征集合features和目标变量target对LassoCV模型进行训练，它会自动进行交叉验证来选择最佳的正则化参数。
lassocv.fit(features, target)
# 使用训练好的LassoCV模型计算在训练数据上的得分，评估模型的性能。
lassocv_score = lassocv.score(features, target)
# 获取LassoCV模型选择的最佳正则化参数值，该值保存在alpha_属性中。
lassocv_alpha = lassocv.alpha_


# 画图找最高点对应的就是最好的参数值
plt.figure(figsize = (10, 4))
plt.plot(alphas, scores, '-ko')
plt.axhline(lassocv_score, color = c)
plt.xlabel(r'$\alpha$')
plt.ylabel('CV Score')
plt.xscale('log', basex = 2)
sns.despine(offset = 15)

print('CV results:', lassocv_score, lassocv_alpha)

# 使用LassoCV模型进行回归分析
# 创建了一个LassoCV模型对象model_l1，它是一个使用交叉验证（cv=10）来选择最佳正则化参数的Lasso回归模型。
# alphas是用于正则化参数的候选值序列，X_train是训练数据集的特征集合，y_train是对应的目标变量（即训练集中的预测目标）。
model_l1 = LassoCV(alphas = alphas, cv = 10, random_state = seed).fit(X_train, y_train)
# 使用训练好的model_l1模型对测试集的特征集合X_test进行预测，得到预测结果y_pred_l1，即模型对测试集样本的目标变量的预测值。
y_pred_l1 = model_l1.predict(X_test)
# 计算模型在测试集上的得分（即拟合优度）。X_test是测试集的特征集合，y_test是测试集中的实际目标变量。
# 模型的得分可以用来评估模型对新样本的预测准确性，得分越高表示模型的拟合能力越好。
model_l1.score(X_test, y_test)

#0.83307445226244159

def MSE(y_true,y_pred):
    mse = mean_squared_error(y_true, y_pred)
    print('MSE: %2.3f' % mse)
    return mse

def R2(y_true,y_pred):    
    r2 = r2_score(y_true, y_pred)
    print('R2: %2.3f' % r2)     
    return r2

MSE(y_test, y_pred_l1); R2(y_test, y_pred_l1);


#MSE: 3870543.789
#R2: 0.833