线性回归 Linear regression

分类classification/class/category
回归regression:predict a number from infinitely many possible numbers

回归：回归(Regression)是对连续型数据进预测，区别于分类（Classification）问题，天气预报预测明天是否下雨，是一个二分类问题；预测明天的降雨量多少，就是一个回归问题。
线性回归是利用回归方程(函数)对一个或多个自变量(特征值)和因变量(目标值)之间关系进行建模的一种分析方式。只有一个自变量的情况称为单变量回归，大于一个自变量情况的叫做多元回归。
线性回归(Linear regression)：是统计学中最基础的数学模型，几乎各个学科的研究中都能看到线性回归的影子，比如量化金融、计量经济学等。
线性回归的一种最直观解法是最小二乘法，其损失函数是误差的平方，具有最小值点，可以通过解矩阵方程求得这个最小值。

线性回归基本原理

线性回归模型 Multiple Variable Linear Regression

参阅：吴恩达，斯坦福CS231课程资料

单变量线性回归模型:

$$f_{w,b} (x) = wx +b $$

$$x \rightarrow f \rightarrow \widehat{y}$$

多元线性模型的方程:

$$f_{\mathbf{w},b}(\mathbf{x}) = w_0x_0 + w_1x_1 +... + w_{n-1}x_{n-1} + b \tag{1}$$

或者向量形式:

$$f_{\mathbf{w},b}(\mathbf{x}) = \mathbf{w} \cdot \mathbf{x} + b \tag{2}$$

这里$\cdot$表示向量点积

如何评价曲线拟合的效果呢？提出损失函数来评价拟合的效果，一般我们使用平方误差损失函数
平方误差损失函数:

$$J_{(w,b)} = \frac{1}{2m} \sum_{i=1}^{m}{(\widehat{y}_{(i)} - y_{(i)})^2}$$

$$\Downarrow$$

$$J_{(w,b)} = \frac{1}{2m} \sum_{i=1}^{m}{(f_{w,b} (x_{(i)})- y_{(i)})^2}$$

其中$\frac{1}{2}$是为了求导便利，$\frac{1}{m}$是为了将损失平均化，消除样本量m的影响。线性回归的损失函数用最小二乘法，等价于当预测值与真实值的误差满足正态分布时的极大似然估计，这里不再展开。下面我们考虑如何求最优的参数w和b
参数 :

$$w,b$$

目标 : 最小化损失函数

$$\underset{w,b}{minimize} J_{(w,b)}$$

多元方程的代价函数$J(\mathbf{w},b)$:

$$J(\mathbf{w},b) = \frac{1}{2m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})^2 \tag{3}$$

这里：

$$f_{\mathbf{w},b}(\mathbf{x}^{(i)}) = \mathbf{w} \cdot \mathbf{x}^{(i)} + b \tag{4}$$

正规方程求解 The Normal Equation

求解析解 Linear Regression, closed-form solution

解释求解过程图

python实现 LinearRegression

class sklearn.linear_model.LinearRegression(*, fit_intercept=True, normalize=False, copy_X=True, n_jobs=None)
LinearRegression使用系数w =（w1，…，wp）拟合线性模型，以最小化数据集中实际目标值与通过线性逼近预测的目标之间的残差平方和。

值得关注的参数：偏差项（ intercept_ ）与特征权重（coef_）

参数图

python


from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
mypath = 'D:\\mypath'    #指定数据存储路径
house_data = fetch_california_housing(data_home= mypath)

# 划分数据集
x_train, x_test, y_train, y_test = train_test_split(house_data.data, house_data.target, test_size=0.3, random_state=24)

# 需要做标准化处理对于特征值处理
std_x = StandardScaler()
x_train = std_x.fit_transform(x_train)
x_test = std_x.fit_transform(x_test)
# 对于目标值进行标准化
std_y = StandardScaler()
# 在新版的sklearn中，所有的数据都应该是二维矩阵，哪怕它只是单独一行或一列，所以需要使用.reshape(1,-1)进行转换
y_train = std_y.fit_transform(y_train.reshape(-1, 1))
y_test = std_y.transform(y_test.reshape(-1, 1))
y_test = std_y.inverse_transform(y_test)

# 使用线性模型进行预测
lr = LinearRegression()
lr.fit(x_train, y_train)

y_lr_predict = std_y.inverse_transform(lr.predict(x_test))

# print("预测的结果为：", y_lr_predict)
print("均方误差为：", mean_squared_error(y_test, y_lr_predict))

LinearRegression类基于 scipy.linalg.lstsq() 函数（名称代表“最小二乘”），可以直接调用它。它计算了矩阵的伪逆（pseudoinverse ），伪逆本身是使用是使用被称为奇异值分解（Singular ValueDecomposition，SVD）的标准矩阵分解技术来计算的，Scikit-Learn的LinearRegression类使用的SVD方法的复杂度约为O（n2）。如果特征数量加倍，那计算时间大约是原来的4倍。这种方法比计算标准方程更有效。如果矩阵XT X是不可逆的（即奇异的），标准方程可能没有解，例如 m < n 或某些特征是多余的，但伪逆总是有定义的。
特征数量比较大（例如100 000）时，标准方程和SVD的计算将极其缓慢。好的一面是，相对于训练集中的实例数量（O（m））来说，两个都是线性的，所以能够有效地处理大量的训练集，只要内存足够。同样，线性回归模型一经训练（不论是标准方程还是其他算法），预测就非常快速：因为计算复杂度相对于想要预测的实例数量和特征数量来说都是线性的。换句话说，对两倍的实例（或者是两倍的特征数）进行预测，大概需要两倍的时间。

梯度下降法 Gradient Descent(GD)

梯度下降是一种非常通用的优化算法，能够为大范围的问题找到最优解。梯度下降的中心思想就是迭代地调整参数从而使成本函数最小化。假设你迷失在山上的浓雾之中，你能感觉到的只有你脚下路面的坡度。快速到达山脚的一个策略就是沿着最陡的方向下坡。这就是梯度下降的做法：通过测量参数向量θ相关的误差函数的局部梯度，并不断沿着降低梯度的方向调整，直到梯度降为0，到达最小值！具体来说，首先使用一个随机的θ值（这被称为随机初始化），然后逐步改进，每次踏出一步，每一步都尝试降低一点成本函数（如MSE），直到算法收敛出一个最小值。

梯度下降法优化示意图：

示意图

在梯度下降的描述中，模型参数被随机初始化并反复调整使成本函数最小化。学习步长与成本函数的斜率成正比，因此，当参数接近最小值时，步长逐渐变小。
梯度下降计算公式 :$(\alpha ：学习率)$

$$ \begin{cases} w = w -\alpha \frac{\partial}{\partial w} J_{(w,b)} \\ b = b- \alpha \frac{\partial}{\partial b} J_{(w,b)} \end{cases} $$

由平方误差公式推导迭代公式(将上述公式分别代入求导即可)，得：

$$ \begin{cases} w = w - \alpha \cdot \frac{1}{m} \sum_{i=1}^{m}{(f_{w,b} (x_{(i)})- y_{(i)})} \cdot x_{(i)} \\ b = b - \alpha \cdot \frac{1}{m} \sum_{i=1}^{m}{(f_{w,b} (x_{(i)})- y_{(i)})} \end{cases} $$

幸好，线性回归模型的MSE（最小化均方误差）成本函数恰好是个凸函数，这意味着连接曲线上任意两点的线段永远不会跟曲线相交。也就是说，不存在局部最小值，只有一个全局最小值。它同时也是一个连续函数，所以斜率不会产生陡峭的变化。这两点保证的结论是：即便是乱走，梯度下降都可以趋近到全局最小值（只要等待时间足够，学习率也不是太高）。
应用梯度下降时，需要保证所有特征值的大小比例都差不多（比如使用Scikit-Learn的StandardScaler类），否则收敛的时间会长很多。训练模型也就是搜寻使成本函数（在训练集上）最小化的参数组合。这是模型参数空间层面上的搜索：模型的参数越多，这个空间的维度就越多，搜索就越难。同样是在干草堆里寻找一根针，在一个三百维的空间里就比在一个三维空间里要棘手得多。幸运的是，线性回归模型的成本函数是凸函数，针就躺在碗底。

对实数集上的函数，可通过求二阶导数来判别：若二阶导数在区间上非负，则称为凸函数；若二阶导数在区间上恒大于0，则称为严格凸函数。

批量梯度下降 Batch Gradient Descent(BGD)

批梯度下降法：每次梯度下降都使用了所有训练数据，基于完整的训练集。每次迭代时都选择整体最优方向，因此较容易得到全局最优解，但是在数据量较大时，计算会非常慢。

随机梯度下降 Stochastic gradient descent(SGD)

随机梯度下降，每一步在训练集中随机选择一个实例，并且仅基于该单个实例来计算梯度。它在一次迭代时只考虑一个训练样本。高效，容易实现，但需要许多超参数：比如正则项参数、迭代数。SGD对于特征标准化是敏感的。

由于算法的随机性质，SGD比批量梯度下降要不规则得多。成本函数将不再是缓缓降低直到抵达最小值，而是不断上上下下，但是从整体来看，还是在慢慢下降。随着时间的推移，最终会非常接近最小值，但是即使它到达了最小值，依旧还会持续反弹，永远不会停止。所以算法停下来的参数值肯定是足够好的，但不是最优的。
当成本函数非常不规则时，随机梯度下降其实可以帮助算法跳出局部最小值，所以相比批量梯度下降，它对找到全局最小值更有优势。随机性的好处在于可以逃离局部最优，但缺点是永远定位不出最小值。要解决这个困境，有一个办法是逐步降低学习率。开始的步长比较大（这有助于快速进展和逃离局部最小值），然后越来越小，让算法尽量靠近全局最小值。这个过程叫作模拟退火，因为它类似于冶金时熔化的金属慢慢冷却的退火过程。确定每个迭代学习率的函数叫作学习率调度。如果学习率降得太快，可能会陷入局部最小值，甚至是停留在走向最小值的半途中。如果学习率降得太慢，你需要太长时间才能跳到差不多最小值附近，如果提早结束训练，可能只得到一个次优的解决方案。

要使用带有Scikit-Learn的随机梯度下降执行线性回归，可以使用SGDRegressor类，该类默认优化平方误差成本函数。以下代码最多可运行1000个轮次，或者直到一个轮次期间损失下降小于0.001为止（max_iter=1000，tol=1e-3）。它使用默认的学习调度以0.1（eta0=0.1）的学习率开始。最后，它不使用任何正则化（penalty=None）

python


from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=1000, tol=1e-3, penalty=None, eta0=0.1)
sgd_reg.fit(X, y.ravel())

参阅文档

WARNING
When using stochastic gradient descent, the training instances must be independent and identically distributed (IID) to ensure that the parameters get pulled toward the global optimum, on average. A simple way to ensure this is to shuffle the instances during training (e.g., pick each instance randomly, or shuffle the training set at the beginning of each epoch). If you do not shuffle the instances—for example, if the instances are sorted by label—then SGD will start by optimizing for one label, then the next, and so on, and it will not settle close to the global minimum.
警告
当使用随机梯度下降时，训练实例必须是独立和同分布的（IID），以确保参数平均被拉向全局最优。确保这一点的一个简单方法是在训练过程中打乱实例（例如，随机挑选每个实例，或在每个历元开始时打乱训练集）。如果你不打乱实例——例如，如果实例按标签排序——那么SGD将从优化一个标签开始，然后优化下一个标签，以此类推，它不会稳定在接近全局最小值。

随机平均梯度法 Stochasitc Average Gradient(SAG)

由于收敛的速度太慢，有人提出SAG等基于梯度下降的算法。

小批量梯度下降 Mini-Batch Gradient Descent

from sklearn.linear_model import SGDClassifier
每次迭代时，随机从训练集中选取一批样本进行计算，综合了随机梯度下降和批量梯度下降两种方法。通常将选取的样本数量称为batch。在时间和硬件配置允许的前提下，优先选择batch_size较大的。

对比图

Scikit-learn has a gradient descent regression model sklearn.linear_model.SGDRegressor. Like your previous implementation of gradient descent, this model performs best with normalized inputs. sklearn.preprocessing.StandardScaler will perform z-score normalization as in a previous lab. Here it is referred to as 'standard score'.

多元变量线性回归的梯度下降

Gradient descent for multiple variables:

$$\begin{align*} \text{repeat}&\text{ until convergence:} \; \lbrace \newline\; & w_j = w_j - \alpha \frac{\partial J(\mathbf{w},b)}{\partial w_j} \tag{5} \; & \text{for j = 0..n-1}\newline &b\ \ = b - \alpha \frac{\partial J(\mathbf{w},b)}{\partial b} \newline \rbrace \end{align*}$$

where, $n$ is the number of features, parameters $w_j$ , $b$ , are updated simultaneously and where

$$ \begin{align} \frac{\partial J(\mathbf{w},b)}{\partial w_j} &= \frac{1}{m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})x_{j}^{(i)} \tag{6} \\ \frac{\partial J(\mathbf{w},b)}{\partial b} &= \frac{1}{m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)}) \tag{7} \end{align} $$

m 是数据集中训练样本数
$f_{\mathbf{w},b}(\mathbf{x}^{(i)})$ is the model's prediction, while $y^{(i)}$ is the target value

mutilplefeartures

梯度下降法是受推荐的求最优w、b的方法

多元特征的向量化

程序示例

房价预测程序示例（多元线性回归，sklearn库）：

python


from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
mypath = 'D:\\mypath'    #指定数据存储路径
house_data = fetch_california_housing(data_home= mypath)

X_train = house_data.data
y_train = house_data.target

scaler = StandardScaler()
X_norm = scaler.fit_transform(X_train)
print(f"Peak to Peak range by column in Raw        X:{np.ptp(X_train,axis=0)}")   
print(f"Peak to Peak range by column in Normalized X:{np.ptp(X_norm,axis=0)}")
sgdr = SGDRegressor(max_iter=1000)
sgdr.fit(X_norm, y_train)
print(sgdr)
print(f"number of iterations completed: {sgdr.n_iter_}, number of weight updates: {sgdr.t_}")
b_norm = sgdr.intercept_
w_norm = sgdr.coef_
print(f"model parameters:                   w: {w_norm}, b:{b_norm}")

# make a prediction using sgdr.predict()
y_pred_sgd = sgdr.predict(X_norm)
# make a prediction using w,b. 
y_pred = np.dot(X_norm, w_norm) + b_norm  
print(f"prediction using np.dot() and sgdr.predict match: {(y_pred == y_pred_sgd).all()}")
print(f"Prediction on training set:\n{y_pred[:8]}" )
print(f"Target values \n{y_train[:8]}")


# 查看拟合效果，可以自行修改参数
import matplotlib.pyplot as plt
# X_features = ['size(sqft)','bedrooms','floors','age']
fig,ax=plt.subplots(1,2,sharey=True)
for i in range(len(ax)):
    ax[i].scatter(X_train[:,i],y_train, label = 'target')
    ax[i].set_xlabel(house_data.feature_names[i])
    ax[i].scatter(X_train[:,i],y_pred,color = "r", label = 'predict')
ax[0].set_ylabel("Price"); ax[0].legend();
fig.suptitle("target versus prediction using z-score normalized model")
plt.show()

过拟合与欠拟合

过拟合 (overfit) ：一个假设在训练数据上能够获得比其他假设更好的拟合，但是在测试数据集上却不能很好地拟合数据，此时认为这个假设出现了过拟合的现象。(模型过于复杂)
欠拟合 (underfit) ：一个假设在训练数据上不能获得更好的拟合，并且在测试数据集上也不能很好地拟合数据，此时认为这个假设出现了欠拟合的现象。(模型过于简单)

欠拟合的解决

（1）添加特征。欠拟合是因为特征数不够，可以通过组合、泛化、相关性等得到新特征；

（2）添加多项式特征。在机器学习算法普遍使用该方法，如将线性模型通过添加二次项或者三次项使模型泛化能力更强；

（3）减少正则化参数。正则化的目的是用来防止过拟合的，但由于模型出现了欠拟合，则需要减少正则化参数。

过拟合的解决

（1）收集更多的训练数据。collect more training examples

（2）选择特征，剔除部分不必要特征。select feztures

（3）采用正则化方法，regularization。正则化方法包括L0正则、L1正则和L2正则。

（4）采用dropout方法。这个方法在神经网络里面很常用。dropout方法是ImageNet中提出的一种方法，通俗一点讲就是dropout方法在训练的时候让神经元以一定的概率不工作。

正则化

L0范数是指向量中非0的元素的个数。那么，在规则化权重矩阵w时，则希望非0元素尽量少，甚至不存在，起到稀疏化矩阵的作用。这比较极端，也是L0较少使用的一个原因。
L1范数是指向量中各个元素绝对值之和，可以使得权重矩阵w中的元素等于0，也称为“稀疏规则算子”。虽然L0范数和L1范数都可以实现稀疏性，但是因为L0范数不连续，很难优化求解（NP问题），而且L1范数是L0范数的最优凸近似，所以L1范数常被使用。

l1图

我们将L1范数称为loss函数的正则项，也是L1正则化，加入L1范数的线性回归称作LASSO回归。L1正则化会使不重要的特征的权重置为0，L1范数是稀疏约束。
L1范数使得参数稀疏主要作用有两方面：
（1）特征选择：实现特征的自动选择。实际问题中，特征都很多，且存在一些不相关或者包含不重要信息的特征。如果将所有特征都加入模型，则会使得模型求解速度慢。使用L1范数可以进行特征自动选择，主要是将与y不相关或含不重要信息的特征的权重等于0，可以达到去除这些特征，保留重要特征的效果。
（2）可解释性：去除不重要特征，留下少量的重要特征，解释性更强。
为了解决L1没有保留局部特征的问题，引入L2范数替代L1范数。L2范数使得权重矩阵w中的元素趋向于0，越小使得模型的复杂度降低。

L2范数是指向量各元素的平方和然后求平方根。可以使得W的每个元素都很小，都接近于0，但与L1范数不同，它不会让它等于0，而是接近于0。L2正则项有防止过拟合的作用，原因是它使得参数w变小加剧，且更小的参数值w意味着模型的复杂度更低，对训练数据的拟合刚刚好（奥卡姆剃刀），不会过分拟合训练数据，从而使得不会过拟合，以提高模型的泛化能力。Ridge回归

l2图

下面再看看吴恩达的课程中介绍的正则化，首先是在代价函数中引入正则化项：

angrengreg.png

and2.png

然后看看正则化线性回归：

and3

and4

拓展

线性回归的变形有很多，这里选择几个介绍一下。

岭回归

岭回归，Ridge，也称为Tikhonov正则化。岭回归成本函数：

$$ J(\theta) = MSE(\theta) + \alpha\frac{1}{2} \sum_{i=1}^n \theta_i^2 $$

公式中偏置项没有正则化，可以参考上面的吴恩达的课件。在执行岭回归之前缩放数据（例如使用StandardScaler）很重要，因为它对输入特征的缩放敏感。大多数正则化模型都需要如此。 class sklearn.linear_model.Ridge(alpha=1.0, *, fit_intercept=True, copy_X=True, max_iter=None, tol=0.0001, solver='auto', positive=False, random_state=None)

Ridge方法相当于SGDRegressor(penalty='l2', loss="squared_loss"),只不过SGDRegressor实现了一个普通的随机梯度下降学习，推荐使用Ridge(实现了SAG)
关注参数：alpha:正则化力度，也叫 λ
必须为正浮点数。正则化改善了问题的状况，并减少了估计的方差。值越大表示正则化更强。

solver {‘auto’, ‘svd’, ‘cholesky’, ‘lsqr’, ‘sparse_cg’, ‘sag’, ‘saga’}, default=’auto’

‘auto’根据数据类型自动选择求解器。
‘svd’使用X的奇异值分解来计算Ridge系数。奇异矩阵比“ Cholesky”更稳定。
'cholesky'使用标准的scipy.linalg.solve函数通过点（XT，X）的Cholesky分解获得封闭的解
'sparse_cg'使用scipy.sparse.linalg.cg中的共轭梯度求解器。作为一种迭代算法，此求解器比' cholesky '更适合用于大规模数据（设置tol和max_iter的可能性）。
‘lsqr’使用专用的正则化最小二乘例程scipy.sparse.linalg.lsqr。它是最快的，并且使用了迭代过程。
‘sag’使用随机平均梯度下降，“ saga”使用其改进的无偏版本SAGA。两种方法都使用迭代过程，并且当n_samples和n_features都很大时，通常比其他求解器更快。注意，只有在比例大致相同的特征上才能保证“ sag”和“ saga”快速收敛。您可以使用sklearn.preprocessing中的缩放器对数据进行预处理。

python


from sklearn import linear_model
reg = linear_model.Ridge(alpha=.5)
reg.fit([[0, 0], [0, 0], [1, 1]], [0, .1, 1])
reg.coef_
reg.intercept_

Lasso回归

最小绝对收缩和选择算子回归，使用L1范数。

$$ J(\theta) = MSE(\theta) + \alpha \sum_{i=1}^n \vert \theta_i\vert $$

Lasso回归的一个重要特点是它倾向于完全消除掉最不重要特征的权重（也就是将它们设置为零）

python


from sklearn import linear_model
reg = linear_model.Lasso(alpha=0.1)
reg.fit([[0, 0], [1, 1]], [0, 1])
reg.predict([[1, 1]])

弹性网络 Elastic Net

弹性网络是介于岭回归和Lasso回归之间的中间地带。正则项是岭和Lasso正则项的简单混合，你可以控制混合比r。当r=0时，弹性网络等效于岭回归，而当r=1时，弹性网络等效于Lasso回归。弹性网络成本函数：

$$ J(\theta) = MSE(\theta) + r\alpha\sum_{i=1}^n \vert \theta_i\vert + \frac{1-r}{2}\alpha\sum_{i=1}^n \theta_i^2 $$

那么什么时候应该使用普通的线性回归（即不进行任何正则化）、岭、Lasso或弹性网络呢？通常来说，有正则化——哪怕很小，总比没有更可取一些。所以大多数情况下，你应该避免使用纯线性回归。岭回归是个不错的默认选择，但是如果你觉得实际用到的特征只有少数几个，那就应该更倾向于Lasso回归或是弹性网络，因为它们会将无用特征的权重降为零。一般而言，弹性网络优于Lasso回归，因为当特征数量超过训练实例数量，又或者是几个特征强相关时，Lasso回归的表现可能非常不稳定。

局部加权线性回归

局部加权线性回归（Locally Weighted Linear Regression，LWLR）。标准的线性回归是一种无偏差估计，在计算所有点的时候都是无偏差的计算误差并通过优化方法优化误差，如果针对不同的点能够对误差进行调整便可以一定程度上避免标准线性回归带来的欠拟合现象。当k取的值比较大时，LWLR就相当于无偏差的线性回归；当k取的值比较小时，LWLR拟合效果比较好，但是容易出现过拟合；

lwlr

python


#coding = utf-8
import numpy as np
import scipy.stats as stats
from math import *
import matplotlib.pyplot as plt


def getw(x0,x,k):
    """
    然后将这个距离的平方除以-2*k**2，并将这个结果作为矩阵w中对应位置的元素值。这里看起来像是在计算高斯核函数的一部分。
    """
    #w = np.matrix(np.zeros((m, m)))
    w = np.zeros([m,m])
    for i in range(m):
        #计算了x0与x中每个数据点的欧几里得距离（使用了numpy库的linalg.norm函数来计算两点间的距离）
        w[i, i] = exp((np.linalg.norm(x0 - x[i])) / (-2 * k ** 2))
    return w


def getyvalue(x1,x,y,k):
    """
    其中x1和x的维度应该是m，y的维度应该是1，而k是一个标量。
    在函数内部，首先创建了一个m维的零向量y_value和m x m的零矩阵w。这个矩阵的形状由常数m确定，而初始时所有的元素都是0。
    然后，函数通过一个for循环对y_value中的每一个元素进行赋值。在这个循环中，对于每一个i，函数getw(x[i],x, k)都会计算出一个w矩阵。
    然后使用这个w矩阵和x1、x、y进行一些矩阵运算，来计算出theta值。
    具体来说，这个代码在进行一种线性回归的计算，其中theta是线性回归的参数。
    """
    y_value = np.zeros(m)
    #w = np.matrix(np.zeros((m, m)))
    w = np.zeros([m,m])

    for i in range(m):
        w = getw(x[i],x, k)
        theta = np.linalg.inv(x1.T.dot(w).dot(x1)).dot(x1.T).dot(w).dot(y)
        y_value[i] = theta[0] + theta[1] * x[i]
    return y_value

def LR(x,y):
    """
    'LR'函数，用于执行线性回归操作。
    这个函数接收两个参数，x和y，其中x是一个二维数组或类似的数据结构（特征矩阵），y是一个一维数组或类似的数据结构（目标变量）。
    在函数内部，首先从sklearn库中导入LinearRegression类，然后创建一个LinearRegression对象lr。
    调用lr的fit方法，将x和y作为参数传入，这个方法用于拟合模型，根据输入的x和y数据来学习线性回归模型。
    接下来，通过计算lr.intercept_（截距）和lr.coef_（斜率）来预测新的x值对应的y值，其中lr.coef_是回归系数向量，lr.intercept_是截距。计算完成后，将结果存储在变量y1中。
    最后，打印出斜率（lr.coef_）和截距（lr.intercept_），并以元组的形式返回这两个值。
    这段代码的主要目的是根据给定的x和y数据，通过线性回归模型来预测新的x值对应的y值。
    作为对比，在图中用黄色yellow线绘制
    """
    from sklearn.linear_model import LinearRegression
    lr = LinearRegression()
    lr.fit(x, y)
    y1 = lr.intercept_ + x * lr.coef_
    print(lr.coef_, lr.intercept_)
    return  lr.intercept_, lr.coef_


if __name__ == "__main__":
    # 生成数据
    x = np.arange(1, 101)
    x = np.array([float(i) for i in x])
    y = x + [10 * sin(0.3 * i) for i in x] + stats.norm.rvs(size=100, loc=0, scale=1.5)
    #plt.figure(figsize=(12, 6))
    #plt.scatter(x, y)
    #plt.show()

    x = x.reshape(-1, 1)
    x1 = np.c_[np.ones((100, 1)), x]
    y = y.reshape(-1, 1)
    m = len(x)

    y_lwlr = np.zeros(m)
    y_lwlr = getyvalue(x1,x,y,0.2)
    a,b = LR(x,y)
    y_lr = a + b*x
    plt.figure(figsize=(12,6))
    plt.scatter(x, y)
    plt.plot(x, y_lwlr, 'r')
    plt.plot(x, y_lr, 'y')
    plt.show()

程序输出为：

lroutput

“回归”一词的来历

史海拾贝。

史海拾贝

参阅：

斯坦福CS231，Andrew NG
《机器学习实战》Peter Harrington，[美]，ISBN: 9787115317957
《机器学习》西瓜书，周志华，南京大学，ISBN：9787302423287
《跟着迪哥学Python数据分析与机器学习实战》，唐宇迪，人民邮电出版社，ISBN: 9787115512444
https://scikit-learn.org/stable/modules/linear_model.html
http://scikit-learn.org.cn/view/4.html
经典算法（一）：线性回归 - Datamoon的文章 - 知乎
Locally Weighted Linear Regression - jason的文章 - 知乎