用于机器学习的线性代数:求解线性方程组
代数是机器学习算法的底线机制
什么是线性代数，为什么要关心？
先说数据科学中常见的图:散点图
从图表获得
上面的图显示了树的直径和高度之间的关系。每个点都是一棵树的样本。我们的任务是找到最佳拟合线来预测提供直径的高度。
我们如何做到这一点？这时候就需要线性代数了。
来自教育学院
线性回归是线性方程组的一个例子。线性代数是关于线性方程组的工作。我们开始使用矩阵和向量，而不是标量。
线性代数是理解你在机器学习中需要的微积分和统计学的关键。如果你能在向量和矩阵的层面上理解机器学习方法，你将提高对它们如何以及何时工作的直觉。更好的线性代数将全面提升你的游戏。
而理解线性代数的最好方法是什么？执行它。解线性方程组有两种方法:直接法和迭代法。在本文中，我们将使用直接方法，特别是高斯方法。
如何解一个线性方程组？
因为大多数情况下，我们会处理具有许多特征(或变量)的数据。我们将通过处理三维数据来使我们的线性方程组更通用。
让我们为上面的图生成一个示例:
其中，系数 x_0、x_1 和 x_2 以及相应的值 8、4、5 是图中点的样本。这些方程可以分解成矩阵 A、x 和 b
其中 A 和 b 是已知常数的矩阵，x 是未知变量的向量。
A = np.array([[2, 1, 5],
[4, 4, -4],
[1, 3, 1]])
b= np.array([8,4,5])
连接矩阵 A 和 b 得到
n = A.shape[0]
C=np.c_[A,b.reshape(-1,1)]
现在，我们准备用两个步骤来解决我们的问题:
- 应用高斯消去法将上述矩阵化为三角矩阵
这可以用下面的等式来表示:
2.应用向后替换来获得结果
让我们从第一步开始
部分旋转高斯消去法
要获得该矩阵:
这个想法很简单:
- 我们从第一行第一列的透视值开始:行=0，列= 0
- 求主元列的最大绝对值。如果该列中的所有值都是 0，我们就停止。
- 否则，我们交换 E0 和 E1
接下来，应用等效转换，通过以下方式将透视下的所有条目转换为 0:
- 求元素 j，I 和 pivot i，I 的比值(即 2/4 = 1/2)。
- 将 E0 中的所有元素乘以 1/2。将第 1 行中的所有元素减去 1/2 E0(即 2-(4 * 1/2)= 2–2 = 0)
#row
for j in range(i+1, n):
c = C[j,i]/C[i,i]
C[j,:] = C[j,:] - c*C[i,:]
每列简化后，我们继续右边的下一列。
重复该过程:
透视:行= 1，列= 1。最大绝对值:第 2 行中的 2。然后，置换 E1 和 E2
应用等效变换将透视下的所有条目转换为 0
把所有东西放在一起
不错！现在我们有了一个方程组:
一旦我们到了这里，这个方程组就很容易用回代法求解了
反向置换
从高斯消去法，我们得到一个三角矩阵
这个想法是通过自底向上求解来求解上面的方程组。使用从上一个等式获得的值来查找其他值
从第三排开始。除以 8 得到 x_3 的值。
X[n-1] = T[n-1,n]/T[n-1,n-1]
现在在第二行，我们有:
x_2 可以很容易地通过下式求解
用 x1 重复
因此，一般来说，反替换可以表示为:
获得结果
厉害！正如我们预测的那样，我们得到了一个解。为了确保在处理更大的矩阵时这是正确的，我们可以使用 NumPy 中的内置函数
>>> np.linalg.solve(A,b)array([1., 1., 1.])
我们得到的是解的向量，其中每个元素对应于 x_0，x_1，x_2
结论
祝贺你走到这一步！希望这篇文章能帮助你理解什么是线性代数，以及解线性方程组的机制之一。我尽量让这篇文章容易理解。但是我知道，如果你不熟悉线性代数，这可能是一个挑战。没关系。一步一步来。你接触线性代数越多，你就越能理解它。
你可以在我的 Github 中玩和试验上面的代码。
我喜欢写一些基本的数据科学概念，并尝试不同的算法和数据科学工具。你可以通过 LinkedIn 和 Twitter 与我联系。
如果你想查看我写的所有文章的代码，请点击这里。在 Medium 上关注我，了解我的最新数据科学文章，例如:
如果您一直在为矩阵运算导入 Numpy，但不知道该模块是如何构建的，本文将展示…
towardsdatascience.com](/how-to-build-a-matrix-module-from-scratch-a4f35ec28b56) [## 高效 Python 代码的计时
如何比较列表、集合和其他方法的性能
towardsdatascience.com](/timing-the-performance-to-choose-the-right-python-object-for-your-data-science-project-670db6f11b8e) [## 当生活不给你喘息的机会，如何学习数据科学
我努力为数据科学贡献时间。但是发现新的策略使我能够提高我的学习速度和…
towardsdatascience.com](/how-to-learn-data-science-when-life-does-not-give-you-a-break-a26a6ea328fd) [## 字典作为 If-Else 的替代
使用字典创建一个更清晰的 If-Else 函数代码
towardsdatascience.com](/dictionary-as-an-alternative-to-if-else-76fe57a1e4af) [## 字典作为 If-Else 的替代
使用字典创建一个更清晰的 If-Else 函数代码
towardsdatascience.com](/dictionary-as-an-alternative-to-if-else-76fe57a1e4af)
数据科学的线性代数向量
线性代数 I —数据科学中的向量介绍第一部分
随着 AI/ML 的民主化和像 Keras、scikit-learn 等开源库的出现，任何具有基本 python 知识的人都可以在不到 5 分钟的时间内建立一个工作的 ML 分类器。虽然这对于开始来说已经足够了，但如果你想了解不同的最大似然算法是如何工作的，或者将最新的 SOTA(最先进的)论文应用到你的特定领域，缺乏数学专业知识很快就会成为一个瓶颈，正如我亲身经历的那样。
在这组文章中，我将尝试为非数学读者一次介绍一个基本的数学概念，并展示它在人工智能领域的实际应用。
我们从最简单的向量开始。
矢量只是有方向的量。向量一些相关的真实世界的例子是力、速度、位移等。
要移动购物车，你需要向你想要移动购物车的方向推(施加力)。你在推车时所用的力可以用两个值来充分描述，推的强度(幅度)和你推车的方向。任何需要大小和方向来完整描述的量叫做矢量。
向量通常用粗体小写字符表示，如 v，w 等。由于用纸笔书写粗体字比较困难，所以在用纸笔书写时，小写字的上方也会出现一个箭头。对于本文，我们将坚持使用粗体表示。
在图形上，矢量用箭头表示，箭头的长度表示矢量的大小(强度),箭头的角度(相对于参照系；在这种情况下是水平的)表示向量的方向，如下所示。
作者图片
请注意，并不要求向量应该从原点(0，0)开始。他们可以从任何一点开始。例如上图中的 u = w 和 v = a ，因为它们大小相同，方向相同。
用数学方法表示矢量有很多种方法。作为数据科学家，我们感兴趣的是将它们表示为一组数字。因此，向量 u 可以表示为(2，2 ),而向量 v 可以表示为(4，1)。矢量 w 和 a 也是如此。
虽然我们很容易在二维和三维空间中看到向量，但是向量的概念并不局限于二维和三维空间。它可以推广到任何数量的维度，这就是向量在机器学习中如此有用的原因。
例如， c = (2，1，0)表示三维空间中的向量，而 d = (2，1，3，4)表示四维空间中的向量。作为人类，虽然我们无法想象高于 3 的维度，但表示向量的数学方式让我们有能力在更高维度的向量空间上执行操作。
到现在为止，你一定很无聊，想知道为什么作为一个 ML 爱好者，你需要学习基础物理和向量。事实证明，向量在机器学习中有多种应用，从构建推荐引擎到自然语言处理的单词的数字表示等，并形成了 NLP 的所有深度学习模型的基础。
让我们从 numpy 和 tensorlfow 中如何实现 vectors 的代码示例开始。
请注意:完整的代码可作为最后一节的要点。为便于说明，插入了相关小节的图片。
作者图片
下面你可以看到一个简单的 word2vec 实现来展示向量在自然语言处理中的实际应用
作者图片
如你所见，以高维向量格式存储单词是向量在自然语言处理中的主要应用之一。这种类型的嵌入保留了单词的上下文。
在下一节中，我们将学习基本的运算，如加法和减法，以及它如何应用于向量。
向量加法
现在我们已经定义了什么是向量，让我们看看如何对它们进行基本的算术运算。
让我们用同样的两个向量 u 和 v，对它们进行向量加法。
作者图片
为了以图形方式添加两个向量 u 和 v ，我们移动向量 v ，使其尾部从向量 u 的头部开始，如上所示(DE 和 EF 行)。两个矢量之和就是始于 u 尾部，止于 v 头部的矢量 b (直线 DF)。
为了更好的直觉，让我们举一个开车去杂货店的真实例子。在路上，你在加油站停下来加油。让我们假设向量 u 代表加油站(点 E)离你家(点 D)有多远。如果矢量 u 代表从加油站到杂货店的距离(位移)，那么从 u 的尾部画到 v 的头部的矢量 b 代表总和 u + v. It 代表杂货店(F 点)离你家(初始起点 D)有多远。
数学上 b 可以通过看图表示为(6，3)。
还有其他计算矢量加法的图形方法，如平行四边形法，你可以自己探索。
现在不可能每次我们想做向量算术的时候都画一个图，特别是当涉及到高维向量的时候。幸运的是，向量的数学表示为我们做向量加法提供了一种简单的方法。
由于每个向量都是一组数字，让我们看看如果我们将每个向量的相应数字相加，会得到什么。
在上面的例子中，
=(2，2)
v = (4，1)
b = (2+4，2+1) = (6，3)，这与图解得到的解相同。
因此，矢量相加可以通过简单地将每个矢量的相应元素相加来完成，正如你可能已经推断出的那样，只有具有相同维数的矢量才能相加在一起。让我们看一个代码实现
作者图片
向量减法
在我们继续学习向量减法之前，让我们快速看一下向量的另一个有用的性质——标量乘法。
标量只不过是一个只有大小没有方向的量。例如，任何整数都是标量。标量的一个真实世界的例子是质量(体重)，一个人的身高等。
让我们看看当我们把一个矢量乘以一个标量时会发生什么。
u = (2，2)
v=(4，1)
如果我们想用标量 C = 3 乘以 u，一种直观的方法是将向量 u (2，2)中的单个数字乘以 3。让我们看看那看起来怎么样
d= C xu =3 xu=(3 x 2，3 x 2) = (6，6)
让我们在图上画出向量。
作者图片
如你所见，将向量 u 乘以一个正标量值会产生一个方向相同的新向量 d ,但是其大小按系数 C = 3 缩放
让我们试着用负值 C = -1 乘以一个向量
e= C xv=-1 xv=(-1 x 4，-1 x 1) = (-4，-1)。
让我们画出图来，看看是什么样子。
作者图片
如您所见，将向量 v 乘以-1 会产生一个大小相同但方向相反的向量，可表示为
e = -v 或 e + v = 0 (空矢量)
向量减法在图形上可以被认为是向量加法的特殊情况，其中 u -v = u + -v
图解求解
作者图片
****w = u
a =-v
b = w+a = u±v = u-v =(-2，1)
从图中可以看出， b = c 或矢量减法 u -v 等于从 v 的头部到 u 的头部(头部之间的距离)所画出的矢量 c
直观上，这是有道理的，也是符合传统数制的。
7 -5 = 2(其中 2 是加 5 等于 7 的数量)
5 + 2 = 7
类似地，如果你看这个图， c 是一个矢量，当它加到 v 上时，得到矢量 v
u = v + c
现在让我们用数学方法来做，在两个向量中减去各个分量。
c=u-v=(2，2) -(4，1) = (2 -4，2–1)=(-2，1)
其结果与图解法相同。下面给出一个代码示例供参考。
作者图片
在下一节中，让我们看一个代码示例，看看如何在单词向量中使用向量算法来查找单词之间的关系。
由于这是我的第一篇关于媒体的文章，我期待着反馈和建议。请在评论区告诉我如何更好地展示这些内容，如果你觉得有用的话。
完整的笔记本可以在我的 GitHub 库下找到
github.com](https://github.com/kishore145/Mathematics)**
参考
[1] JonKrohn，LML 基础课程-线性代数(2020)，O’Riley
带有 XTensor 的 C++中的线性代数就像 Numpy
使用张量库快速浏览 C++中的代数运算
(图片由作者提供)
介绍
对科学计算和数据科学非常重要的一个概念是在给定的编程语言中操作、创建和利用矩阵作为类型的能力。对于 Python，数据科学家使用 Numpy 库，它提供了一种奇妙且易于理解的方法来在 Python 内部创建和操作矩阵。对于 Julia 来说，矩阵和向量的概念从语言一开始就包含在内。对于 C++，有一个很棒的新兴库叫做 XTensor。
XTensor 是一个 C++库，它允许程序员使用与 Numpy 语法相似的矩阵和数组。与 Numpy 的相似之处是 XTensor 的真正优势，因为两者之间的转换可以相对容易地完成。为了使用 XTensor，您当然需要安装它。幸运的是，对于数据科学家来说，这可以在您的环境中通过 Conda 包管理器轻松完成
conda install -c conda-forge xtensor
或者，如果您使用的是 Debian 或基于 Debian 的 Linux 发行版，可以通过 Apt 软件包管理器获得 XTensor:
sudo apt-get install xtensor-dev
最后但同样重要的是，如果您在 FreeBSD 或另一个类似 Unix 的操作系统上，并且有适当的基本构建头文件，您总是可以从源代码构建:
cmake -DCMAKE_INSTALL_PREFIX=path_to_prefix ..
make install
或者，对于 Windows 用户:
mkdir build
cd build
cmake -G "NMake Makefiles" -DCMAKE_INSTALL_PREFIX=path_to_prefix ..
nmake
nmake install
基本矩阵
如果您是 Python 和 Numpy 用户，那么您很幸运使用了 XTensor。这里有一个绝对无价的 XTensor 数字备忘单的链接，你可以用它来应用你的数字熟悉度:
[## 从 numpy 到 xtensor - xtensor 文档
惰性辅助函数返回张量表达式。返回类型不包含任何值，而是在访问或…
xtensor.readthedocs.io](https://xtensor.readthedocs.io/en/latest/numpy.html)
对于基于 Python 的数据科学家和最终用户来说，这将使从 Numpy 到 XTensor 的过渡变得极其简单和流畅。安装 XTensor 后，您可以编译并执行以下代码，以确保头文件现在就在您的机器上:
#include <iostream>
// Base Arrays:
#include "xtensor/xarray.hpp"
// XTensor IO:
#include "xtensor/xio.hpp"
// XTensor View:
#include "xtensor/xview.hpp"
int main()
{
return 0;
}
如果编译和运行该代码时什么都没有发生，那么
恭喜你！
您的安装成功，XTensor 现在在您的机器上可用。如果输出是一个错误，您很可能需要尝试不同的安装方法，或者从源代码构建。当然，需要确保的另一件重要事情是，您处于安装了 XTensor 的 Conda 环境中。
接下来，我们当然会通过添加我们在那个小测试中添加的内容来开始我们的新代码。
#include <iostream>
// Base Arrays:
#include "xtensor/xarray.hpp"
// XTensor IO:
#include "xtensor/xio.hpp"
// XTensor View:
#include "xtensor/xview.hpp"
包含这些头文件将给我们一个新的名称空间 xt。首先，为了创建一个矩阵，我们将使用 xt 名称空间中的 xarray 对象模板。对于这个例子，我们将有一个充满整数的矩阵，所以每当我们创建新对象时，我们将提供“int”参数。我把它命名为 array，我们会把一些整数放进去。
xt::xarray<int>array
({{500, 500, 300}, {500, 600, 800}})
有两种不同类型的 x 传感器矩阵:
- xarray 设置了动态维的矩阵。
- x tensor——具有静态维度集的矩阵。
接下来，如果我们在笔记本中，那么我们可以简单地输入“array”来显示我们的新矩阵。
(图片由作者提供)
或者，如果您正在编译和执行代码，那么您将需要使用 cout 并将这些代码放入 main()函数中。在本例中，我不打算这样做，因为我是在笔记本中工作，但下面是您应该如何构建此代码，例如:
#include <iostream>
// Base Arrays:
#include "xtensor/xarray.hpp"
// XTensor IO:
#include "xtensor/xio.hpp"
// XTensor View:
#include "xtensor/xview.hpp"
int main()
{
xt::xarray<int>array
({{500, 500, 300}, {500, 600, 800}});
std::cout << array;
}
在某些情况下，你也可以这样做:
std::cout << array;
这就是包含 iostream off the bat 将派上用场的地方，即使您正在笔记本上工作。
我们可以使用 shape()函数重塑新数组，它是 array 对象的子对象。
array.reshape({3, 2})
(图片由作者提供)
我们也可以将一个新类型投射到数组的所有 dim 中。这是很有用的，例如，如果我们想要计算小数部分，但最初创建了一个充满整数的矩阵。我们可以使用 cast 函数来实现这一点，该函数在 xt 名称空间中。
xt::cast<double>(array)
(图片由作者提供)
我们还可以使用所有典型的 Numpy 初始化器/生成器来创建新矩阵。机器学习中常用的一个非常有用的版本是 np.zeroes()。在 xtensor 中，相当于 XT 名称空间中的零，这当然是一个从我们的常规 xarray 继承的对象。注意，它是从 xarray 而不是 xtensor 继承的，所以我们新发现的矩阵的维度将是可变的。
xt::zeros<double>({3, 4})
(图片由作者提供)
这也适用于 Numpy 中的其他生成器，比如 eye、ones 和其他许多生成器。
(图片由作者提供)
当然，我们也可以像使用 Numpy 一样索引数组。一个微小的区别是圆括号用于调用索引，而不是方括号。
(图片由作者提供)
操作
虽然能够访问生成器和创建矩阵是件好事，但是 XTensor 当然不仅仅是等同于 Numpy。我们基本上可以完成 Numpy 中任何可用的操作，比如将两个矩阵连接在一起:
// Our second array
xt::xarray<int>arraytwo
({{500, 500, 300}, {500, 600, 800}})// Reshape it to fit first array's dimensions
arraytwo.reshape({3, 2})// Concatenate it to the first array
xt::concatenate(xtuple(array, arraytwo))
(图片由作者提供)
我们也可以挤压，例如:
xt::squeeze(array)
我们还可以访问您可能需要对矩阵执行的标准运算，例如求和:
xt::sum(array, {0, 1})
// ^^ Array ^^^ Axis
平均:
xt::mean(array)
和标准:
xt::stddev(array, {0,1})
(图片由作者提供)
此外，我们可以访问线性代数运算，例如点。要访问这些，您需要访问 linalg 名称空间，它存储在 xt 名称空间中，如下所示:
xt::linalg::dot(array, arraytwo)
为了获得这个名称空间，您可能还需要使用 XTensor-blas。
[## xtensor-stack/xtensor-blas
xtensor-blas 是 xtensor 库的扩展，通过 cxxblas 和…提供 blas 和 LAPACK 库的绑定。
github.com](https://github.com/xtensor-stack/xtensor-blas)
结论
XTensor 是一个非常棒的工具，它将 Numpy 的易用性和标准化带到了 C++的前沿，用于科学计算和处理矩阵。一般来说，C++是一种很好的语言，因为它的可靠性——在这一点上已经有几十年的历史了，并且在它的生命周期中得到了很好的使用。
此外，我认为 XTensor 为有兴趣在数据科学和统计应用程序中使用 C++的数据科学家创造了一个很好的出口。如果没有别的，XTensor 就是一个非常酷的库，它有可能在 Python 用户和 C++用户之间架起一座桥梁。
Pytorch 中的线性分类
利用二元分类法检测乳腺癌
这篇中型文章将探索 Pytorch 库以及如何实现线性分类算法。我们将在经典且易于理解的数据集上应用该算法。在本文中，您将更加熟悉 Pytorch 的基础知识。之后，你应该准备好深入更高级的问题(例如时间序列)。
文章的结构如下:
- 加载数据集
- 拆分数据
- 预处理数据
- 创建模型
- 全梯度下降
- 模型评估
加载数据集
将使用的数据集是乳腺癌威斯康辛数据集，它非常适合演示二元分类。
问题很简单。科学家们收集了患有(恶性)或没有(良性)乳腺癌的患者的特征。我们的工作是辨别是非。这可以帮助医生在日常工作中处理数据。乳腺癌作为一个例子的原因是因为它的相关性。2018 年，据估计全球有 627，000 名女性死于乳腺癌[ 世卫组织，2020]。
from sklearn.datasets import load_breast_cancerdata = load_breast_cancer()
X, Y = data.data, data.target
在该数据集中，有 30 个特征(例如，平均半径、平均纹理)。
print(data.feature_names)['mean radius' 'mean texture' 'mean perimeter' 'mean area' 'mean smoothness' 'mean compactness' 'mean concavity' 'mean concave points' 'mean symmetry' 'mean fractal dimension' 'radius error' 'texture error' 'perimeter error' 'area error' 'smoothness error' 'compactness error' 'concavity error' 'concave points error' 'symmetry error' 'fractal dimension error' 'worst radius' 'worst texture' 'worst perimeter' 'worst area' 'worst smoothness' 'worst compactness' 'worst concavity' 'worst concave points' 'worst symmetry' 'worst fractal dimension']
目标映射很简单。这些特征表明你患有恶性或良性肿瘤。大多数肿瘤是良性的，这意味着它们不是癌性的，不会杀死你。
print(data.target_names)['malignant' 'benign']
拆分数据
from sklearn.model_selection import train_test_splitX_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.3)
首先，使用训练数据创建模型，然后使用测试数据进行测试。如果你熟悉这个领域，这应该是一个惊喜。因此，上面的代码使用 sklearn 将数据分成随机的训练和测试子集。一行代码，但是非常重要。
预处理数据
from sklearn.preprocessing import StandardScalerscaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)
为了获得更好的结果，一点缩放是很重要的。我们不希望一个特性支配另一个特性。我们能做的是确保数据被很好地分配。这里没有什么新东西，如果你不熟悉标准化功能，机器学习课程可能会有所帮助(例如 deeplearning.ai )。
你可以自己做，但为什么要重新发明轮子呢？对于这个任务，可以使用 sklearn 的 StandardScaler 。这是一个微妙的细节，但以与训练集相同的方式扩展测试集是很重要的。因为训练集比测试集大得多，而且两者都是随机选择的，所以它应该是对实际分布的更好估计。
创建模型
class BinaryClassification(torch.nn.Module):def __init__(self, input_dimension):
super().__init__()
self.linear = torch.nn.Linear(input_dimension, 1)def forward(self, input_dimension):
return self.linear(input_dimension)_, input_dimension = X_train.shape
model = BinaryClassification(input_dimension)
感兴趣的模型是一个非常简单的线性模型。在示例中，它展示了如何创建自己的 Pytorch“模块”,但是您也可以使用一行程序来完成。同样的事情也可以这样做。因此，如果您更喜欢这样，请随意使用一行程序。但是，本文意在告知。在 Pytorch 中创建一个定制的神经网络就是这么容易。
model = torch.nn.Linear(input_dimension, 1)
还有一点你可能已经注意到了，我们没有使用 sigmoid 激活函数，但这是可以解释的。不使用与 sigmoid 函数结合的二进制交叉熵损失函数，这里更倾向于使用具有 logits 损失的函数的二进制交叉熵。后者在数值上更稳定，从而导致更好的结果。
def configure_loss_function():
return torch.nn.BCEWithLogitsLoss()def configure_optimizer(model):
return torch.optim.Adam(model.parameters())
最后，使用了 Adam 优化器。可能有更好的优化算法，但它的表现相当好。还不需要用别人做实验。
模型评估
二元分类示例
剩下的代码只是一个完整的梯度下降循环以及训练和测试精度的计算。在每个时期都会发生以下步骤:
- 向前通过二进制分类模型。
- 损失函数测量 BCEWithLogits 损失。
- 损耗的梯度被重置为零。
- 损失函数向后传播(计算损失的梯度)。
- Adam optimizer 朝着“正确”的方向前进。
此时不给出完整的代码是愚蠢的。你可以自己随意摆弄算法和数据集。
提示:使用谷歌 Colab，它就像类固醇上的木星笔记本——也非常容易上手。
结论
二进制分类算法的结果在训练集和测试集上都有大约 98%的准确率——很神奇，不是吗？这也不需要太多的努力。如果你对这篇文章有任何改进，请在下面的评论中告诉我。这将有助于我和未来的读者，提前谢谢你。
线性判别分析——导论
探索 LDA 作为降维和分类技术
LDA 是一种降维算法，类似于 PCA。然而，虽然 PCA 是一种无监督算法，其重点在于最大化数据集中的方差，但是 LDA 是一种监督算法，其最大化类别之间的可分性。
杰米·麦卡弗里拍摄，版权所有。
在建立分类问题时，目标是确保类别的最大可分性或区分度。
假设我们有一个包含两列的数据集—一个解释变量和一个二进制目标变量(值为 1 和 0)。二进制变量的分布如下:
目标变量的分布(在 draw.io 上创建)
绿点代表 1，红点代表 0。因为只有一个解释变量，所以用一个轴(X)表示。但是，如果我们试图放置一个线性分割线来划分数据点，我们将无法成功地做到这一点，因为这些点分散在整个轴上。
划分点的分隔线(在 draw.io 上创建)
因此，似乎一个解释变量不足以预测二元结果。因此，我们将引入另一个特征 X2，并检查点在 2 维空间的分布。
二维空间中点的分布(在 draw.io 上创建)
在二维空间中，输出的分界似乎比以前更好。但是，在已经有多个要素的数据集中增加维度可能不是一个好主意。在这些情况下，LDA 通过最小化维度来拯救我们。在下图中，目标类被投影到一个新轴上:
目标值在变换轴上的投影(在 draw.io 上创建)
这些阶层现在很容易划分。LDA 将原始特征变换到一个新的轴上，称为线性判别式(LD ),从而降低维数并确保类的最大可分性。
为了将这种可分性用数字表示，我们需要一个度量可分性的标准。计算这两个等级的平均值之间的差异可以是这样一种度量。差值越大，则表明两点之间的距离越大。然而，这种方法没有将数据的扩散考虑在内。因此，即使较高的平均值也不能确保某些类不会相互重叠。
第二项措施是考虑类内的均值和方差。为了确保最大的可分性，我们将最大化均值之间的差异，同时最小化方差。
类内方差(S1 和 S2)和两个类的均值(M1 和 M2)
分数的计算方法是(M1-M2) /(S1 +S2)。
如果有三个解释变量- X1、X2、X3，LDA 将把它们转换成三个轴-LD1、LD2 和 LD3。这三个轴将根据计算出的分数排名第一、第二和第三。
因此，LDA 帮助我们降低维度和分类目标值。
我们将通过一个例子来看看 LDA 是如何实现这两个目标的。下面的数据显示了 IBM 的一个虚构数据集，它记录了员工数据和减员情况。目标是根据不同的因素，如年龄、工作年限、出差性质、教育程度等，预测员工的流失情况。
sns.countplot(train['Attrition'])
plt.show()
大约有 1470 份记录，其中 237 名员工离开了组织，1233 名员工没有离开。是被编码为 1，否被编码为 0。
如果没有正确预测，员工流失会导致失去有价值的人，导致组织效率降低，团队成员士气低落等。因此，有必要正确预测哪个员工可能会离开。换句话说，如果我们预测某个员工会留下来，但实际上该员工离开了公司，那么漏报的数量就会增加。我们的目标是最小化假阴性，从而提高召回率(TP/(TP+FN))。
我们将尝试使用 KNN 对类别进行分类:
knn=KNeighborsClassifier(n_neighbors=10,weights='distance',algorithm='auto', p=3)
start_time = time.time()
knn.fit(X_train_sc,y_train)
end_time = time.time()
time_knn = end_time-start_time
print(time_knn)
pred = knn.predict(X_test_sc)
print(confusion_matrix(y_test,pred))
print(classification_report(y_test,pred))
适应 KNN 所用的时间:0 . 46360 . 38386838661
对于 0.05 离开的员工来说，回忆是非常差的。
我们现在将使用 LDA 作为分类算法并检查结果。
lda_0 = LDA()
lda_0.fit(X_train_sc, y_train)
y_test_pred_0 = lda_0.predict(X_test_sc)
print(confusion_matrix(y_test, y_test_pred_0))
print(classification_report(y_test, y_test_pred_0))
召回率急剧上升至 0.36
最后，我们将使用 LDA 转换训练集，然后使用 KNN。
#Transformation by LDA
lda_1 = LDA(n_components = 1, solver='eigen', shrinkage='auto')
X_train_lda = lda_1.fit_transform(X_train_lda, y_train_lda)
X_test_lda = lda_1.transform(X_test_lda)
接下来是标准缩放。
sc3 = StandardScaler()
X_train_lda_sc = sc3.fit_transform(X_train_lda)
X_test_lda_sc = sc3.transform(X_test_lda)
现在我们对转换后的数据应用 KNN。
knn=KNeighborsClassifier(n_neighbors=8,weights='distance',algorithm='auto', p=3)
start_time = time.time()
knn.fit(X_train_lda_sc,y_train_lda)
end_time = time.time()
time_lda = end_time-start_time
print(time_lda)
print(confusion_matrix(y_test_lda,pred))
print(classification_report(y_test_lda,pred))
对转换后的数据运行 KNN 程序所用的时间:0 . 46866 . 46868686661
召回率现在提高了 6%。
此外，KNN 拟合 LDA 转换数据所用的时间是 KNN 一个人所用时间的 50%。
全部代码可在这里获得。
希望我已经演示了 LDA 的使用，包括分类和将数据转换成不同的轴！
敬请关注更多内容！一如既往，感谢您的任何反馈。
线性判别分析，已解释
里面的艾
直觉、插图和数学:它如何不仅仅是一个降维工具，以及为什么它对于现实世界的应用是健壮的。
通过混合判别分析(MDA)学习的边界(蓝线)成功地分离了三个混合类。MDA 是 LDA 的强大扩展之一。
关键要点
- 线性判别分析不仅是一种降维工具，也是一种稳健的分类方法。
- 无论有没有数据正态性假设，我们都可以得到相同的 LDA 特征，这解释了它的鲁棒性。
线性判别分析被用作分类、降维和数据可视化的工具。它已经存在了一段时间了。尽管 LDA 很简单，但它经常产生健壮、体面且可解释的分类结果。在处理现实世界的分类问题时，LDA 通常是在使用其他更复杂和更灵活的方法之前的基准方法。
使用 LDA(及其变体)的两个突出例子包括:
- 破产预测 : Edward Altman 的 1968 年模型使用训练好的 LDA 系数预测公司破产的概率。根据 31 年的数据评估，准确率在 80%到 90%之间。
- 面部识别:虽然从主成分分析(PCA)中学习到的特征被称为特征脸，但从 LDA 中学习到的特征被称为鱼脸，以统计学家罗纳德·费雪爵士的名字命名。我们稍后解释这种联系。
本文首先介绍了经典的 LDA，以及为什么它作为一种分类方法根深蒂固。接下来，我们看到这种方法中固有的降维，以及它如何导致降秩 LDA。在那之后，我们看到费希尔是如何熟练地得出同样的算法，而没有对数据做任何假设。一个手写数字分类问题被用来说明 LDA 的性能。最后总结了该方法的优缺点。
本文改编自我的一篇博文。如果你喜欢 LaTex 格式的数学和 HTML 风格的页面，你可以在我的博客上阅读这篇文章。此外，还有一组相应的幻灯片，我在其中展示了大致相同的材料，但解释较少。
判别分析分类
让我们看看 LDA 作为一种监督分类方法是如何导出的。考虑一个一般的分类问题:一个随机变量 X 来自 K 类中的一个，带有一些特定于类的概率密度 f ( x )。判别规则试图将数据空间划分为代表所有类别的 K 个不相交区域(想象棋盘上的盒子)。对于这些区域，通过判别分析进行分类仅仅意味着如果 x 在区域 j 中，我们将 x 分配到类别 j 中。问题是，我们如何知道数据 x 属于哪个区域？自然，我们可以遵循两个分配规则:
- 最大似然法则:如果我们假设每一类都以相等的概率发生，那么分配 x 给类 j 如果
- 贝叶斯法则:如果我们知道类先验概率 π ，那么分配 x 给类 j 如果
线性和二次判别分析
如果我们假设数据来自多元高斯分布，即 X 的分布可以用其均值(【μ)和协方差(【σ)来表征，就可以得到上述分配规则的显式形式。遵循贝叶斯规则，如果对于 i = 1，…， K ，数据 x 在所有 K 类中具有最高的可能性，我们将数据 j 分类:
上述函数称为判别函数。注意这里对数似然的使用。换句话说，判别函数告诉我们数据 x 来自每个类的可能性有多大。因此，分离任意两个类的判定边界 k 和 l 是两个判别函数具有相同值的 x 的集合。因此，任何落在决策边界上的数据都同样可能来自这两个类(我们无法决定)。
LDA 出现在我们假设在 K 类中协方差相等的情况下。也就是说，不是每个类一个协方差矩阵，而是所有类都有相同的协方差矩阵。那么我们可以得到下面的判别函数:
注意这是 x 中的线性函数。由此可见，任意一对类之间的决策边界也是 x 中的线性函数，其得名原因:线性判别分析。在没有相等协方差假设的情况下，似然中的二次项不会抵消，因此得到的判别函数是 x 中的二次函数:
在这种情况下，判定边界在 x 中是二次的。这被称为二次判别分析(QDA)。
哪个更好？艾达还是 QDA？
在实际问题中，总体参数通常是未知的，只能通过训练数据作为样本均值和样本协方差矩阵来估计。虽然与 LDA 相比，QDA 适应更灵活的决策边界，但是需要估计的参数数量也比 LDA 增加得更快。对于 LDA，需要 (p+1) 个参数来构造(2)中的判别函数。对于具有 K 个类的问题，我们将只需要 (K-1) 这样的判别函数，通过任意选择一个类作为基类(从所有其他类中减去基类可能性)。因此，LDA 的估计参数的总数是 (K-1)(p+1) 。
另一方面，对于(3)中的每个 QDA 判别函数，需要估计均值向量、协方差矩阵、类先验:
-均值: p
-协方差: p(p+1)/2
-类先验: 1
同样，对于 QDA，需要估计***(K-1){ p(p+3)/2+1 }***参数。
因此，LDA 中估计的参数数量随 p 线性增加，而 QDA 的参数数量随 p 二次增加。当问题维数较大时，我们预计 QDA 的性能会比 LDA 差。
两全其美？LDA 和 QDA 之间的妥协
我们可以通过正则化各个类别的协方差矩阵来找到 LDA 和 QDA 之间的折衷。正则化意味着我们对估计的参数施加一定的限制。在这种情况下，我们要求各个协方差矩阵通过惩罚参数(例如α:
公共协方差矩阵也可以通过罚参数例如β朝着单位矩阵正则化:
在输入变量的数量大大超过样本数量的情况下，协方差矩阵的估计可能很差。收缩有望提高估计和分类精度。下图说明了这一点。
有收缩和无收缩 LDA 的性能比较。归功于 scikit-learn 。
生成上图的脚本。
LDA 的计算
从(2)和(3)可以看出，如果我们先对角化协方差矩阵，鉴别函数的计算可以简化。也就是说，数据被转换为具有相同的协方差矩阵(无相关性，方差为 1)。对于 LDA，我们是这样进行计算的:
步骤 2 将数据球形化，以在变换的空间中产生单位协方差矩阵。步骤 4 是通过下面的(2)获得的。
我们举两个类的例子，看看 LDA 到底在做什么。假设有两个类， k 和 l 。我们将 x 归类到 k 类，如果
上述条件意味着类 k 比类 l 更有可能产生数据 x 。按照上面概述的四个步骤，我们编写
也就是说，我们将数据 x 分类到类别 k 如果
导出的分配规则揭示了 LDA 的工作原理。等式的左侧(l.h.s .)是 x* 在连接两类平均值的线段上的正交投影长度。右手边是由类别先验概率校正的片段的中心位置。本质上，LDA 将球形数据分类到最接近的类均值。这里我们可以做两个观察:
- 当类别先验概率不相同时，判定点偏离中间点，即，边界被推向具有较小先验概率的类别。
- 数据被投影到由类平均值(规则的 l . h . s .**x ***的乘法和平均值减法)跨越的空间上。然后在该空间中进行距离比较。
降秩 LDA
我刚才描述的是用于分类的 LDA。LDA 还因其能够找到少量有意义的维度而闻名，使我们能够可视化并处理高维问题。我们所说的有意义是什么意思，LDA 是如何找到这些维度的？我们将很快回答这些问题。
首先，看看下面的情节。对于具有三种不同类型的葡萄酒和 13 个输入变量的葡萄酒分类问题，该图将数据可视化在 LDA 找到的两个判别坐标中。在这个二维空间中，各个阶层可以很好地区分开来。相比之下，使用主成分分析发现的前两个主成分，类别没有被清楚地分开。
LDA 固有的降维
在上面的葡萄酒例子中，一个 13 维的问题在 2d 空间中被可视化。为什么会这样？这是可能的，因为 LDA 有一个固有的降维。从上一节中我们已经观察到，LDA 在不同类均值所跨越的空间中进行距离比较。两个不同的点位于一维直线上；三个不同的点位于 2d 平面上。同样， K 类意味着位于一个维数最多为 (K-1) 的超平面上。特别地，该装置所跨越的子空间是
当在该空间中进行距离比较时，与该子空间正交的距离不会添加任何信息，因为它们对每个类别的贡献是相等的。因此，通过将距离比较限制到这个子空间，不会丢失任何对 LDA 分类有用的信息。这意味着，通过将数据正交投影到这个子空间上，我们可以安全地将我们的任务从一个 p 维问题转换为一个 (K-1) 维问题。当 p 远大于 K 时，这是维度数量的一个相当大的下降。
如果我们想进一步减少尺寸，从 p 到 L ，比如二维的 L = 2 怎么办？我们可以尝试从 (K-1) 维空间构造一个 L 维子空间，并且这个子空间在某种意义上对于 LDA 分类是最优的。
最佳子空间是什么？
Fisher 提出，当空间中球形数据的类均值在方差方面具有最大分离时，小得多的 L 维子空间是最佳的。根据这个定义，通过对球形类均值进行 PCA，可以简单地找到最佳子空间坐标，因为 PCA 找到最大方差的方向。计算步骤总结如下:
通过这个过程，我们将数据从 X 替换到 Z ，并将问题维度从 p 减少到 L 。通过设置 L = 2 ，通过该程序找到前一葡萄酒图中的判别坐标 1 和 2。使用新数据重复先前的分类 LDA 过程被称为降秩 LDA。
费希尔的 LDA
如果之前降秩 LDA 的推导看起来与你之前知道的非常不同，你并不孤单！启示来了。Fisher 根据他的最优性定义以不同的方式导出了计算步骤。他执行降秩 LDA 的步骤后来被称为 Fisher 判别分析。
费希尔没有对数据的分布做任何假设。相反，他试图找到一个“合理的”规则，以便分类任务变得更容易。特别地，费希尔找到了原始数据的线性组合，其中类间方差 B = cov( M )相对于类内方差 W 被最大化，如(6)中所定义的。
下面的图来自 ESL，显示了为什么这个规则有直观的意义。该规则开始寻找一个方向， a ，其中，在将数据投影到该方向上之后，类意味着它们之间具有最大间隔，并且每个类在其内部具有最小方差。在这个规则下找到的投影方向，如右图所示，使得分类更加容易。
在费雪的“合理法则”下找到的投影方向如右图所示。
寻找方向:费希尔的方式
使用费希尔的合理规则，找到最佳投影方向相当于解决一个优化问题:
回想一下，我们希望找到一个方向，使类间方差最大化(分子)，类内方差最小化(分母)。这可以转化为广义特征值问题。对于那些感兴趣的人，你可以在我最初的博客文章中找到这个问题是如何解决的。
求解后得到最优子空间坐标，也称为判别坐标，作为 inv(W)∏B的特征向量。可以看出，这些坐标与上述降秩 LDA 公式中从 X 到 Z 的坐标相同。
令人惊讶的是，与降秩 LDA 公式不同，Fisher 在没有对总体进行任何高斯假设的情况下得出这些判别坐标。人们希望，有了这个合理的规则，即使数据不完全符合高斯分布，LDA 也能表现良好。
手写数字问题
下面的例子展示了 Fisher’s LDA(简称 LDA)的可视化和分类能力。我们需要使用 64 个变量(来自图像的像素值)来识别 10 个不同的数字，即 0 到 9。数据集取自这里的。首先，我们可以将训练图像可视化，它们看起来像这样:
接下来，我们在前半部分数据上训练 LDA 分类器。解决前面提到的广义特征值问题给了我们一个最佳投影方向的列表。在这个问题中，我们保留了前四个坐标，转换后的数据如下所示。
上面的图允许我们解释训练好的 LDA 分类器。例如，坐标 1 有助于对比 4 和 2/3，而坐标 2 有助于对比 0 和 1。随后，坐标 3 和 4 有助于区分坐标 1 和 2 中没有很好分开的数字。我们使用数据集的另一半来测试训练好的分类器。下面的报告总结了结果。
最高精度为 99%，最低精度为 77%，这是一个相当不错的结果，因为该方法是在大约 70 年前提出的。此外，我们没有做任何事情来改善这个具体问题的程序。例如，输入变量中存在共线性，收缩参数可能不是最佳的。
LDA 概述
LDA 的优点:
- 简单的原型分类器:使用到类均值的距离，这很容易解释。
- 决策边界是线性的:实现简单，分类稳健。
- 降维:它提供了关于数据的信息性低维视图，这对可视化和特征工程都很有用。
LDA 的缺点:
- 线性决策边界可能无法充分区分类别。希望支持更通用的边界。
- 在高维设置中，LDA 使用了太多的参数。LDA 的正则化版本是期望的。
- 希望支持更复杂的原型分类。
谢谢你一直读到最后！在下一篇文章中，我们将介绍灵活的、惩罚的和混合判别分析来解决 LDA 的三个缺点。有了这些概括，LDA 可以处理更加困难和复杂的问题，例如特征图像中显示的问题。
如果你对更多统计学习的东西感兴趣，可以随意看看我的其他文章:
CNN 有什么独特之处，卷积到底是做什么的？这是一个无数学介绍的奇迹…
towardsdatascience.com](/a-math-free-introduction-to-convolutional-neural-network-ff38fbc4fc76) [## 引擎盖下:什么联系线性回归，岭回归，主成分分析？
从挽救病态回归问题，使快速计算的正则化路径，这是…
towardsdatascience.com](/under-the-hood-what-links-ols-ridge-regression-and-pca-b64fcaf37b33)
参考
- R.a .、费希尔、、分类问题中多重测量的使用、 (1936)，《优生学年鉴》，7 卷 2 期，179–188 页。
- J.Friedman，T. Hastie 和 R. Tibshirani， 统计学习的要素 (2001)，统计中的斯普林格系列。
- K.V. Mardia，J. T. Kent 和 J. M. Bibby，多元分析 (1979)，概率与数理统计，学术出版社。
Python 中的线性插值一行代码
推导实用算法
利用numpy简化每个人工具箱中应该有的技术
插值对数线性和反转(线性对数)值
介绍
线性插值从离散数据中创建一个连续函数。这是梯度下降算法的基础构建块，该算法用于几乎所有机器学习技术的训练。尽管深度学习非常复杂，但没有它就无法工作。这里我们先来看看构造线性插值算法的理论基础，然后实现一个实用的算法并将其应用到一个案例中。
衍生物
如果你想跳过推导，向下滚动到标题为的章节，一句话作为笑点。不然继续看！这不会花那么长时间…
在之前的文章中，我们通过使用一些线性代数讨论了线性回归的基础。特别是，我们用这个结果来确定回归系数:
我们也有这个方程来预测这些系数:
如果我们将第一个方程代入第二个方程，我们将得到预测值的直接估计值:
所有这些 X 项形成了一个叫做投影矩阵的东西，也叫做帽子矩阵:
投影矩阵 将 (或 投影 )的观测值转换为预测值。为此，它必须包含对任意值 Y 进行预测所需的所有信息。**
如此看来，如果我们知道了，我们就知道了如何预测任何值！为了说明这一点，让我们创建一个玩具示例，看看帽子矩阵是什么样子的。我们将使用位于坐标(0，1)和(5，2)的两个数据点。**
发展
让我们用 Python 演示一下:
***import numpy as np
np.set_printoptions(suppress=True) # no scientific formatting# np.r_ and np.c_ are short-hands for concatenation of arrays
# np.r_ performs along the first-axis (rows)
# np.c_ performs along the second-axis (columns)# our x-values
x = np.c_[1., 1.,], [0., 5.]]# our y-values
y = np.r_[1., 2.]print(x)
array([[1., 0.],
[1., 5.]])***
到目前为止很简单。我们可以使用第一个等式求解回归系数:
***# @ is the numpy operator for matrix multiplication
c = np.linalg.inv(x.T @ x) @ x.T @ Yprint(c)
array([1., 0.2])***
正如我们所料，斜率为 1/5，截距为 1。下一个是投影矩阵:
***# @ is the numpy operator for matrix multiplication
h = x @ np.linalg.inv(x.T @ x) @ x.Tprint(h)
array([[1., 0.],
[0., 1.]])***
回想一下，当我们将 H 乘以 Y 时，我们得到了预测值【ŷ】。因为我们在对角线上有 1，在其他地方有 0，在这种情况下， H 实际上是单位矩阵，所以我们简单地返回实际 y 值(1，2)的预测值(1，2)。当我们的点不像几乎所有线性回归问题那样形成一条直线时，情况就不一样了。然而，对于插值，我们总是可以做这样的假设，因为只有两个点限制了我们想要插值的点！
插入文字
实际插值呢？如果我们回忆一下第二个方程，我们可以为X 提供任何值以便创建一个预测值。让我们利用这个事实重写我们的投影矩阵，用 X₀ 代表我们的实际值，用 X₁ 代表我们要插值的值。我们将使用 W 将其表示为权重矩阵，原因我们稍后会看到:
如果我们在 0 和 5 的值之间尝试 X₁ 的各种值会发生什么？
**for i in np.linspace(0, 5, 6):
w = np.c_[1., i] @ np.linalg.inv(x.T @ x) @ x.T
y_hat = w @ y
print(f'x1 = {i:.1f}, w = {w}, y_hat = {y_hat}')x1 = 0.0, w = [[1.0 0.0]], y_hat = [1.0]
x1 = 1.0, w = [[0.8 0.2]], y_hat = [1.2]
x1 = 2.0, w = [[0.6 0.4]], y_hat = [1.4]x1 = 3.0, w = [[0.4 0.6]], y_hat = [1.6]
x1 = 4.0, w = [[0.2 0.8]], y_hat = [1.8]x1 = 5.0, w = [[0.0 1.0]], y_hat = [2.0]**
有趣的事情发生了！在每种情况下，权重矩阵的和是 1。当我们将这个权重矩阵乘以 Y 值时，我们得到的是两个值的加权平均值。第一个权重与两个值之间插值的的距离成正比，第二个权重正好是第一个权重的补码。换句话说:**
一句话
这给了我们一行的线性插值:
***new_y = np.c_[1., new_x] @ np.linalg.inv(x.T @ x) @ x.T @ y***
当然，这有点噱头。我们必须准确地知道原始 x 值数组中的两个值，新的插值 x 值位于这两个值之间。我们需要一个函数来确定这两个值的索引。值得庆幸的是，numpy恰恰包含了这样一个恰好的函数:np.searchsorted。让我们用它把线性代数变成一个矢量化的函数。我们不是为每个插值计算权重矩阵，而是存储相关信息:权重本身，以及 X 值的索引(当我们计算插值时，它们也将由 Y 值的索引决定)。
**from typing import NamedTupleclass Coefficient(NamedTuple):
lo: np.ndarray # int
hi: np.ndarray # int
w: np.ndarray # floatdef interp_coef(*x0:* np.ndarray, *x*: np.ndarray) -> Coefficient:
# find the indices into the original array
hi = np.minimum(len(x0) - 1, np.searchsorted(x0, x, 'right'))
lo = np.maximum(0, hi - 1)
# calculate the distance within the range
d_left = x - x0[lo]
d_right = x0[hi] - x
d_total = d_left + d_right # weights are the proportional distance
w = d_right / d_total # correction if we're outside the range of the array
w[np.isinf(w)] = 0.0 # return the information contained by the projection matrices
*return* Coefficient(lo, hi, w)def interp_with(*y0*: np.ndarray, *coef*: Coefficient) -> np.ndarray:
*return* coef.w * y0[coef.lo] + (1 - coef.w) * y0[coef.hi]**
这个算法的方便之处在于我们只需要计算一次我们的权重。有了权重，我们就可以对与 x 值向量对齐的任何 y 值向量进行插值。
示例应用程序
在工程问题中，我们正在建模的物理过程通常是指数函数或幂律函数。当我们为这些过程收集数据时，我们可以以对数间隔的时间间隔进行。许多解释方法被构造成以均匀的测井间隔时间间隔工作。然而，一些数值方法更便于处理均匀的线性间距值。所以来回转换的能力很方便。线性插值是执行这种变换的一种工具。
**事实上，我们并不局限于线性插值。我们可以首先线性化我们的值，然后进行插值。让我们来看一个径向系统中流量、 、q 、压力、 **p 、的例子。这些数据是用均匀对数间隔的样本合成的，当我们将它绘制在双对数图上时，我们会看到大多数数据的不同幂律行为。在后期，流速转变为指数下降，但这不会对我们的插值产生太大影响。
径流系统的压力和流量数据
让我们使用一组等距值进行插值。我们将假设幂律函数作为插值的基础，首先对我们的值应用对数变换。我们还将反转插值，以验证插值是否正确执行:
**# construct out linear time series
tDi = np.linspace(1e-3, 1e4, 100_001)# interpolate from log to linear
coef = interp_coef(np.log(time), np.log(time_interp))
pi = np.exp(interp_with(np.log(pressure), coef))
qi = np.exp(interp_with(np.log(rate), coef))# reverse from linear to log
coef_i = interp_coef(np.log(time_interp), np.log(time))
pr = np.exp(interp_with(np.log(pi), coef_i))
qr = np.exp(interp_with(np.log(qi), coef_i))**
插值对数线性和反转(线性对数)值
摘要
线性方法是处理数据的重要工具。变量的变换，例如通过使用示例中所示的np.log，允许我们在不期望数据是线性的情况下应用这些技术。
虽然我们推导并实现了自己的算法，但标准算法确实存在:scipy.interpolated.interp1d。然而，为了在一些包中使用，比如当使用pymc3时，可能需要为自己实现类似这样的东西。
还存在更一般的插值方法。其中之一是 B 样条，我们将在本系列的下一篇文章中讨论它。b 样条也可以执行线性插值，以及二次，三次等。b 样条通过扩展我们在这里开发的概念来工作。也就是说，将权重 向量 扩展为权重 矩阵 。我们无需编写函数来查找每个插值的索引，只需用系数和新 x 值的矩阵乘法来计算插值的 y 值即可。
从零开始的线性混合模型
生命科学的数理统计和机器学习
使用最大似然法推导和编码 LMM
这是来自专栏 生命科学的数理统计和机器学习 的第十八篇文章，我试图以简单的方式解释生物信息学和计算生物学中使用的一些神秘的分析技术。线性混合模型(也称为线性混合效应模型)在生命科学中广泛使用，有许多教程显示如何在 R 中运行该模型，但是有时不清楚在似然最大化过程中随机效应参数是如何优化的。在我之前的文章 线性混合模型如何工作 中，我介绍了模型的概念，在本教程中，我们将应用最大似然(ML) 方法从头开始推导并编码线性混合模型(LMM)，即我们将使用普通 R 编码 LMM，并将输出与来自 lmer 和 lme R 函数的输出进行比较。本教程的目标是解释 LMM“喜欢我的祖母”，暗示没有数学背景的人应该能够理解 LMM 在幕后做什么。
玩具数据集
让我们考虑一个玩具数据集，它非常简单，但仍然保留了线性混合建模(LMM)的典型设置的所有必要元素。假设我们只有 4 个数据点 /样本 : 2 个来自个体#1 ，另外 2 个来自个体#2 。此外，这 4 个点分布在两种状态之间:未处理和已处理。假设我们测量了每个个体对治疗的反应(或)，并且想要说明治疗是否导致研究中个体的显著反应。换句话说，我们的目标是实施类似于 配对 t 检验 的**，并评估治疗的意义。稍后，我们将把来自 LMM 和配对 t 检验的输出联系起来，并表明它们确实是相同的**。在玩具数据集中， Treat 栏中的 0 表示“未处理”，1 表示“已处理”。首先，我们将使用简单的普通最小二乘法(OLS) 线性回归，它不考虑数据点之间的相关性。
****
从技术上来说，这是可行的，但是，这不是一个很好的匹配，我们有一个问题 这里。普通最小二乘(OLS)线性回归假设所有观察值(图上的数据点)都是独立的**，这将导致不相关的，因此 正态分布残差 。然而，我们知道图上的数据点属于 2 个个体，即每个个体 2 个点。原则上，我们可以分别为每个个体拟合一个线性模型。然而，这也不是一个很好的契合。我们对每个人都有两点，所以太少了，无法对每个人进行合理的拟合。此外，正如我们之前在和中看到的，个体拟合并不能说明总体/群体概况，因为与其他个体拟合相比，其中一些拟合可能具有相反的行为。******
相比之下，如果我们想要将所有四个数据点 拟合在一起，我们将需要以某种方式说明它们是而非独立的事实，即它们中的两个属于个体#1，两个属于个体#2。这可以在线性混合模型(LMM)或配对检验中完成，例如 配对 t 检验 (参数)或 Wilcoxon 符号秩检验(非参数)。
具有 Lmer 和 Lme 的线性混合模型
当观测值之间存在非独立性时，我们使用 LMM。在我们的例子中，观察集中在个体内部。让我们对斜率和截距应用具有固定效应的 LMM，对截距应用随机效应，这将导致在 Resp～Treat 公式中增加一个 (1 | Ind) 项:
这里 REML=FALSE 仅仅意味着我们正在使用传统的最大似然(ML)优化，而不是受限最大似然(我们将在另一个时间讨论 REML)。在 lmer 输出的随机效应部分，我们看到最小化的 2 个参数的估计值:剩余方差对应于标准偏差(标准差)4.243，以及随机效应(个体间共享)方差与标准偏差为 5.766 的截距相关。类似地，在 lmer 输出的固定效果部分，我们可以看到两个估计值:1)截距等于 6.5，以及 2)斜率/ Treat 等于 9。因此，我们有 4 个优化参数，对应于 4 个数据点。如果我们查看玩具数据集部分中的第一个数字，并意识到未处理样本的两个值的平均值为(3 + 10) / 2 =6.5，则固定效应的值是有意义的，我们将其表示为 β 1 ，而处理过的样本的平均值为(6 + 25) / 2 = 15.5，我们将其表示为 β 2 。后者相当于 6.5 + 9，即截距固定效应的估计值(=6.5)加上斜率固定效应的估计值(=9)。这里，我们注意随机和固定效应的精确值，因为我们将在以后推导和编码 LMM 时再现它们。
默认情况下， lme4 R 包和 lmer R 函数不提供具有统计显著性的测量，例如 p 值，但是，如果您仍然希望您的 LMM 拟合具有 p 值，则可以使用 nlme R 包中的 lme 函数:
同样，这里我们有截距( StdDev = 5.766264 )和残差( StdDev = 4.242649 )的随机效果，以及截距(值= 6.5)和斜率/处理(值= 9)的固定效果。非常有趣的是，固定效应的标准误差及其 t 值(Treat 的 t 值=1.5)在 lmer 和 lme 之间并不一致。然而，如果我们在 lmer 函数中要求 REML=TRUE，则包括 t 值的固定效应统计在 lme 和 lmer 之间是相同的**，然而随机效应统计是不同的。**
这就是最大似然法(ML)和限制最大似然法(REML)** 之间的区别，我们将在下次讨论。**
LMM 与配对 T 检验的关系
之前，我们说过 LMM 是简单配对 t 检验的一种更复杂的形式。让我们证明，对于我们的玩具数据集，它们确实给出了相同的输出。在途中，我们还将了解配对 t 检验和非配对 t 检验之间的技术差异。让我们首先在处理组和未处理组样本之间进行配对 t 检验，考虑它们之间的非独立性:
我们可以看到，配对 t 检验报告的 t 值=1.5 和 p 值= 0.3743 与 LMM 使用 nlme 函数或 REML = TRUE 的 lmer 获得的结果相同。配对 t 检验统计报告的“差异均值= 9”也与来自 lmer 和 nlme 的固定效应估计值一致，记住我们有处理估计值= 9，它只是处理和未处理样本均值之间的差异。****
现在，配对 t 检验到底在做什么？嗯，配对 t 检验的想法是使设置看起来像一个单样本 t 检验**，其中一组中的值被检验是否与零有显著偏差，这是第二组的一种平均值。换句话说，我们可以查看配对 t 检验，就好像我们将单个拟合的截距(见第一个图)或未处理组的平均值向下移动到零。简单地说，这相当于从已处理的 Resp 值中减去未处理的 Resp 值，即如下所示，将 Resp 变量转换为 Resp_std (标准化响应)，然后对 Resp_std 变量而不是 Resp 执行非配对 t 检验:**
我们观察到，对于 Treat = 0，即未治疗组，响应值变为 0，而治疗组(Treat=1)的响应值减少了未治疗组的值。然后，我们简单地使用新的 Resp_std 变量并运行非配对 t-test，结果相当于对原始 Resp 变量运行配对 t-test。因此，我们可以总结，LMM 复制了配对 t 检验的结果，但允许更大的灵活性，例如，不仅两个(如 t 检验)，但多组比较等。
线性混合模型背后的数学
现在让我们尝试使用我们的玩具数据集来推导几个 LMM 方程。我们将再次查看这 4 个数据点，并对治疗效果进行一些数学记法，*，这是除了固定效果**，以及由于两个个体之间的点聚集而导致的块式结构*，实际上是随机效果贡献。我们要用 β 和 u 参数来表示响应(Resp)坐标 y 。****
这里， β 1 是个体在未处理状态下的反应，而β2是对处理的反应**。也可以说， β 1 是未处理样本的平均值，而 β 2 是处理过的样本的平均值。变量 u 1 和 u 2 是分别说明个体#1 和个体#2 特定效果的块变量。最后，ϵij∽n(0，∑)是残差**，即我们无法建模的误差，只能试图将其最小化作为最大似然优化问题的目标。因此，我们可以把响应变量 y 记为参数 β ， u ，即固定和随机效果，以及 ϵ 作为 Eq。(1).在一般形式下，这个代数方程组可以重写为方程。其中指数 i = 1，2 对应于治疗效果，j = 1，2 描述个体效果。我们也可以用矩阵形式 Eq 来表示这个方程组。(3).因此，我们得出以下著名的 LMM 矩阵形式，它出现在所有的教科书中，但并不总是得到正确的解释。(4).******
这里， X 被称为设计矩阵，而 K 被称为块矩阵**，它编码了数据点之间的关系，即它们是否来自相关的个体，或者甚至像我们的情况一样来自同一个个体。值得注意的是，治疗被建模为固定效应，因为治疗-未治疗的水平用尽了所有可能的治疗结果。相比之下，数据的块式结构被建模为随机效应，因为个体是从群体中采样的，可能无法正确代表个体的整个群体。换句话说，存在与随机效应相关联的误差，即uj∾N(0， σs )** ，而固定效应被假定为无误差。例如，性别通常被建模为固定效应，因为它通常被假设为只有两个水平(男性，女性)，而生命科学中的批量效应应被建模为随机效应，因为潜在的额外实验方案或实验室会产生更多的，即许多水平的样本之间的系统差异，从而混淆数据分析。根据经验，人们可能会认为固定效应不应该有很多水平，而随机效应通常是多水平的分类变量，其中的水平只是所有可能性的一个样本，而不是全部。********
让我们继续推导数据点 Y 的均值和方差。由于随机效应误差和残差均来自均值为零的正态分布，而 E [ Y ]中的非零分量来自于固定效应，我们可以将 Y 的期望值表示为等式。(5).接下来，固定效应项的方差为零，因为固定效应被假设为无误差，因此，对于 Y 的方差，我们获得等式。(6).
情商。(6)是考虑到 var(ku)=k var(u) k**^t 和 var(ϵ)=σ*** I和 var(u)= σ s I 而得到的，其中 I 是一个 4×4的恒等式矩阵*。这里， σ 是残差方差(未建模/未缩减误差)，而 σ s 是随机截距效应(跨数据点共享)方差。 σ s 前面的矩阵称为亲属矩阵，由等式给出。(7).亲属矩阵编码了数据点之间的所有关联。例如，一些数据点可能来自基因相关的人，地理位置非常接近，一些数据点可能来自技术复制。这些关系被编码在亲属矩阵中。因此，等式中数据点的方差-协方差矩阵。(6)采用 Eq 的最终形式。(8).一旦我们获得了方差-协方差矩阵，我们就可以继续进行需要方差-协方差的最大似然函数的优化过程。****
最大似然(ML)原理的 LMM
为什么我们花这么多时间推导方差-协方差矩阵，它与线性回归有什么关系？原来，拟合线性模型的整个概念，以及许多其他传统频率统计的概念，如果不是全部的话，都来自于最大似然(ML) 原理。为此，我们需要最大化多元高斯 分布函数关于参数 β 1 、 β 2 、 σs 和 σ 、Eq。(9).
这里|σy|表示方差-协方差矩阵的行列式。我们看到方差-协方差矩阵的逆矩阵和行列式被明确地包含在似然函数中，这就是为什么我们必须通过随机效应方差 σs 和残差方差 σ 来计算其表达式。似然函数的最大化等同于对数似然函数的最小化。(10).**
我们将需要对方差-协方差矩阵的行列式、逆方差-协方差矩阵和逆方差-协方差矩阵与 Y-**Xβ--项的乘积执行冗长的符号推导。根据我的经验，这在 R / Python 中很难做到，但是我们可以使用**(或者类似的**Mathematica或者Matlab)进行符号计算，并推导出方差协方差矩阵的行列式和逆的表达式:****
复杂的符号推导在 Maple / Mathematica / Matlab 环境下变得容易
使用 Maple，我们可以获得方差-协方差矩阵的行列式，如等式。(11).接下来，Eq 中的最后一项。对于，对数似然**采用等式的形式。(12).**
现在，我们准备使用 optim R 函数，针对 β 1 、 β 2 、 σs 和 σ 执行对数似然函数的数值最小化:
我们可以看到，最小化算法已经成功地收敛，因为我们得到了“收敛= 0”消息。在输出中，∑= 4.242640687是残差标准差，它恰好再现了来自 lme 和 lmer 的结果(REML = FALSE)。以此类推，∑s= 5.766281297是共享标准差，其再次精确地再现了来自 lme 和 lmer(REML = FALSE)函数的相应随机效应截距输出。不出所料，固定效应 β 1 = 6.5 是协议中未处理样本的平均值，具有来自 lmer 和 lme 的截距固定效应估计值。接下来， β 2 = 15.5 是处理样本的平均值，它是来自 lmer 和 lme R 函数的截距固定效应估计值(= 6.5)加上斜率/处理固定效应估计值(= 9)。
出色的工作！通过从零开始推导和编码线性混合模型(LMM)，我们已经成功地从 lmer / lme 函数中再现了固定效果和随机效果输出！
摘要
在这篇文章中，我们学习了如何在一个 玩具数据集上导出并编码一个线性混合模型(LMM)**。我们介绍了 LMM 和配对 t 检验之间的关系，并从 lmer 和 lme R 函数中复制了固定和随机效应参数。******
在下面的评论中，让我知道哪些来自生命科学的分析技术对你来说是特别神秘的，我会在以后的文章中介绍它们。在我的 Github 上查看帖子中的代码。请在媒体关注我，在 Twitter @NikolayOskolkov 关注我，在 Linkedin 关注我。在下一篇文章中，我们将讨论最大似然法(ML)和受限最大似然法(REML) 的区别，敬请关注。
使用 Python 进行线性编程
一步一步的介绍使用 Python 中的纸浆库公式化和解决线性优化问题。
安托万·道特里在 Unsplash 上拍摄的照片
文章目的
本文的主要目的是向读者介绍一种最简单也是最常用的工具，使用 PuLP 库在 Python 中编写线性优化问题的代码。它还对优化和线性编程进行了快速介绍，这样，即使那些对优化、规范分析或运筹学知之甚少或一无所知的读者也可以很容易地理解文章的上下文以及文章将要讨论的内容。我们也将触及如何用数学符号来表达一个 LP。
照片由digity Marketing在 Unsplash 上拍摄
先决条件
给定的前提条件有是好的，没有必要。
- Python 中的基本编码知识。
- 对线性规划、目标函数、约束和决策变量有基本的了解。
照片由hello queue在 Unsplash 上拍摄
优化简介
优化是通过控制约束环境中的一组决策来寻找给定目标的最大值或最小值的过程。简而言之，优化问题包括通过从允许的集合中系统地选择输入值并计算函数值来最大化或最小化真实函数。
真实函数(目标函数)可以是从仓库向客户交付货物的成本，在给定有限数量的驾驶员和时间(约束)的情况下，我们希望通过选择最优路线和最优车辆组(决策变量)来最小化该成本。这是运筹学和最优化领域中路线最优化的一般情况。
计算机科学领域的另一个非常著名的问题是 TSP 或旅行推销员问题，其中我们希望找到最短的路线或最便宜的路线来穿越所有城市，给定它们之间的成对距离。在这种情况下，我们的目标函数变成最小化旅行的总距离(或总成本),决策变量变成二元变量，告知旅行者是否应该从城市 I 旅行到城市 j，并且应用约束，使得旅行者覆盖所有城市并且不会两次访问一个城市。我们今天还将处理一个更简单但类似的问题。
斯蒂芬·门罗在 Unsplash 上拍摄的照片
线性规划导论
线性规划基本上是最优化的一个子集。线性规划或线性优化是一种优化技术，其中我们试图使用一组变化的决策变量为线性约束系统找到线性目标函数的最优值。
由 Inductiveload —自有作品，公共领域，【https://commons.wikimedia.org/w/index.php?curid=6666051
理解手头的问题
您希望将货物从两个不同的仓库运送到四个不同的客户手中的成本降到最低。每个仓库的供应量是有限的，每个客户都有一定的需求。我们需要通过从给定的仓库装运产品来满足客户的需求，这样运输的总成本是最低的，并且我们还能够使用每个仓库可用的有限供应来满足客户的需求。
数据
假设这家公司是只供应鞋类的 crocs，这里的客户是它的经销商，他们需要大量的 Crocs。供应的产品在性质上是一致的。
- 仓库 i 到客户 j 的运输矩阵成本如下。每个值也可以表示为 Cij，表示从仓库 I 运送到客户 j 的成本 C。
成本矩阵
2.客户需求和仓库可用性如下。
需求(必需)和供应(可用性)矩阵
公式化问题
让我们开始用数学方程式来表述这个问题。我们需要确定我们的 LP 的 3 个主要组成部分，即
1)决策变量
我们将我们的决策变量定义为 Xij，它基本上告诉我们 X 个产品应该从仓库 I 交付给客户 j。
决策变量定义
2)目标函数
我们的目标函数被定义为运输这些产品的总成本，我们需要最小化这个总成本。因此，目标函数定义为:-
目标函数
3)制约因素
对于给定的问题，我们将有两种主要类型的约束:-
**3.1)仓库约束条件或供应约束条件:**这些约束条件基本上表示每个仓库在所有 4 个客户之间完成的总供应量小于或等于该仓库的最大可用性/容量。
供应限制
**3.2)客户约束或需求约束:**这些约束基本上是说，对于每个客户，跨两个仓库完成的供应应该等于(或大于等于)该客户的需求。我们可以用≥代替=因为我们的目标函数总是试图最小化成本，因此永远不会提供超过需要的。这样做是因为在一些优化问题中，我们可能无法通过严格的等式约束得到可行的解决方案。虽然，这里的情况并非如此。
需求约束
问题的表述到此结束。我们现在进一步理解如何用 Python 来编码这个问题，并找到供应货物的最小成本。我们还将获得最佳答案，该答案将建议应该由哪个仓库向哪个客户供应多少货物。
纸浆快速入门
PuLP 是一个用 Python 编写的免费开源软件。它用于将优化问题描述为数学模型。然后，PuLP 可以调用众多外部 LP 解算器(CBC、GLPK、CPLEX、Gurobi 等)中的任何一个来求解该模型，然后使用 python 命令来操作和显示解决方案。请务必阅读它的文档，它非常有用。以下链接也有助于您了解如何在 Python 环境中安装库 PuLP 和任何所需的求解器。
来源:https://coin-or . github . io/pulp/main/installing _ pulp _ at _ home . htm
钻研代码
你可以在下面的 Github repo 中找到将在下面解释的完整代码(Jupyter notebook)。
使用 Python (PuLP)进行线性编程的快速指南。-mni PS/线性编程-Python-1
github.com](https://github.com/mnips/Linear-Programming-Python-1)
导入所需的库
from pulp import *
import pandas as pd
import numpy as np
第一条语句从 PuLP 库中导入我们将使用的所有必需的函数。Pandas 是一个数据操作库，Numpy 是一个主要用于在 Python 中处理多维数组的库。我通常只导入这些库，因为它们在几乎所有的数据分析项目中都会用到。
数据定义
让我们定义数据并将其分配给变量，然后这些变量可用于输入模型、目标函数和约束条件。
所有变量都是直观的，易于解释。
模型初始化
我们可以通过调用LpProblem()函数来初始化模型。函数中的第一个参数表示我们想要给模型起的名字。第二个参数告诉我们的模型，我们是否希望最小化或最大化我们的目标函数。
model = LpProblem("Supply-Demand-Problem", LpMinimize)
如果你想最大化你的目标函数，你可以使用LpMaximize。
定义决策变量
您可以在您的模型中定义变量名，以使您的模型对于稍后阅读它的人来说看起来更直观。因此，我们为我们的决策变量创建指数，这将在后面定义。
variable_names = [str(i)+str(j) for j in range(1, n_customers+1) for i in range(1, n_warehouses+1)]variable_names.sort()print("Variable Indices:", variable_names)
输出如下所示:-
此外，我们使用LpVariables.matrix来定义变量。在定义决策变量时，我们也可以使用字典或单例变量，但在这种情况下，这似乎是最好的方法，因为对于更大的问题，仓库或客户的数量可能会增加。我们将决策变量命名为X，并使用上面定义的索引作为第二个参数，这有助于 PuLP 理解我们想要一个 2*4 的矩阵。第三个参数是一个类别，它告诉我们决策变量只能取Integer值。默认为Continuous。这也告诉我们，我们的线性规划问题其实是一个整数 LP 。如果我们还有可以取连续值的决策变量，我们称之为 MILP 或混合整数 LP 。在第四个也是最后一个参数中，我们设置了一个0的lower bound，表示我们的决策变量≥ 0。为了利用 Numpy 数组操作，我们可以将决策变量转换为 Numpy 数组。
DV_variables = LpVariable.matrix("X", variable_names, cat = "Integer", lowBound= 0 )allocation = np.array(DV_variables).reshape(2,4)print("Decision Variable/Allocation Matrix: ")print(allocation)
输出如下所示:-
目标函数
现在让我们定义我们的目标函数，它基本上是提供产品的总成本。换句话说，它是成本矩阵和上面定义的分配矩阵的和积。我们可以如下定义我们的目标函数。lpSum在 Python 中与sum函数交替使用，因为它在执行与 PuLP 变量相关的操作时速度更快，并且还能很好地汇总变量。我们使用+=速记操作符进一步将目标函数添加到模型中。
obj_func = lpSum(allocation*cost_matrix)print(obj_func)model += obj_funcprint(model)
让我们看看我们的模型现在是什么样子。
正如我们所看到的，我们已经给我们的问题起了一个名字。在目标函数中，我们试图最小化成本，并且我们所有的决策变量都就位了。最好在创建模型时打印出来，以便了解我们是否遗漏了什么。现在我们继续向我们的模型添加约束。
限制
我们需要添加两种主要类型的约束:-
1)仓库或供应限制
如前所述，这些约束表明，对于一个给定的仓库或i-th仓库，在所有客户之间完成的总分配或供应的产品应该不违反该仓库的可用性。
此约束的输出:
2)客户或需求约束
这些约束表明，为每个客户或j-th客户完成的分配应该满足该客户的需求。
此约束的输出:
检查模型
现在我们已经完成了所有需要的公式，让我们检查模型看起来怎么样。这可以通过打印模型来实现:print(model)。
完全优化模型
我们还可以将这个模型保存在一个.lp文件中，任何不熟悉我们模型的人都可以参考这个文件。它基本上就像一个文本文件，包含上面打印的优化模型的确切细节。
model.writeLP(“Supply_demand_prob.lp”)
运行模型并检查状态
既然我们已经检查了模型看起来很好，我们现在应该运行模型并检查我们是否得到了问题的可行/最优解决方案。默认情况下，纸浆使用 CBC 解算器，但我们也可以启动其他解算器，如 GLPK，古罗比等。我已经明确地在这里调用了 CBC。类似地，我们可以调用任何其他求解器来代替 CBC。
上面代码的输出是Optimal，它告诉我们我们的模型已经能够找到问题的最优解。
输出目标函数值和决策变量值
现在让我们看看公司通过打印出我们的问题的最优解(即目标函数值)而必须承担的最小成本，并看看将产品从仓库运送到客户的最优安排。
输出是:
此外，我们可以检查每个仓库需要供应多少产品，以及每个仓库需要多少容量。
输出:
因此，我们在 2 号仓库只需要 45000 件，而不是 80000 件。虽然在这种情况下非常幼稚，但是我们可以从优化问题的输出中做很多类似的分析，做出相关的商业决策。比如，如果每个仓库都有运营成本。这个供求问题可以有许多不同的形式。探索！
摘要
我们简要地看了最优化和线性规划。我们还学习了如何用数学方程来表述一个问题。此外，我们通过利用 Python 和 PuLP 库并分析其结果，深入研究了编码 LP 问题。就这样，我们到了这篇文章的结尾。希望你觉得这有用！
我还要感谢我亲爱的朋友 Karan Bhanot，他通过他的文章激励了我，也激励我与世界分享我的知识！
注: 我用过 Python 3 . 7 . 6 版和纸浆 2.1 版。您没有必要使用相同的版本，但有时由于 PuLP 库中的一些更新，可能会有微小的差异导致错误(主要是由于语法变化)，因此将此作为快速注释添加进来。
用 Python 进行线性编程
面向工业工程师的 Python
探索 SciPy 的“linprog”功能
由 Kim 拍摄的图片可在 Unsplash 获得
运筹学
运筹学是一种科学的决策方法，通常在需要分配稀缺资源的条件下，寻求系统的最佳设计和操作。制定决策的科学方法需要使用一个或多个数学/优化模型(即实际情况的表示)来做出最佳决策。
优化模型寻求在满足给定 约束 的决策变量的所有值的集合中找到优化(最大化或最小化) 目标函数 的 决策变量 的值。它的三个主要组成部分是:
- 目标函数:要优化的函数(最大化或最小化)
- 决策变量:影响系统性能的可控变量
- 约束:决策变量的一组限制(即线性不等式或等式)。非负约束限制决策变量取正值(例如，您不能产生负的项目数 x 1、 x 2 和 x 3)。
优化模型的解称为最优可行解。
建模步骤
精确地模拟一个运筹学问题是最重要的——有时也是最困难的——任务。一个错误的模型将导致一个错误的解决方案，因此，不会解决原来的问题。以下步骤应该由具有不同专业领域的不同团队成员执行，以获得模型的准确和更好的视图:
- 问题定义:定义项目的范围，确定结果是三个要素的确定:决策变量的描述、目标的确定和限制条件(即约束条件)的确定。
- 模型构建:将问题定义转化为数学关系。
- 模型求解:使用标准优化算法。在获得解决方案后，应进行敏感性分析，以找出由于某些参数的变化而导致的解决方案的行为。
- 模型有效性:检查模型是否如预期的那样工作。
- 实施:将模型和结果转化为解决方案的建议。
线性规划
线性规划(也称为 LP)是一种运筹学技术，当所有目标和约束都是线性的(在变量中)并且所有决策变量都是连续的时使用。在层次结构中，线性规划可以被认为是最简单的运筹学技术。
Python 的 SciPy 库包含了 linprog 函数来解决线性编程问题。使用 linprog 时，在编写代码时需要考虑两个因素:
- 该问题必须被公式化为最小化问题
- 不等式必须表示为≤
最小化问题
让我们考虑下面要解决的最小化问题:
让我们来看看 Python 代码:
结果
最大化问题
由于 Python 的 SciPy 库中的 linprog 函数被编程为解决最小化问题，因此有必要对原始目标函数进行转换。通过将目标函数的系数乘以-1(即通过改变它们的符号)，可以将每个最小化问题转化为最大化问题。
让我们考虑下面要解决的最大化问题:
让我们来看看 Python 代码:
结果
总结想法
线性规划代表了一个更好的决策制定的伟大的优化技术。Python 的 SciPy 库中的linprog函数允许用几行代码解决线性编程问题。虽然有其他免费优化软件(如 GAMS、AMPL、TORA、LINDO)，但使用 linprog 函数可以节省大量时间，因为您不必从头开始编写单纯形算法，也不必检查每个操作直到达到最优值。
— —
如果你觉得这篇文章有用，欢迎在 GitHub 上下载我的个人代码。你也可以直接在 rsalaza4@binghamton.edu 给我发邮件，在LinkedIn上找到我。有兴趣了解工程领域的数据分析、数据科学和机器学习应用的更多信息吗？通过访问我的媒体 个人资料 来探索我以前的文章。感谢阅读。
罗伯特
线性回归 101
学习线性回归的外行指南
Alexander Schimmeck 在 Unsplash 上的照片
在这篇文章中，我想谈谈最简单的机器学习算法之一——线性回归。我将尝试使用一个真实的例子，用简单的术语解释这个模型。
让我们深入研究一下。
介绍
让我们假设你在学校参加一个高中家长会。你看到学校里的孩子们都和他们的父母一起来，坐在一个巨大的大厅里。
当孩子们和他们的父母一起起床时，你会注意到一些奇怪的事情。孩子们的身高和他们父母的身高直接相关。
即如果孩子的父母是高个子，那么孩子就是高个子，如果父母是矮个子，那么孩子也是矮个子。
你决定深入研究这个问题，从收集一些关于父亲和孩子身高的数据开始。
您将获得以下信息。
形象化
现在，您开始考虑以图表形式显示这些数据的最佳方式。
首先绘制一个散点图，以了解这两组身高之间的关系，Y 轴表示孩子的身高，X 轴表示父亲的身高
描述孩子和父亲身高之间关系的情节
画出图后，你可以直观地看到父亲的身高和孩子的身高之间有明确的关系。
随着父亲身高的增加，孩子的身高也会增加。
接下来，我们将尝试通过散点图拟合一条线，使点和线之间的距离最小化。直线和每个点之间的距离称为误差。
你现在考虑如何用方程的形式表达这条线。我们可以这样写一个方程。
孩子的身高=常数+直线的斜率*(父亲的身高)+误差项
上述方程被称为简单线性回归方程。
我们可以通过延长这条线，并查看它与 Y 轴相交的位置，来获得常数的值。
常数项和直线斜率也称为该方程的参数，称为β0 和β1。
所以我们可以把上面的等式写成
y = β0 + β1 * x + e
其中 y 是孩子的身高，x 是父亲的身高，β0 是常数，β1 是斜率，e 是误差项。
如果您有足够的关于父亲身高和孩子身高的数据(训练数据),我们可以使用它来创建一个简单的线性模型，如果我们知道父亲的身高，就可以用它来预测孩子的身高。
模型解释
任何建模活动的一个重要部分是如何解释模型参数。
我们可以这样解释这个模型。
β1 —该参数告诉我们，如果父亲的身高增加 1 个单位，孩子的身高会增加多少，反之亦然。
β0 —这可以解释为当父亲的身高为零时孩子的身高，这在现实世界中不可能发生，因为父亲的身高永远不会为零。
总结
简单线性回归虽然是一个非常幼稚的模型，因为它只假设线性关系，但仍然非常重要，因为它有助于我们理解因变量和自变量之间的关系。
我们现在还可以在这个模型的基础上添加更多的特征，比如母亲的身高、孩子的性别等等。
这个添加了多个独立特征的模型被称为多元线性回归，我将在以后的帖子中讨论。
机器学习中线性回归的一种实用方法
线性回归初学者实践指南
在之前的博客文章中，我试图给你一些关于机器学习基础的直觉。在本文中，我们将从我们的第一个机器学习算法开始，那就是线性回归。
首先，我们将学习线性回归的数学方面，然后我将尝试阐明一些重要的回归术语，如假设和成本函数，最后我们将通过构建我们自己的回归模型来实施我们所学的内容。
什么是线性回归？
回归模型是监督学习模型，通常在要预测的值具有离散或定量性质时使用。使用回归模型的一个最常见的例子是通过训练该地区的房屋销售数据来预测房屋价格。
图片来自 CDOT 维基
线性回归模型背后的思想是获得最适合数据的直线。所谓最佳拟合，意思是所有点离回归线的总距离应该是最小的。通常这些点离我们回归线的距离被称为误差，尽管从技术上来说并不是。我们知道直线方程的形式是:
其中 y 是因变量，x 是自变量，m 是直线的斜率，c 是系数(或 y 截距)。这里，y 被认为是因变量，因为它的值取决于自变量和其他参数的值。
该方程是任何线性回归问题的基础，并被称为线性回归的假设函数。大多数机器学习算法的目标是构建一个模型，即一个假设，以根据我们的自变量估计因变量。
这个假设将我们的输入映射到输出*。*线性回归的假设通常表示为:
在上面提到的表达式中， hθ(x) 是我们的假设，θ0 是截距，而 θ1 是模型的系数。
了解成本函数
成本函数用于计算模型的执行情况。通俗地说，代价函数就是所有误差的总和。在构建我们的 ML 模型时，我们的目标是最小化成本函数。
图片来自JMP.com
回归问题中经常使用的一个常见函数是均方误差或 MSE ，它测量已知值和预测值之间的差异。
事实证明，取上述方程的根是更好的选择，因为这些值不太复杂，因此通常使用均方根误差或 RMSE 。我们还可以使用平均绝对误差等其他参数来评估回归模型。
RMSE 告诉我们数据点离回归线有多近。现在，我们将通过构建我们自己的线性回归模型来预测房子的价格，从而实现我们到目前为止所学到的知识。
来自介质的图像
在这里，我将使用 Google Colab 来构建这个模型。您还可以使用 Jupyter 笔记本之类的其他 ide 来玩这个模型。
用于这个线性回归项目的代码可以在 这里 找到。
第一步:导入库并加载数据
我们的第一步是导入构建模型可能需要的库。没有必要在一个地方导入所有的库。为了开始，我们将进口熊猫， Numpy ， Matplotlib 等。
*#Import Required Libraries* import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns*#Read the Dataset*df=pd.read_csv('kc_house_data.csv')
一旦这些库被导入，我们的下一步将是获取数据集和加载我们的数据。在加载数据时，需要注意的是文件名应该有它的格式(。csv/。xls)在末尾指定。
我用于这个模型的数据可以直接从 这里 下载。
CSV 文件最常用于此目的，尽管 excel 表格也可以使用。唯一的区别是，在使用 excel 表作为数据集时，我们将不得不使用 read_excel() ，而不是 read_csv() 。
第二步:可视化数据
成功加载数据后，我们的下一步是可视化这些数据。数据可视化是数据科学家角色的重要组成部分。建议将数据可视化，以便找出不同参数之间的相关性。
*#Visualising the data using heatmap* plt.figure()
sns.heatmap(df.corr(),cmap='coolwarm')
plt.show()
Matplotlib 和search是优秀的库，可以用来可视化我们在各种不同地块上的数据。
第三步:特征工程
在可视化我们的数据的同时，我们发现两个参数之间有很强的相关性: sqft_living 和 price 。因此，我们将使用这些参数来构建我们的模型。
*#Selecting the required parameters* area = df[‘sqft_living’]
price = df['price']x = np.array(area).reshape(-1,1)
y = np.array(price)
更多的参数也可以添加到模型中，尽管这可能会影响其准确性。使用各种特征来预测响应变量结果的模型被称为多元回归模型。
第四步:拟合线性回归模型
选择所需参数后，下一步是从 sklearn 库中导入方法 train_test_split 。这用于将我们的数据分为训练和测试数据。通常 70–80%的数据作为训练数据集，而剩余的数据构成测试数据集。
*#Import LinearRegression and split the data into training and testing dataset* from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(x,y,test_size=0.2,random_state = 0)y_train = y_train.reshape(-1,1)
y_test = y_test.reshape(-1,1)*#Fit the model over the training dataset* from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X_train,y_train)
在此之后，从sk learn . model _ selection导入线性回归，并且模型适合训练数据集。我们模型的截距和系数可以计算如下:
*#Calculate intercept and coefficient* print(model.intercept_)
print(model.coef_)pred=model.predict(X_test)
predictions = pred.reshape(-1,1)*#Calculate root mean squared error to evaluate model performance* from sklearn.metrics import mean_squared_error
print('MSE : ', mean_squared_error(y_test,predictions)
print('RMSE : ', np.sqrt(mean_squared_error(y_test,predictions)))
可以通过找到模型的均方根误差来评估模型的性能。RMSE 越小，模型越好。
使用梯度下降的线性回归
梯度下降是一种迭代优化算法寻找一个函数的最小值。为了理解这个算法，想象一个没有方向感的人想要到达谷底。
他走下斜坡，在斜坡陡的时候迈大步，在斜坡不那么陡的时候迈小步。他根据当前位置决定下一个位置，当他到达他的目标山谷底部时停下来。梯度下降也是如此。
*#Initializing the variables* m = 0
c = 0
L = 0.001
epochs = 100n = float(len(x))
梯度下降法逐步应用于我们的 m 和 c。最初让 m = 0，c = 0。设 L 为我们的学习率。这控制了 m 值随每步变化的程度。
for i in range(epochs):
Y_pred=m*x+c
Dm = (-2/n)*sum(x*(y-Y_pred))
Dc = (-2/n)*sum(y-Y_pred)
m = m-L*Dm
c = c-L*Dc
print(m,c)*#Predicting the values* y_pred = df['sqft_living'].apply(lambda a:c+m*a)
y_pred.head()
优选地，L 被赋予一个小的值，以便提高精度。我们的下一步是计算损失函数相对于 m 和 c 的偏导数。一旦完成，我们更新 c 和 m 的值，并重复该过程，直到我们的损失函数非常小。
至此，我们已经到了这篇文章的结尾。我希望这篇文章能帮助你了解线性回归算法背后的思想。如果你有任何问题，或者如果你认为我犯了任何错误，请联系我！可以通过:邮箱或 LinkedIn 与我联系。
线性回归算法——面向非数学家的隐蔽数学
内部 AI
它是机器学习中借用的数学和统计学领域最古老的算法之一。
在计算机出现之前，线性回归是在不同领域中使用的最流行的算法之一。今天有了强大的计算机，我们可以解决多维线性回归，这在以前是不可能的。在一元或多维线性回归中，基本的数学概念是相同的。
今天有了机器学习库，像 Scikit -learn ，就可以在建模中使用线性回归，而不需要理解它背后的数学概念。在我看来，对于一个数据科学家和机器学习专业人士来说，在使用算法之前，了解算法背后的数学概念和逻辑是非常必要的。
我们大多数人可能没有学习过高等数学和统计学，看到算法背后的数学符号和术语时，我们会感到害怕。在本文中，我将用简化的 python 代码和简单的数学来解释线性回归背后的数学和逻辑，以帮助您理解
概述
我们将从一个简单的一元线性方程开始，没有任何截距/偏差。首先，我们将学习像 Scikit-learn 这样的包所采取的逐步解决线性回归的方法。在本演练中，我们将理解梯度下降的重要概念。此外，我们将看到一个带有一个变量和截距/偏差的简单线性方程的示例。
步骤 1: 我们将使用 python 包 NumPy 来处理样本数据集，并使用 Matplotlib 来绘制各种可视化图形。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
第二步:让我们考虑一个简单的场景，其中单个输入/自变量控制结果/因变量的值。在下面的代码中，我们声明了两个 NumPy 数组来保存自变量和因变量的值。
Independent_Variable=np.array([1,2,3,12,15,17,20,21,5,7,9,10,3,12,15,17,20,7])
Dependent_Variable=np.array([7,14,21,84,105,116.1,139,144.15,32.6,50.1,65.4,75.4,20.8,83.4,103.15,110.9,136.6,48.7])
*第三步:*让我们快速画一个散点图来了解数据点。
plt.scatter(Independent_Variable, Dependent_Variable, color='green')
plt.xlabel('Independent Variable/Input Parameter')
plt.ylabel('Dependent Variable/ Output Parameter')
plt.show()
我们的目标是制定一个线性方程，它可以预测自变量/输入变量的因变量值，误差最小。
因变量=常数*自变量
用数学术语来说，Y =常数*X
就可视化而言，我们需要找到最佳拟合线，以使点的误差最小。
在机器学习领域，最小误差也被称为损失函数。
损失函数公式(作者用 word 写的然后截图)
我们可以用所有独立的数据点计算方程 Y =常数*X 中每个假定常数值的迭代损失。目标是找到损耗最小的常数，并建立方程。请注意，在损失函数等式中,“m”代表点数。在当前示例中，我们有 18 个点，因此 1/2m 转化为 1/36。不要被损失函数公式吓到。我们将损失计算为每个数据点的计算值和实际值之差的平方和，然后除以两倍的点数。我们将在下面的文章中借助 python 中的代码一步一步地破译它。
步骤 4: 为了理解识别方程背后的核心思想和数学，我们将考虑下面代码中提到的有限的一组常数值，并计算损失函数。
在实际的线性回归算法中，特别是间隙，常数被考虑用于损失函数计算。最初，考虑用于损失函数计算的两个常数之间的差距较大。随着我们越来越接近实际的解决方案常数，考虑更小的差距。在机器学习的世界中，学习率是损失函数计算中常数增加/减少的差距。
m=[-5,-3,-1,1,3,5,6.6,7,8.5,9,11,13,15]
步骤 5: 在下面的代码中，我们计算所有输入和输出数据点的每个常量值(即在前面步骤中声明的列表 m 中的值)的损失函数。
我们将每个常数的计算损失存储在一个 Numpy 数组“errormargin”中。
errormargin=np.array([])
for slope in m:
counter=0
sumerror=0
cost=sumerror/10
for x in Independent_Variable:
yhat=slope*x
error=(yhat-Dependent_Variable[counter])*(yhat-Dependent_Variable[counter])
sumerror=error+sumerror
counter=counter+1
cost=sumerror/18
errormargin=np.append(errormargin,cost)
*第六步:*我们将绘制常量的计算损失函数，以确定实际常量值。
plt.plot(m,errormargin)
plt.xlabel("Slope Values")
plt.ylabel("Loss Function")
plt.show()
曲线在最低点的常数的值是真正的常数，我们可以用它来建立直线的方程。
在我们的例子中，对于常数 6.8 的值，曲线处于最低点。
该值为 Y=6.8*X 的直线可以以最小的误差最好地拟合数据点。
这种绘制损失函数并在损失曲线的最低点识别方程中固定参数的真实值的方法被称为 梯度下降 。作为一个例子，为了简单起见，我们考虑了一个变量，因此损失函数是一个二维曲线。在多元线性回归的情况下，梯度下降曲线将是多维的。
我们已经学习了计算自变量系数的内功。接下来，让我们一步一步地学习线性回归中计算系数和截距/偏差的方法。
*第一步:*和前面一样，让我们考虑一组自变量和因变量的样本值。这些是可用的输入和输出数据点。我们的目标是制定一个线性方程，它可以预测自变量/输入变量的因变量值，误差最小。
因变量=(系数*自变量)+常数
用数学术语来说，y=(系数*x)+ c
请注意，系数也是一个常数项乘以方程中的自变量。
Independent_Variable=np.array([1,2,4,3,5])
Dependent_Variable=np.array([1,3,3,2,5])
*第二步:*我们将假设系数的初始值和常数“m”和“c”分别为零。我们将在误差计算的每一次迭代之后以 0.001 的小学习率增加 m 和 c 的值。Epoch 是我们希望在整个可用数据点上进行这种计算的次数。随着历元数量的增加，解会变得更加精确，但这会消耗时间和计算能力。基于业务案例，我们可以决定计算值中可接受的误差，以停止迭代。
LR=0.001
m=0
c=0
epoch=0
步骤 3: 在下面的代码中，我们在可用的数据集上运行 1100 次迭代，并计算系数和常数值。
对于每个独立的数据点，我们计算相关值(即 yhat)，然后计算计算的和实际的相关值之间的误差。
基于该误差，我们改变系数和常数的值用于下一次迭代计算。
新系数=当前系数—(学习率*误差)
新常数=当前常数-(学习率误差独立变量值)
while epoch<1100:
epoch=epoch+1
counter=0
for x in Independent_Variable:
yhat=(m*x)+c
error=yhat-Dependent_Variable[counter]
c=c-(LR*error)
m=m-(LR*error*x)
counter=counter+1
在对可用数据集进行 1100 次迭代后，我们检查系数和常数的值。
print("The final value of m", m)
print("The final value of c", c)
数学上可以表示为 y=(0.81*x)+0.33
最后，让我们将之前的输出与 Scikit -learn 线性回归算法的结果进行比较
from sklearn.linear_model import LinearRegression
reg = LinearRegression().fit(Independent_Variable.reshape(-1,1), Dependent_Variable)
print(reg.coef_)
print(reg.intercept_)
通过对可用数据集进行 1100 次迭代，系数和常数/偏差的计算值非常接近于 Scikit -learn 线性回归算法的输出。
我希望这篇文章能让你对线性回归的幕后数学计算和概念有一个明确的理解。此外，我们已经看到了梯度下降被应用于寻找最优解的方式。在多元线性回归的情况下，数学和逻辑保持不变，只是在更多的维度上进一步扩展。
Python 中房价的线性回归分析
如今，我们都变得如此习惯于谈论深度学习、大数据、神经网络……我们似乎忘记了，即使这些大话题正在蓬勃发展，但不是每个企业都需要它们，至少现在是这样。因此，我想分享一点我在常见统计话题上的经验。在这篇和下一篇博客中，我将分别演示如何用 Python 和 R 进行线性回归分析。
项目摘要
这个由 Kaggle 提供的项目包括一个由 79 个解释变量组成的数据集，描述了爱荷华州 Ames 住宅的各个方面。你可以在这里找到项目的所有数据来源和描述。
这项竞赛的目标是预测每所挂牌房屋的最终销售价格。我们做这个问题是因为我们希望应用我们的回归技术和数据探索性分析技巧。这个习题集允许我们使用以上所有的练习。此外，数据集包含许多缺失值，这使我们能够获得处理缺失数据的经验。
我们运行了多重回归技术，如 XGB 回归器、 SGD 回归器、 MLP 回归器、决策树回归器、随机森林回归器、 CatBoost 回归器、轻型 GBM 和 SVR 。我们把一些放在一起，实验看哪一个的均方根误差最小。
我们提交的预测有两列(图 1.1)。第一列是房子的 ID，第二列是我们预测的销售价格。根据预测值对数和观察到的销售价格日志之间的均方根误差对预测进行评分。最小的 RMSE 是最好的预测。在运行基础模型和优化模型后，我们发现我们的 tubed CatBoost 回归模型为我们获得了最低的 RMSE 分数 0.12392，在我们提交时，它在排行榜上排名第 863 位(共 4525 位)。
图 1.1 —提交样本
数据
-概述
竞赛为我们提供了两个数据集。其中之一是有 1460 个观察值和 81 列的训练数据，它包含每栋房子的 ID 和销售价格。另一个数据集是维持文件，它包含 1459 个观察值。
我们首先研究了训练数据集的特征类型。图 2.1 显示整个数据集有 43 个分类特征和 36 个数字变量，包括“ID”。
图 2.1 —数字的数量&分类特征
然后，我们努力检测数据集中是否有任何丢失的值。在训练数据集中，有 19 个变量缺少值。维持数据集包含 33 个缺失值的要素。有些要素缺失值的百分比非常高。例如，PoolQC(池质量)特性的缺失值高达 99.5%。像 PoolQC 这样的特性可能会损害预测模型的准确性，并导致我们得出无效的结论。因此，我们删除了至少有 80%的值丢失的特征。详情见“数据处理”部分。
在我们开始处理数据之前，我们还研究了 38 个数字特征的必要统计信息。在图 2.2 中，我们发现一些特性的范围相对较大，比如销售价格。我们将要使用的许多回归模型需要我们在处理之前对这些特征进行某种类型的转换。
图 2.2 —所有数字变量的统计信息
-数据处理
缺失值总是会损害预测模型的准确性。在进行任何转换之前，我们决定首先检查至少有 80%的值丢失的列。根据图 2.3，有四个变量有大量的缺失值，它们是“Alley”、“PoolQC”、“Fence”和“MiscFeature”。
图 2.3——至少有 80%缺失值的变量
在处理数据时，确保对定型数据集和维持数据集进行相同的更改是至关重要的。每个数据集中列数的差异、不一致的数据格式以及分类变量中不匹配的值数目都可能给我们的模型带来麻烦。因此，为了更好地准备数据并确保任何转换都将反映在定型数据集和维持数据集上，我们使用图 2.4 中的代码将这两个数据集临时合并为标签 0 和 1，并且还提供了将数据集拆分回原始数据集的代码。
图 2.4 —合并和取消合并数据集的代码
对于所有的分类变量，我们将它们转换成虚拟变量来处理缺失值。我们用数字变量中所有缺失值的当前变量的平均值来填充它们。
-特征工程
我们还在现有变量的基础上创建了一些新变量。例如，我们使用关于房子何时装修和出售的数据，生成了一个名为“AgeofHouse”的新变量，这有助于我们直观地了解房价和房子质量之间的关系。
图 2.5 —新功能“房屋时代”
我们想知道是否需要对“销售价格”进行对数转换，因为它包含一些相对较大的值，这些值会对模型的准确性产生负面影响。对数变换前后的残差与拟合图(图 2.6 和图 2.7)证实了我们的假设，即变换使模型更加线性。
图 2.6——对数变换前的残差与拟合值
图 2.7——对数变换后的残差与拟合值
-相关图
我们创建了一些图表来研究数据。我们首先查看目标变量和其余变量之间的相关性。图 2.8 显示了与房屋销售价格绝对相关的前 30 个特征。深色表示变量与销售价格负相关。
图 2.8 —前 30 个“最”相关变量
我们还详细研究了“销售价格”和个体变量之间的关系。我们选择了几个有趣的图表来呈现在这份报告中。例如，我们制作了一个散点图(图 2.9)关于“销售价格”和“房屋年龄”，这是我们在上一步中创建的特征。我们可以看到，房子越新，房子越有可能以更高的价格出售。因此，我们确信价格超过 50 万美元的房屋是在最近 15 年内建造或翻新的。
图 2.9——房龄与销售价格
另一个例子是“总体平等”和“销售价格”的情节。“总体质量”是与“销售价格”最正相关的变量，该图(图 2.10)清楚地表明，评级越高，价格越高。
图 2.10 —总体质量与销售价格
我们的最终训练数据集有 1460 个观察值和 289 个变量，而维持数据集有 1459 个观察值和 289 个变量。
建模
-单个型号的性能(基本型号与调整型号)
在这种情况下，选择了九个回归模型:
OLS，XGBRegressor，SGD regressor
DecisionTreeRegressor，RandomForestRegressor，SVR
CatBoostRegressor，LightGBM，MLPRegressor 。
我们使用 OLS 作为基础模型，并生成了图 3.1 中模型的回归结果:
#original
#generate OLS model
sm_model = sm.OLS(y_train, sm.add_constant(X_train))
sm_model_fit = sm_model.fit()print(sm_model_fit.summary())
图 3.1 — OLS 回归结果
根据 OLS 报告，我们了解到 R 平方为 0.943，这意味着该模型可以解释 94.3%的“销售价格”。该模型的 kaggle 得分为 0.15569。
下表(图 3.2)列出了所有模型的 kaggle 分数，包括未调整和调整的参数。在所有其他模型中，具有调谐参数的 CatBoostRegressor 具有最佳性能， 0.12392 。在我们撰写本报告时，这个分数使我们在 4525 个团队中排名第 863 位(相当于前 19%)。
# Example codemlp=MLPRegressor(random_state=42)#paramters for gridsearch
mlp_param = {
'hidden_layer_sizes': [(273,230,30), (273,230,20), (273,230,50)],
'activation': ['tanh'],
'solver': ['sgd'], #, 'adam'
'alpha': [0.0001],
'learning_rate': ['adaptive'], #'constant',
}
----------------------#Apply GridSearchCV
mlp_grid = GridSearchCV(mlp,
mlp_param,
cv = 5,
n_jobs = 3,
verbose=True)mlp_grid.fit(X_train,y_train)
print(mlp_grid.best_score_)
print(mlp_grid.best_params_)output:
0.8218847474014372
{'activation': 'tanh', 'alpha': 0.0001, 'hidden_layer_sizes': (273, 230, 30), 'learning_rate': 'adaptive', 'solver': 'sgd'}
----------------------#Fit the model with bets paramters
mlp_best_model=MLPRegressor(activation= 'tanh',
alpha=0.0001,
hidden_layer_sizes=(273, 230, 30),
learning_rate= 'adaptive',
solver= 'sgd')
mlp_best_model.fit(X_train,y_train)
----------------------#Make Predictions and submit files
mlp_prediction = mlp_best_model.predict(test.drop(["Id"],1))
mlp_DATA = pd.DataFrame({"ID":holdout["Id"],
"SalePrice":np.exp(mlp_prediction)})
mlp_DATA.to_csv("mlpstandard1.csv",index=False)
图 3.2——卡格尔得分表
-堆叠模型和平均系综
我们还创建了几个堆积模型和平均系综模型。这些类型的模型通常可以通过让多个模型一起工作来帮助减少误差。然而，这些模型都不能胜过 CatboostRegressor。输出如图 3.3 所示:
图 3.3——叠加和平均系综模型的 Kaggle 得分
结论
总的来说，我们在“房价高级回归技术”问题上的主要挑战是有大量的缺失数据。我们用多种解决方案测试了数据的缺失值，但是仍然很难找到一种方法来显著提高模型的准确性。此外，我们认为，如果我们可以有一个更大的训练数据集，它也可能有助于改进模型。
在提交了各种模型之后，我们注意到我们调整的 CatBoost 回归模型是房价的最佳预测器。我们还运行了由三个或更多调整的回归模型组成的四个堆叠和平均集成模型，但没有一个比 CatBoost 回归器得分更高。
我们发现这个项目很有趣，因为我们使用许多不同的解释变量来预测单一价格。这也很有趣，因为我们可以在其他城市或州的数据集上使用相同的方法来预测房价，并比较不同位置的平均价格。我们也可以使用这样的预测模型来帮助预测未来的房屋市场销售价格。
本博客来源于蒋金航和马人·奥尔森的班级项目报告。感谢安德鲁阿比莱斯和张子萱的贡献。最初的报告发布在我的网站上，我写这篇博客时做了一些小改动。你可以在这里找到建模的代码脚本。
以前的文章:
线性回归及其假设
来源:stocksnap.io
举例说明了线性回归的假设，如共线性、多元正态性、自相关、同异方差
上周，我在帮助我的朋友准备一个数据科学家职位的面试。我在网上搜索了相关问题，发现问题“数据科学涉及的假设有哪些？“T1”在我的搜索中出现得非常频繁。尽管大多数博客都提供了这个问题的答案，但细节仍然缺失。我研究了基本假设，并希望与您分享我的发现。首先，我会简单地告诉你这些假设，然后举例说明。
线性回归模型的基本假设如下:
- 自变量(X)和因变量(y)之间存在线性关系
- 不同特征之间很少或没有多重共线性
- 残差应呈正态分布(多变量正态性)
- 残基之间很少或没有自相关
- 误差的同方差
现在，让我们来看看如何验证一个假设，以及在假设不成立的情况下应该做些什么。让我们逐一关注这些要点。我们将获取一个包含不同葡萄酒特征的数据集。该数据集已被数据科学家同事用于多个示例，并由 UCI 机器学习知识库公开提供( Wine_quality data 或来自此处的 CSV 文件)。我将使用的另一个数据集是温度数据集(可从此处获得)。我在纳格什·辛格·肖汉的一篇文章中偶然发现了这些数据集。我会建议你下载数据，用它来寻找行数，列数，是否有 NaN 值的行等。熊猫是一个阅读 CSV 和处理数据的非常好的工具。如果您不熟悉熊猫，请尝试使用以下方式阅读该文件:
现在可以用，dataset . head()/dataset . tail()/dataset . describe 等来玩数据了。
- 线性关系
自变量和因变量之间应该存在线性关系。这很容易用散点图来验证。我将使用温度数据集来显示线性关系。您可以绘制 Tmax vrs T_min(或 T_avg vrs T_min ),如图 1 所示。如果你已经处理过数据，你可能已经观察到有一个月列，因此我们甚至可以根据月份标记(颜色代码)散点图，只是为了看看不同月份的温度是否有明显的区别(图 1b)。
输出应该类似于:
图 1
这里可以观察到 T_max 和 T_min 遵循线性趋势。你可以想象一条直线穿过数据。因此，我们已经确保我们的数据遵循第一个假设。
如果数据不是线性的呢？
- 也许用线性模型拟合数据是一个错误的想法。它可能更适合多项式模型(非线性回归)。
- 数据转换，如果 y 似乎是 x 的指数，那么在 y 和 log(x)之间画一条曲线怎么样(或者 y 对 x 的平方)。现在你可以对这些数据进行线性回归。但是你为什么要这么做呢？因为你(我也是)比非线性回归更懂线性回归。我们已经为线性回归、假设验证等建立了许多工具，这些工具对于非线性回归可能并不容易获得。此外，一旦你拟合了 y 和变换后的 x 之间的线性回归，就不难回到原来的 y 对 x 的关系。
2。不同特征之间没有多重共线性
为什么‘无多重共线性’？
当我们进行线性回归分析时，我们在寻找 y = mx + c 类型的解，其中 c 是截距，m 是斜率。“m”的值决定了将 x 改变 1 时 y 将改变多少。对于多元线性回归，同样的关系适用于以下等式:y = m1x1 +m2x2 +m3x3 … + c。理想情况下，m1 表示 y 在改变 x1 时会改变多少，但如果 x1 的改变会改变 x2 或 x3 呢？在这种情况下，y 和 m1(或 m2、m3 等)之间的关系将非常复杂。
如何检查‘多重共线性’？
Seaborn 提供了一个 pairplot 函数，它可以绘制变量之间的属性。这些图是散点图，我们需要看看这些属性是否呈现线性关系。这是可视化和感受不同属性的线性关系的最简单的工具，但只有当涉及的要素数量限制在 10-12 个时才是好的。在下一节中，我们将讨论如果涉及到更多的特性该怎么办。让我们先画出我们的配对图。为此，我将使用 Wine_quality 数据，因为它具有高度相关的特性(图 2)。
输出应该是 11x11 的图形，如下所示:
图 2
如果你观察像 pH 值和固定酸度显示线性相关性(具有负协方差)。如果你观察完整的情节，你会发现
- 固定酸度与柠檬酸、密度和酸度相关
- 挥发性酸度与柠檬酸相关
我让你去寻找其他的共线性关系。让我们绕道了解一下这种共线性的原因。如果你还记得你的高中化学，pH 值被定义为
pH =-log[H+]=-log(酸的浓度)
由此很直观的得出 pH 值和柠檬酸或挥发酸度是负相关的。
本练习还提供了一个关于数据的领域知识如何帮助更有效地处理数据的示例。大概这就是数据科学对所有科学领域的科学家开放的原因。
不同特征的共线程度如何？
在上一节中，我们绘制了不同的要素，以检查它们是否共线。在本节中，我们将回答什么是共线性的度量？
我们可以测量相关性(注意“相关性”不是“共线性”)，如果两个特征之间的绝对相关性很高，我们可以说这两个特征是共线的。为了测量不同特征之间的相关性，我们使用相关矩阵/热图。为了绘制热图，我们可以使用 seaborn 的热图函数(图 3)。
输出应为彩色编码矩阵，并在网格中标注相关性:
图 3
现在，根据您的统计知识，您可以决定一个阈值，如 0.4 或 0.5，如果相关性大于此阈值，则认为是一个问题。与其在这里给出一个明确的答案，我不如向你提一个问题。这个相关阈值的理想值应该是多少？现在，假设您的数据集包含 10，000 个示例(或行),如果数据集包含 100，000 个或 1000 个示例，您会改变答案吗？也许你可以在评论中给出你的答案。
相关性高怎么办？
假设您已经列出了不同特征之间的共线关系。现在怎么办？首先要考虑的是一个特性是否可以被删除。如果两个特征直接相关，例如酸度和 pH 值，我会毫不犹豫地删除其中一个。因为 pH 只不过是酸量的负对数。这就像在两个不同的尺度上拥有相同的信息。但更大的问题是:
- 哪个功能可以去除 pH 值或酸量？
- 接下来应该删除哪个功能？
这个我没有明确的答案。我所学的是计算方差通货膨胀系数 VIF。它被定义为公差的倒数，而公差是 1- R2。虽然我将讨论 VIF，但通常有以下方法可用于处理共线性:
a)贪婪淘汰
b)递归特征消除
c)套索正则化(L1 正则化)
d)主成分分析
我们将使用 VIF 值来查找应该首先消除的要素。然后，我们将重新计算 VIF，以检查是否有任何其他功能需要消除。如果我们在相关标度上工作，不同变量之间的相关性在消除前后不会改变。VIF 给出了这个优势来衡量淘汰的效果。我们使用 statsmodels，oulier_influence 模块来计算 VIF。此模型要求我们在模型中添加一个常量变量来计算 VIF，因此在代码中我们使用“add_constant(X)”,其中 X 是包含所有要素的数据集(质量列被删除，因为它包含目标值)。
我采用的方法是消除具有最高 VIF 的要素，然后重新计算 VIF。如果 VIF 值大于 10，则移除 VIF 次高的要素，否则我们将不再处理多重共线性。
为了检验其他假设，我们需要进行线性回归。我已经使用 scikit 学习线性回归模块来做同样的事情。我们将模型分为测试和训练模型，使用训练数据拟合模型，并使用测试数据进行预测。
3。 多元常态
残差应该是正态分布的。这一点可以通过绘制 QQ 图很容易地检查出来。我们将使用 statsmodels，qqplot 来绘制它。我们首先导入 qqplot 属性，然后向它提供残差值。我们得到的 Q-Q 图如图 4 所示。
图 4
理想情况下，它应该是一条直线。这种模式表明我们的模型有严重的问题。我们不能依赖这个回归模型。让我们检查一下其他假设是否成立。
4。无自相关
Durbin Watson 的 d 检验可以帮助我们分析残基之间是否存在任何自相关。既然网上有很多关于这个测试的资料，我就给你提供另一种方式。在图 5 中显示了每个属性的残基图，以检查残基是否显示任何相关性。在下面的代码中，dataset2 是 X_test 的 pandas 数据帧。
输出将是一系列图(测试数据集的 1 个图/列)
图 5
图 5 显示了数据是如何在没有任何特定模式的情况下很好地分布的，从而验证了残留物没有自相关。我们需要在所有的图中验证这一点(X 轴是特征，所以有多少个特征就有多少个图)。
5。同质性
用拟合线绘制误差散点图将显示残留物是否与该线形成任何模式。如果是，那么数据不是异方差的或者数据是异方差的。而如果散点图不形成任何模式，并且随机分布在拟合线周围，则残差是均方的。
为了绘制关于拟合线的残差(y _ test-y _ pred ),可以写出拟合线的方程(通过使用*。coeff_ 和*。拦截)。使用此等式获得 y 值，但将这些 y 值绘制在 X 轴上，因为我们想要绘制关于拟合线的残差(X 轴应该是拟合线)。现在可以观察到残留物的模式。下面是相同的代码:
如果您是 python 新手，并且希望远离编写代码，您可以使用 yellowbrick regressor 的“Redidualsplot”模块来执行相同的任务。在这里，您只需要适应测试和训练数据，其余的将由模型本身完成。在这里，我使用 scikit learn 的 LinearRegression()模型，您可以选择使用不同的模型。
无论如何，使用上述任何一种方法都会得到相同的结果，如图 6 所示。
这里的残留物显示了一个清晰的模式，表明我们的模型有问题。
因此，我们的模型无法支持多元正态性和同方差假设(分别见图 4 和图 6)。这表明，要么数据不适合线性回归，要么给定的特征不能基于给定的特征真正预测葡萄酒的质量。但这是展示线性回归基本假设的一个很好的练习。如果你对这个问题有更好的解决办法，请告诉我。
欢迎建设性的批评/建议。
本主题的其他好读物:
手工线性回归
线性回归是数据科学家最基本也是最强大的工具。让我们仔细看看最小二乘直线和相关系数。
线性回归的发明
Johannes Plenio 在 Unsplash 上拍摄的照片
线性回归是线性代数的一种形式，据称是由卡尔·弗里德里希·高斯(1777-1855)发明的，但最早发表在阿德里安·玛丽·勒让德(1752-1833)的一篇科学论文中。高斯用最小二乘法猜测谷神星小行星何时何地会出现在夜空中(统计回归的发现，2015)。这不是一个爱好项目，这是一个资金充足的研究项目，目的是为了海洋导航，这是一个高度竞争的领域，对技术中断非常敏感。
线性回归原理
线性回归是一种从 x 预测 y 的方法 在我们的例子中， y 是因变量，x 是自变量。 我们想要预测给定的x 值的 y 值。现在，如果数据是完全线性的，我们可以简单地根据 y = mx+ b 计算直线的斜率截距形式。要预测 y ，我们只需插入给定的 x 和 b 的值。 在现实世界中，我们的数据不会是完全线性的。它很可能以散点图上的数据点簇的形式出现。从散点图中，我们将确定， 描述数据线性质量的最佳拟合线 是什么，以及 该线与点群的拟合程度如何？
线性回归试图通过将线性方程拟合到观察到的数据来模拟两个变量之间的关系( 线性回归 ，n.d .)。
散点图
让我们编造一些数据作为例子。黑猩猩狩猎团体的规模和成功狩猎的百分比之间的关系已经被很好的记录了。(Busse，1978)我将从 Busse 获取一些数据点用于本文，并使用 seaborn 散点图绘制数据。注意到我在数据中画的线并不完全符合它，但是这些点近似于一个线性模式吗？我通过数据画出的线是最小二乘法线**，用来预测给定 x 值 的 y 值。仅使用手工绘制的基本最小二乘法线，我们可以预测 4 只黑猩猩的狩猎队将有大约 52%的成功率。我们不是 100%准确，但随着更多的数据，我们可能会提高我们的准确性。数据与最小二乘法直线的拟合程度就是相关系数*。*****
最小平方线
在上面的图表中，我只是通过我判断为最佳拟合的数据手工画了一条线。我们要用斜率截距形式 y = mx + b 来计算这条线，才能做出真正的预测。我们寻求的是一条直线，在这条直线上，直线和每个点之间的差异尽可能小。这是最佳拟合线。
最小二乘直线被定义为从数据点到直线的垂直距离的平方和尽可能小的直线(Lial，格伦威尔和 Ritchey，2016)。
最小二乘法直线有两个分量:斜率 m、 和 y 截距 b. 我们将首先求解 m ，然后求解*b .*m和 b 的方程为:
在 MS Word 公式编辑器中创建
那可是好多 Sigmas)！。不过不用担心，适马只是表示“的总和”，比如“x 的总和”，用∑x 来象征，也就是 x 列的总和，“黑猩猩的数量”我们需要计算∑x、∑y、∑xy、∑x 和∑y。然后，将每一个部分输入到等式中，用于计算和b。根据我们的原始数据集创建下表。
现在很简单，将我们的适马值插入到 m 和 b 的公式中。n 是数据集中值的数量，在我们的例子中是 8。
你有它！你可以根据给定的 x 的值来预测 y ，使用你的等式: y = 5.4405x + 31.6429。 这意味着我们的线从 31.6429 开始，每有一只黑猩猩加入狩猎队，y 值就会增加 5.4405 个百分点。为了验证这一点，让我们来预测 4 只黑猩猩的狩猎成功率。
****y = 5.4405(4)+31.6429***, which results in **y=53.4***
我们刚刚预测了黑猩猩狩猎队狩猎成功的百分比，这仅仅是基于对他们群体大小的了解，这是相当惊人的！
让我们使用 python 在之前的散点图上绘制最小二乘直线，以展示它如何拟合数据。Seaborn.regplot()是在这种情况下使用的一个很好的图表，但是出于演示的目的，我将手动创建 y=mx+b 线 并将其放置在 seaborn 图表上。
然而，现在您可以进行预测了，您需要用相关系数来限定您的预测，相关系数描述了数据与您的计算线的吻合程度。**
相关系数
我们使用相关系数来确定最小二乘法是否是我们数据的好模型。如果数据点不是线性的，那么直线就不是正确的预测模型。卡尔·皮尔逊发明了相关系数 r ，介于 1 和-1 之间，衡量两个变量之间线性关系的强度(Lial，格伦威尔和 Ritchey，2016)。如果 r 正好是-1 或 1，则表示数据正好符合、线，没有偏离线。 r=0 表示没有线性相关。当 r 值 接近零时，意味着关联度也降低。
相关系数由以下公式描述
幸运的是，这些适马值已经在前面的表格中计算过了。我们简单地把它们代入我们的方程。
我们的值接近正 1，这意味着数据是高度相关的，并且是正的。你可以通过观察散点图上的最小二乘方线来确定这一点，但是相关系数给了你科学的证据！
结论
线性回归是数据科学家或统计学家可用的最佳机器学习方法之一。有许多方法可以使用您的编程技能来创建机器学习模型，但让您自己熟悉模型使用的数学绝对是一个好主意。
参考
***布塞博士(1978 年)。黑猩猩会合作捕猎吗？美国博物学家， 112 (986)，767–770。【https://doi.org/10.1086/283318 ***
利亚尔，格伦威尔和里奇(2016)。有限数学与微积分及应用，第 10 版。纽约州纽约市:皮尔森[ISBN-13 9780133981070]。
线性回归。(未注明)。检索于 2020 年 4 月 11 日，来自http://www.stat.yale.edu/Courses/1997-98/101/linreg.htm
统计回归的发现。(2015 年 11 月 6 日)。价格经济学。http://priceonomics . com/the-discovery-of-statistical-regression/
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
