A. python 求解线性方程
问题有问题吧,不过我把应该做的过程给你。迹卜
# 结果
r = []
# 步长
step = 0.00000000001
# 初始化x
x = step
while x < 10:
# 计算y
y = 0.5 * (10 - x)
# 判断符合条件的x, y
if y - x - 0.1 >= 0 and x - y - 0.5 >= 0:
r.append((x, y))
x += step
结果
设置步长,步长姿链穗越小越精准。
先把方程转换一下 10*x+20*y=100 -> y = 0.5 * (10 - x)
题中,y >唤念 0, 所以 10 - x > 0, 所以 x < 10, 遍历 0 < x < 10。记录所有符合条件的x,y
B. 万字教你如何用 Python 实现线性规划
想象一下,您有一个线性方程组和不等式系统。这样的系统通常有许多可能的解决方案。线性规划是一组数学和计算工具,可让您找到该系统的特定解,该解对应于某些其他线性函数的最大值或最小值。
混合整数线性规划是 线性规划 的扩展。它处理至少一个变量采用离散整数而不是连续值的问题。尽管乍一看混合整数问题与连续变量问题相似,但它们在灵活性和精度方面具有显着优势。
整数变量对于正确表示自然用整数表示的数量很重要,例如生产的飞机数量或服务的客户数量。
一种特别重要的整数变量是 二进制变量 。它只能取 零 或 一 的值,在做出是或否的决定时很有用,例如是否应该建造工厂或者是否应该打开或关闭机器。您还可以使用它们来模拟逻辑约束。
线性规划是一种基本的优化技术,已在科学和数学密集型领域使用了数十年。它精确、相对快速,适用于一系列实际应用。
混合整数线性规划允许您克服线性规划的许多限制。您可以使用分段线性函数近似非线性函数、使用半连续变量、模型逻辑约束等。它是一种计算密集型工具,但计算机硬件和软件的进步使其每天都更加适用。
通常,当人们试图制定和解决优化问题时,第一个问题是他们是否可以应用线性规划或混合整数线性规划。
以下文章说明了线性规划和混合整数线性规划的一些用例:
随着计算机能力的增强、算法的改进以及更多用户友好的软件解决方案的出现,线性规划,尤其是混合整数线性规划的重要性随着时间的推移而增加。
解决线性规划问题的基本方法称为,它有多种变体。另一种流行的方法是。
混合整数线性规划问题可以通过更复杂且计算量更大的方法来解决,例如,它在幕后使用线性规划。这种方法的一些变体是,它涉及使用 切割平面 ,以及。
有几种适用于线性规划和混合整数线性规划的合适且众所周知的 Python 工具。其中一些是开源的,而另一些是专有的。您是否需要免费或付费工具取决于问题的规模和复杂性,以及对速度和灵活性的需求。
值得一提的是,几乎所有广泛使用的线性规划和混合整数线性规划库都是以 Fortran 或 C 或 C++ 原生和编写的。这是因为线性规划需要对(通常很大)矩阵进行计算密集型工作。此类库称为求解器。Python 工具只是求解器的包装器。
Python 适合围绕本机库构建包装器,因为它可以很好地与 C/C++ 配合使用。对于本教程,您不需要任何 C/C++(或 Fortran),但如果您想了解有关此酷功能的更多信息,请查看以下资源:
基本上,当您定义和求解模型时,您使用 Python 函数或方法调用低级库,该库执行实际优化工作并将解决方案返回给您的 Python 对象。
几个免费的 Python 库专门用于与线性或混合整数线性规划求解器交互:
在本教程中,您将使用SciPy和PuLP来定义和解决线性规划问题。
在本节中,您将看到线性规划问题的两个示例:
您将在下一节中使用 Python 来解决这两个问题。
考虑以下线性规划问题:
你需要找到X和Ÿ使得红色,蓝色和黄色的不平等,以及不平等X 0和ÿ 0,是满意的。同时,您的解决方案必须对应于z的最大可能值。
您需要找到的自变量(在本例中为 x 和 y )称为 决策变量 。要最大化或最小化的决策变量的函数(在本例中为 z) 称为 目标函数 、 成本函数 或仅称为 目标 。您需要满足的 不等式 称为 不等式约束 。您还可以在称为 等式约束 的约束中使用方程。
这是您如何可视化问题的方法:
红线代表的功能2 X + Ý = 20,和它上面的红色区域示出了红色不等式不满足。同样,蓝线是函数 4 x + 5 y = 10,蓝色区域被禁止,因为它违反了蓝色不等式。黄线是 x + 2 y = 2,其下方的黄色区域是黄色不等式无效的地方。
如果您忽略红色、蓝色和黄色区域,则仅保留灰色区域。灰色区域的每个点都满足所有约束,是问题的潜在解决方案。该区域称为 可行域 ,其点为 可行解 。在这种情况下,有无数可行的解决方案。
您想最大化z。对应于最大z的可行解是 最优解 。如果您尝试最小化目标函数,那么最佳解决方案将对应于其可行的最小值。
请注意,z是线性的。你可以把它想象成一个三维空间中的平面。这就是为什么最优解必须在可行区域的 顶点 或角上的原因。在这种情况下,最佳解决方案是红线和蓝线相交的点,稍后您将看到。
有时,可行区域的整个边缘,甚至整个区域,都可以对应相同的z值。在这种情况下,您有许多最佳解决方案。
您现在已准备好使用绿色显示的附加等式约束来扩展问题:
方程式 x + 5 y = 15,以绿色书写,是新的。这是一个等式约束。您可以通过向上一张图像添加相应的绿线来将其可视化:
现在的解决方案必须满足绿色等式,因此可行区域不再是整个灰色区域。它是绿线从与蓝线的交点到与红线的交点穿过灰色区域的部分。后一点是解决方案。
如果插入x的所有值都必须是整数的要求,那么就会得到一个混合整数线性规划问题,可行解的集合又会发生变化:
您不再有绿线,只有沿线的x值为整数的点。可行解是灰色背景上的绿点,此时最优解离红线最近。
这三个例子说明了 可行的线性规划问题 ,因为它们具有有界可行区域和有限解。
如果没有解,线性规划问题是 不可行的 。当没有解决方案可以同时满足所有约束时,通常会发生这种情况。
例如,考虑如果添加约束x + y 1会发生什么。那么至少有一个决策变量(x或y)必须是负数。这与给定的约束x 0 和y 0相冲突。这样的系统没有可行的解决方案,因此称为不可行的。
另一个示例是添加与绿线平行的第二个等式约束。这两行没有共同点,因此不会有满足这两个约束的解决方案。
一个线性规划问题是 无界的 ,如果它的可行区域是无界,将溶液不是有限。这意味着您的变量中至少有一个不受约束,可以达到正无穷大或负无穷大,从而使目标也无限大。
例如,假设您采用上面的初始问题并删除红色和黄色约束。从问题中删除约束称为 放松 问题。在这种情况下,x和y不会在正侧有界。您可以将它们增加到正无穷大,从而产生无限大的z值。
在前面的部分中,您研究了一个与任何实际应用程序无关的抽象线性规划问题。在本小节中,您将找到与制造业资源分配相关的更具体和实用的优化问题。
假设一家工厂生产四种不同的产品,第一种产品的日产量为x ₁,第二种产品的产量为x 2,依此类推。目标是确定每种产品的利润最大化日产量,同时牢记以下条件:
数学模型可以这样定义:
目标函数(利润)在条件 1 中定义。人力约束遵循条件 2。对原材料 A 和 B 的约束可以从条件 3 和条件 4 中通过对每种产品的原材料需求求和得出。
最后,产品数量不能为负,因此所有决策变量必须大于或等于零。
与前面的示例不同,您无法方便地将其可视化,因为它有四个决策变量。但是,无论问题的维度如何,原理都是相同的。
在本教程中,您将使用两个Python 包来解决上述线性规划问题:
SciPy 设置起来很简单。安装后,您将拥有开始所需的一切。它的子包 scipy.optimize 可用于线性和非线性优化。
PuLP 允许您选择求解器并以更自然的方式表述问题。PuLP 使用的默认求解器是COIN-OR Branch and Cut Solver (CBC)。它连接到用于线性松弛的COIN-OR 线性规划求解器 (CLP)和用于切割生成的COIN-OR 切割生成器库 (CGL)。
另一个伟大的开源求解器是GNU 线性规划工具包 (GLPK)。一些着名且非常强大的商业和专有解决方案是Gurobi、CPLEX和XPRESS。
除了在定义问题时提供灵活性和运行各种求解器的能力外,PuLP 使用起来不如 Pyomo 或 CVXOPT 等替代方案复杂,后者需要更多的时间和精力来掌握。
要学习本教程,您需要安装 SciPy 和 PuLP。下面的示例使用 SciPy 1.4.1 版和 PuLP 2.1 版。
您可以使用pip以下方法安装两者:
您可能需要运行pulptest或sudo pulptest启用 PuLP 的默认求解器,尤其是在您使用 Linux 或 Mac 时:
或者,您可以下载、安装和使用 GLPK。它是免费和开源的,适用于 Windows、MacOS 和 Linux。在本教程的后面部分,您将看到如何将 GLPK(除了 CBC)与 PuLP 一起使用。
在 Windows 上,您可以下载档案并运行安装文件。
在 MacOS 上,您可以使用 Homebrew:
在 Debian 和 Ubuntu 上,使用apt来安装glpk和glpk-utils:
在Fedora,使用dnf具有glpk-utils:
您可能还会发现conda对安装 GLPK 很有用:
安装完成后,可以查看GLPK的版本:
有关详细信息,请参阅 GLPK 关于使用Windows 可执行文件和Linux 软件包进行安装的教程。
在本节中,您将学习如何使用 SciPy优化和求根库进行线性规划。
要使用 SciPy 定义和解决优化问题,您需要导入scipy.optimize.linprog():
现在您已经linprog()导入,您可以开始优化。
让我们首先解决上面的线性规划问题:
linprog()仅解决最小化(而非最大化)问题,并且不允许具有大于或等于符号 ( ) 的不等式约束。要解决这些问题,您需要在开始优化之前修改您的问题:
引入这些更改后,您将获得一个新系统:
该系统与原始系统等效,并且将具有相同的解决方案。应用这些更改的唯一原因是克服 SciPy 与问题表述相关的局限性。
下一步是定义输入值:
您将上述系统中的值放入适当的列表、元组或NumPy 数组中:
注意:请注意行和列的顺序!
约束左侧和右侧的行顺序必须相同。每一行代表一个约束。
来自目标函数和约束左侧的系数的顺序必须匹配。每列对应一个决策变量。
下一步是以与系数相同的顺序定义每个变量的界限。在这种情况下,它们都在零和正无穷大之间:
此语句是多余的,因为linprog()默认情况下采用这些边界(零到正无穷大)。
注:相反的float("inf"),你可以使用math.inf,numpy.inf或scipy.inf。
最后,是时候优化和解决您感兴趣的问题了。你可以这样做linprog():
参数c是指来自目标函数的系数。A_ub和b_ub分别与不等式约束左边和右边的系数有关。同样,A_eq并b_eq参考等式约束。您可以使用bounds提供决策变量的下限和上限。
您可以使用该参数method来定义要使用的线性规划方法。有以下三种选择:
linprog() 返回具有以下属性的数据结构:
您可以分别访问这些值:
这就是您获得优化结果的方式。您还可以以图形方式显示它们:
如前所述,线性规划问题的最优解位于可行区域的顶点。在这种情况下,可行区域只是蓝线和红线之间的绿线部分。最优解是代表绿线和红线交点的绿色方块。
如果要排除相等(绿色)约束,只需删除参数A_eq并b_eq从linprog()调用中删除:
解决方案与前一种情况不同。你可以在图表上看到:
在这个例子中,最优解是红色和蓝色约束相交的可行(灰色)区域的紫色顶点。其他顶点,如黄色顶点,具有更高的目标函数值。
您可以使用 SciPy 来解决前面部分所述的资源分配问题:
和前面的例子一样,你需要从上面的问题中提取必要的向量和矩阵,将它们作为参数传递给.linprog(),然后得到结果:
结果告诉您最大利润是1900并且对应于x ₁ = 5 和x ₃ = 45。在给定条件下生产第二和第四个产品是没有利润的。您可以在这里得出几个有趣的结论:
opt.statusis0和opt.successis True,说明优化问题成功求解,最优可行解。
SciPy 的线性规划功能主要用于较小的问题。对于更大和更复杂的问题,您可能会发现其他库更适合,原因如下:
幸运的是,Python 生态系统为线性编程提供了几种替代解决方案,这些解决方案对于更大的问题非常有用。其中之一是 PuLP,您将在下一节中看到它的实际应用。
PuLP 具有比 SciPy 更方便的线性编程 API。您不必在数学上修改您的问题或使用向量和矩阵。一切都更干净,更不容易出错。
像往常一样,您首先导入您需要的内容:
现在您已经导入了 PuLP,您可以解决您的问题。
您现在将使用 PuLP 解决此系统:
第一步是初始化一个实例LpProblem来表示你的模型:
您可以使用该sense参数来选择是执行最小化(LpMinimize或1,这是默认值)还是最大化(LpMaximize或-1)。这个选择会影响你的问题的结果。
一旦有了模型,就可以将决策变量定义为LpVariable类的实例:
您需要提供下限,lowBound=0因为默认值为负无穷大。该参数upBound定义了上限,但您可以在此处省略它,因为它默认为正无穷大。
可选参数cat定义决策变量的类别。如果您使用的是连续变量,则可以使用默认值"Continuous"。
您可以使用变量x和y创建表示线性表达式和约束的其他 PuLP 对象:
当您将决策变量与标量相乘或构建多个决策变量的线性组合时,您会得到一个pulp.LpAffineExpression代表线性表达式的实例。
注意:您可以增加或减少变量或表达式,你可以乘他们常数,因为纸浆类实现一些Python的特殊方法,即模拟数字类型一样__add__(),__sub__()和__mul__()。这些方法用于像定制运营商的行为+,-和*。
类似地,您可以将线性表达式、变量和标量与运算符 ==、=以获取表示模型线性约束的纸浆.LpConstraint实例。
注:也有可能与丰富的比较方法来构建的约束.__eq__(),.__le__()以及.__ge__()定义了运营商的行为==,=。
考虑到这一点,下一步是创建约束和目标函数并将它们分配给您的模型。您不需要创建列表或矩阵。只需编写 Python 表达式并使用+=运算符将它们附加到模型中:
在上面的代码中,您定义了包含约束及其名称的元组。LpProblem允许您通过将约束指定为元组来向模型添加约束。第一个元素是一个LpConstraint实例。第二个元素是该约束的可读名称。
设置目标函数非常相似:
或者,您可以使用更短的符号:
现在您已经添加了目标函数并定义了模型。
注意:您可以使用运算符将 约束或目标附加到模型中,+=因为它的类LpProblem实现了特殊方法.__iadd__(),该方法用于指定 的行为+=。
对于较大的问题,lpSum()与列表或其他序列一起使用通常比重复+运算符更方便。例如,您可以使用以下语句将目标函数添加到模型中:
它产生与前一条语句相同的结果。
您现在可以看到此模型的完整定义:
模型的字符串表示包含所有相关数据:变量、约束、目标及其名称。
注意:字符串表示是通过定义特殊方法构建的.__repr__()。有关 的更多详细信息.__repr__(),请查看Pythonic OOP 字符串转换:__repr__vs__str__ .
最后,您已准备好解决问题。你可以通过调用.solve()你的模型对象来做到这一点。如果要使用默认求解器 (CBC),则不需要传递任何参数:
.solve()调用底层求解器,修改model对象,并返回解决方案的整数状态,1如果找到了最优解。有关其余状态代码,请参阅LpStatus[]。
你可以得到优化结果作为 的属性model。该函数value()和相应的方法.value()返回属性的实际值:
model.objective持有目标函数model.constraints的值,包含松弛变量的值,以及对象x和y具有决策变量的最优值。model.variables()返回一个包含决策变量的列表:
如您所见,此列表包含使用 的构造函数创建的确切对象LpVariable。
结果与您使用 SciPy 获得的结果大致相同。
注意:注意这个方法.solve()——它会改变对象的状态,x并且y!
您可以通过调用查看使用了哪个求解器.solver:
输出通知您求解器是 CBC。您没有指定求解器,因此 PuLP 调用了默认求解器。
如果要运行不同的求解器,则可以将其指定为 的参数.solve()。例如,如果您想使用 GLPK 并且已经安装了它,那么您可以solver=GLPK(msg=False)在最后一行使用。请记住,您还需要导入它:
现在你已经导入了 GLPK,你可以在里面使用它.solve():
该msg参数用于显示来自求解器的信息。msg=False禁用显示此信息。如果要包含信息,则只需省略msg或设置msg=True。
您的模型已定义并求解,因此您可以按照与前一种情况相同的方式检查结果:
使用 GLPK 得到的结果与使用 SciPy 和 CBC 得到的结果几乎相同。
一起来看看这次用的是哪个求解器:
正如您在上面用突出显示的语句定义的那样model.solve(solver=GLPK(msg=False)),求解器是 GLPK。
您还可以使用 PuLP 来解决混合整数线性规划问题。要定义整数或二进制变量,只需传递cat="Integer"或cat="Binary"到LpVariable。其他一切都保持不变:
在本例中,您有一个整数变量并获得与之前不同的结果:
Nowx是一个整数,如模型中所指定。(从技术上讲,它保存一个小数点后为零的浮点值。)这一事实改变了整个解决方案。让我们在图表上展示这一点:
如您所见,最佳解决方案是灰色背景上最右边的绿点。这是两者的最大价值的可行的解决方案x和y,给它的最大目标函数值。
GLPK 也能够解决此类问题。
现在你可以使用 PuLP 来解决上面的资源分配问题:
定义和解决问题的方法与前面的示例相同:
在这种情况下,您使用字典 x来存储所有决策变量。这种方法很方便,因为字典可以将决策变量的名称或索引存储为键,将相应的LpVariable对象存储为值。列表或元组的LpVariable实例可以是有用的。
上面的代码产生以下结果:
如您所见,该解决方案与使用 SciPy 获得的解决方案一致。最有利可图的解决方案是每天生产5.0第一件产品和45.0第三件产品。
让我们把这个问题变得更复杂和有趣。假设由于机器问题,工厂无法同时生产第一种和第三种产品。在这种情况下,最有利可图的解决方案是什么?
现在您有另一个逻辑约束:如果x ₁ 为正数,则x ₃ 必须为零,反之亦然。这是二元决策变量非常有用的地方。您将使用两个二元决策变量y ₁ 和y ₃,它们将表示是否生成了第一个或第三个产品:
除了突出显示的行之外,代码与前面的示例非常相似。以下是差异:
这是解决方案:
事实证明,最佳方法是排除第一种产品而只生产第三种产品。
就像有许多资源可以帮助您学习线性规划和混合整数线性规划一样,还有许多具有 Python 包装器的求解器可用。这是部分列表:
其中一些库,如 Gurobi,包括他们自己的 Python 包装器。其他人使用外部包装器。例如,您看到可以使用 PuLP 访问 CBC 和 GLPK。
您现在知道什么是线性规划以及如何使用 Python 解决线性规划问题。您还了解到 Python 线性编程库只是本机求解器的包装器。当求解器完成其工作时,包装器返回解决方案状态、决策变量值、松弛变量、目标函数等。
C. python多元线性回归怎么计算
1、什么是多元线性回归模型?
当y值的影响因素不唯一时,采用多元线性回归模型。
y =y=β0+β1x1+β2x2+...+βnxn
例如商品的销售额可能不电视广告投入,收音机广告投入,报纸广告投入有关系,可以有 sales =β0+β1*TV+β2* radio+β3*newspaper.
2、使用pandas来读取数据
pandas 是一个用于数据探索、数据分析和数据处理的python库
[python]view plain
importpandasaspd
[html]view plain
<prename="code"class="python">#
data=pd.read_csv('/home/lulei/Advertising.csv')
#displaythefirst5rows
data.head()
上面代码的运行结果:
上面显示的结果类似一个电子表格,这个结构称为Pandas的数据帧(data frame),类型全称:pandas.core.frame.DataFrame.
pandas的两个主要数据结构:Series和DataFrame:
Series类似于一维数组,它有一组数据以及一组与之相关的数据标签(即索引)组成。
DataFrame是一个表格型的数据结构,它含有一组有序的列,每列可以是不同的值类型。DataFrame既有行索引也有列索引,它可以被看做由Series组成的字典。
[python]view plain
#displaythelast5rows
data.tail()
只显示结果的末尾5行
[html]view plain
#checktheshapeoftheDataFrame(rows,colums)
data.shape
查看DataFrame的形状,注意第一列的叫索引,和数据库某个表中的第一列类似。
(200,4)
3、分析数据
特征:
TV:对于一个给定市场中单一产品,用于电视上的广告费用(以千为单位)
Radio:在广播媒体上投资的广告费用
Newspaper:用于报纸媒体的广告费用
响应:
Sales:对应产品的销量
在这个案例中,我们通过不同的广告投入,预测产品销量。因为响应变量是一个连续的值,所以这个问题是一个回归问题。数据集一共有200个观测值,每一组观测对应一个市场的情况。
注意:这里推荐使用的是seaborn包。网上说这个包的数据可视化效果比较好看。其实seaborn也应该属于matplotlib的内部包。只是需要再次的单独安装。
[python]view plain
importseabornassns
importmatplotlib.pyplotasplt
#ots
sns.pairplot(data,x_vars=['TV','Radio','Newspaper'],y_vars='Sales',size=7,aspect=0.8)
plt.show()#注意必须加上这一句,否则无法显示。
[html]view plain
这里选择TV、Radio、Newspaper作为特征,Sales作为观测值
[html]view plain
返回的结果:
[python]view plain
sns.pairplot(data,x_vars=['TV','Radio','Newspaper'],y_vars='Sales',size=7,aspect=0.8,kind='reg')
plt.show()
结果显示如下:
4、线性回归模型
优点:快速;没有调节参数;可轻易解释;可理解。
缺点:相比其他复杂一些的模型,其预测准确率不是太高,因为它假设特征和响应之间存在确定的线性关系,这种假设对于非线性的关系,线性回归模型显然不能很好的对这种数据建模。
线性模型表达式:y=β0+β1x1+β2x2+...+βnxn其中
y是响应
β0是截距
β1是x1的系数,以此类推
在这个案例中:y=β0+β1∗TV+β2∗Radio+...+βn∗Newspaper
(1)、使用pandas来构建X(特征向量)和y(标签列)
scikit-learn要求X是一个特征矩阵,y是一个NumPy向量。
pandas构建在NumPy之上。
因此,X可以是pandas的DataFrame,y可以是pandas的Series,scikit-learn可以理解这种结构。
[python]view plain
#
feature_cols=['TV','Radio','Newspaper']
#
X=data[feature_cols]
#
X=data[['TV','Radio','Newspaper']]
#printthefirst5rows
printX.head()
#checkthetypeandshapeofX
printtype(X)
printX.shape
输出结果如下:
[python]view plain
#selectaSeriesfromtheDataFrame
y=data['Sales']
#
y=data.Sales
#printthefirst5values
printy.head()
输出的结果如下:
(2)、构建训练集与测试集
[html]view plain
<prename="code"class="python"><spanstyle="font-size:14px;">##构造训练集和测试集
fromsklearn.cross_validationimporttrain_test_split#这里是引用了交叉验证
X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=1)
#default split is 75% for training and 25% for testing
[html]view plain
printX_train.shape
printy_train.shape
printX_test.shape
printy_test.shape
输出结果如下:
注:上面的结果是由train_test_spilit()得到的,但是我不知道为什么我的版本的sklearn包中居然报错:
处理方法:1、我后来重新安装sklearn包。再一次调用时就没有错误了。
2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)
(3)sklearn的线性回归
[html]view plain
fromsklearn.linear_modelimportLinearRegression
linreg=LinearRegression()
model=linreg.fit(X_train,y_train)
printmodel
printlinreg.intercept_
printlinreg.coef_
输出的结果如下:
[html]view plain
#
zip(feature_cols,linreg.coef_)
输出如下:
y=2.668+0.0464∗TV+0.192∗Radio-0.00349∗Newspaper
如何解释各个特征对应的系数的意义?
对于给定了Radio和Newspaper的广告投入,如果在TV广告上每多投入1个单位,对应销量将增加0.0466个单位。就是加入其它两个媒体投入固定,在TV广告上每增加1000美元(因为单位是1000美元),销量将增加46.6(因为单位是1000)。但是大家注意这里的newspaper的系数居然是负数,所以我们可以考虑不使用newspaper这个特征。这是后话,后面会提到的。
(4)、预测
[python]view plain
y_pred=linreg.predict(X_test)
printy_pred
[python]view plain
printtype(y_pred)
输出结果如下:
5、回归问题的评价测度
(1) 评价测度
对于分类问题,评价测度是准确率,但这种方法不适用于回归问题。我们使用针对连续数值的评价测度(evaluation metrics)。
这里介绍3种常用的针对线性回归的测度。
1)平均绝对误差(Mean Absolute Error, MAE)
(2)均方误差(Mean Squared Error, MSE)
(3)均方根误差(Root Mean Squared Error, RMSE)
这里我使用RMES。
[python]view plain
<prename="code"class="python">#计算Sales预测的RMSE
printtype(y_pred),type(y_test)
printlen(y_pred),len(y_test)
printy_pred.shape,y_test.shape
fromsklearnimportmetrics
importnumpyasnp
sum_mean=0
foriinrange(len(y_pred)):
sum_mean+=(y_pred[i]-y_test.values[i])**2
sum_erro=np.sqrt(sum_mean/50)
#calculateRMSEbyhand
print"RMSEbyhand:",sum_erro
最后的结果如下:
(2)做ROC曲线
[python]view plain
importmatplotlib.pyplotasplt
plt.figure()
plt.plot(range(len(y_pred)),y_pred,'b',label="predict")
plt.plot(range(len(y_pred)),y_test,'r',label="test")
plt.legend(loc="upperright")#显示图中的标签
plt.xlabel("thenumberofsales")
plt.ylabel('valueofsales')
plt.show()
显示结果如下:(红色的线是真实的值曲线,蓝色的是预测值曲线)
直到这里整个的一次多元线性回归的预测就结束了。
6、改进特征的选择
在之前展示的数据中,我们看到Newspaper和销量之间的线性关系竟是负关系(不用惊讶,这是随机特征抽样的结果。换一批抽样的数据就可能为正了),现在我们移除这个特征,看看线性回归预测的结果的RMSE如何?
依然使用我上面的代码,但只需修改下面代码中的一句即可:
[python]view plain
#
feature_cols=['TV','Radio','Newspaper']
#
X=data[feature_cols]
#
#X=data[['TV','Radio','Newspaper']]#只需修改这里即可<prename="code"class="python"style="font-size:15px;line-height:35px;">X=data[['TV','Radio']]#去掉newspaper其他的代码不变
# print the first 5 rowsprint X.head()# check the type and shape of Xprint type(X)print X.shape
最后的到的系数与测度如下:
LinearRegression(_X=True, fit_intercept=True, normalize=False)
然后再次使用ROC曲线来观测曲线的整体情况。我们在将Newspaper这个特征移除之后,得到RMSE变小了,说明Newspaper特征可能不适合作为预测销量的特征,于是,我们得到了新的模型。我们还可以通过不同的特征组合得到新的模型,看看最终的误差是如何的。
备注:
之前我提到了这种错误:
注:上面的结果是由train_test_spilit()得到的,但是我不知道为什么我的版本的sklearn包中居然报错:
处理方法:1、我后来重新安装sklearn包。再一次调用时就没有错误了。
2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)
这里我给出我自己写的函数:
D. 使用Python的线性回归问题,怎么解决
本文中,我们将进行大量的编程——但在这之前,我们先介绍一下我们今天要解决的实例问题。
1) 预测房子价格
闪电侠是一部由剧作家/制片人Greg Berlanti、Andrew Kreisberg和Geoff Johns创作,由CW电视台播放的美国电视连续剧。它基于DC漫画角色闪电侠(Barry Allen),一个具有超人速度移动能力的装扮奇特的打击犯罪的超级英雄,这个角色是由Robert Kanigher、John Broome和Carmine Infantino创作。它是绿箭侠的衍生作品,存在于同一世界。该剧集的试播篇由Berlanti、Kreisberg和Johns写作,David Nutter执导。该剧集于2014年10月7日在北美首映,成为CW电视台收视率最高的电视节目。
绿箭侠是一部由剧作家/制片人 Greg Berlanti、Marc Guggenheim和Andrew Kreisberg创作的电视连续剧。它基于DC漫画角色绿箭侠,一个由Mort Weisinger和George Papp创作的装扮奇特的犯罪打击战士。它于2012年10月10日在北美首映,与2012年末开始全球播出。主要拍摄于Vancouver、British Columbia、Canada,该系列讲述了亿万花花公子Oliver Queen,由Stephen Amell扮演,被困在敌人的岛屿上五年之后,回到家乡打击犯罪和腐败,成为一名武器是弓箭的神秘义务警员。不像漫画书中,Queen最初没有使用化名”绿箭侠“。
由于这两个节目并列为我最喜爱的电视节目头衔,我一直想知道哪个节目更受其他人欢迎——谁会最终赢得这场收视率之战。 所以让我们写一个程序来预测哪个电视节目会有更多观众。 我们需要一个数据集,给出每一集的观众。幸运地,我从维基网络上得到了这个数据,并整理成一个.csv文件。它如下所示。
闪电侠
闪电侠美国观众数
绿箭侠
绿箭侠美国观众数
1 4.83 1 2.84
2 4.27 2 2.32
3 3.59 3 2.55
4 3.53 4 2.49
5 3.46 5 2.73
6 3.73 6 2.6
7 3.47 7 2.64
8 4.34 8 3.92
9 4.66 9 3.06
观众数以百万为单位。
解决问题的步骤:
首先我们需要把数据转换为X_parameters和Y_parameters,不过这里我们有两个X_parameters和Y_parameters。因此,把他们命名为flash_x_parameter、flash_y_parameter、arrow_x_parameter、arrow_y_parameter吧。然后我们需要把数据拟合为两个不同的线性回归模型——先是闪电侠,然后是绿箭侠。 接着我们需要预测两个电视节目下一集的观众数量。 然后我们可以比较结果,推测哪个节目会有更多观众。
步骤1
导入我们的程序包:
Python
1
2
3
4
5
6
7
# Required Packages
import csv
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import datasets, linear_model
步骤2
写一个函数,把我们的数据集作为输入,返回flash_x_parameter、flash_y_parameter、arrow_x_parameter、arrow_y_parameter values。
Python
1
2
3
4
5
6
7
8
9
10
11
12
13
# Function to get data
def get_data(file_name):
data = pd.read_csv(file_name)
flash_x_parameter = []
flash_y_parameter = []
arrow_x_parameter = []
arrow_y_parameter = []
for x1,y1,x2,y2 in zip(data['flash_episode_number'],data['flash_us_viewers'],data['arrow_episode_number'],data['arrow_us_viewers']):
flash_x_parameter.append([float(x1)])
flash_y_parameter.append(float(y1))
arrow_x_parameter.append([float(x2)])
arrow_y_parameter.append(float(y2))
return flash_x_parameter,flash_y_parameter,arrow_x_parameter,arrow_y_parameter
现在我们有了我们的参数,来写一个函数,用上面这些参数作为输入,给出一个输出,预测哪个节目会有更多观众。
Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Function to know which Tv show will have more viewers
def more_viewers(x1,y1,x2,y2):
regr1 = linear_model.LinearRegression()
regr1.fit(x1, y1)
predicted_value1 = regr1.predict(9)
print predicted_value1
regr2 = linear_model.LinearRegression()
regr2.fit(x2, y2)
predicted_value2 = regr2.predict(9)
#print predicted_value1
#print predicted_value2
if predicted_value1 > predicted_value2:
print "The Flash Tv Show will have more viewers for next week"
else:
print "Arrow Tv Show will have more viewers for next week"
把所有东西写在一个文件中。打开你的编辑器,把它命名为prediction.py,复制下面的代码到prediction.py中。
Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
# Required Packages
import csv
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import datasets, linear_model
# Function to get data
def get_data(file_name):
data = pd.read_csv(file_name)
flash_x_parameter = []
flash_y_parameter = []
arrow_x_parameter = []
arrow_y_parameter = []
for x1,y1,x2,y2 in zip(data['flash_episode_number'],data['flash_us_viewers'],data['arrow_episode_number'],data['arrow_us_viewers']):
flash_x_parameter.append([float(x1)])
flash_y_parameter.append(float(y1))
arrow_x_parameter.append([float(x2)])
arrow_y_parameter.append(float(y2))
return flash_x_parameter,flash_y_parameter,arrow_x_parameter,arrow_y_parameter
# Function to know which Tv show will have more viewers
def more_viewers(x1,y1,x2,y2):
regr1 = linear_model.LinearRegression()
regr1.fit(x1, y1)
predicted_value1 = regr1.predict(9)
print predicted_value1
regr2 = linear_model.LinearRegression()
regr2.fit(x2, y2)
predicted_value2 = regr2.predict(9)
#print predicted_value1
#print predicted_value2
if predicted_value1 > predicted_value2:
print "The Flash Tv Show will have more viewers for next week"
else:
print "Arrow Tv Show will have more viewers for next week"
x1,y1,x2,y2 = get_data('input_data.csv')
#print x1,y1,x2,y2
more_viewers(x1,y1,x2,y2)
可能你能猜出哪个节目会有更多观众——但运行一下这个程序看看你猜的对不对。
3) 替换数据集中的缺失值
有时候,我们会遇到需要分析包含有缺失值的数据的情况。有些人会把这些缺失值舍去,接着分析;有些人会用最大值、最小值或平均值替换他们。平均值是三者中最好的,但可以用线性回归来有效地替换那些缺失值。
这种方法差不多像这样进行。
首先我们找到我们要替换那一列里的缺失值,并找出缺失值依赖于其他列的哪些数据。把缺失值那一列作为Y_parameters,把缺失值更依赖的那些列作为X_parameters,并把这些数据拟合为线性回归模型。现在就可以用缺失值更依赖的那些列预测缺失的那一列。
一旦这个过程完成了,我们就得到了没有任何缺失值的数据,供我们自由地分析数据。
为了练习,我会把这个问题留给你,所以请从网上获取一些缺失值数据,解决这个问题。一旦你完成了请留下你的评论。我很想看看你的结果。
个人小笔记:
我想分享我个人的数据挖掘经历。记得在我的数据挖掘引论课程上,教师开始很慢,解释了一些数据挖掘可以应用的领域以及一些基本概念。然后突然地,难度迅速上升。这令我的一些同学感到非常沮丧,被这个课程吓到,终于扼杀了他们对数据挖掘的兴趣。所以我想避免在我的博客文章中这样做。我想让事情更轻松随意。因此我尝试用有趣的例子,来使读者更舒服地学习,而不是感到无聊或被吓到。
谢谢读到这里——请在评论框里留下你的问题或建议,我很乐意回复你。
E. python时间序列模型预测为什么时一条直线
python时间序列模型预测时一条直线是因为是线性模型的原因。线性模型也称作趋势模型,它表示一个时间序列可以用一条直线来表示。它的基本等式:以一个公司的销售总额为例,一开始的初始是5000,每隔一个时间步长增加2500。指数平滑法是时间序列分析方法中的一种。它是一种用于预测未来发展趋势的建模算法。它有三种不同形式:一次指数平滑法、二次指数平滑法、及三次指数平滑法。三种指数平滑法都要更新上一时间步长的计算结果,并使用当前时间步长的数据中包含的新信息。通过混合新信息和旧信息来实现。
F. python解决jacob迭代法求解线性方程组
题主好. 经典的 Jacobi 迭代算法如下:
设 A = D - E, 则 x = D^{-1}*b + D^{-1}*E*x = C + T*x
可以参考如下代码(复制代码后请注意缩进):
import numpy as np
def linalg_solve_jacobi(A, b, x0, max_it, tol=1.0e-7):
# 判断 A, b 的维数是否正确
if A.shape[0]!=A.shape[1] or A.shape[0]!= b.shape[0]:
raise Exception('A must be square or A and b must be compatible!')
D = np.diag(A)
# 判断对角元素是否含零
if np.amin(np.absolute(D)) < 1e-14:
raise Exception('Diagonal elements of A must not be zeros!')
# 设置初始迭代步数为0
n = 0
#
# Jacobi 算法:
# 设 A = D - E, 则 x = D^{-1}*b + D^{-1}*E*x = C + T*x
#
invD = np.diag(1.0/D) # invD is inv(D)
C = np.dot(invD,b) # C is inv(D)*b
T = np.dot(invD, np.diag(D)-A) # T is inv(D)*E
while n < max_it:
x = C + np.dot(T, x0)
if np.linalg.norm(x-x0)<tol:
break
x0[:] = x
n+=1
# 如果超过最大迭代步数, 迭代失败
if n>max_it:
raise Exception('Failed to converge within {} steps!'.format(max_it))
# 成功, 返回
return x, n
if __name__ == "__main__":
A = np.array([[10, -1, 2, 0],[-1, 11, -1, 3],
[2, -1, 10, -1], [0, 3, -1, 8]], dtype=np.float64)
b = np.array([6, 25, -11, 15], dtype=np.float64)
x0 = np.array([0,0,0,0], dtype=np.float64)
max_it = 100
tol = 1.0e-7
x,n=linalg_solve_jacobi(A, b, x0, max_it, tol)
print(x,n)
G. 根号a-+5的最小值和a的值
梯度下降是非常常用的优化算法。作为机器学习的基础知识,这是一个必须要掌握的算法。借助本文,让我们来一起详细了解一下这个算法。
前言
本文的代码可以到我的Github上获取:
https://github.com/paulQuei/gradient_descent
本文的算法示例通过Python语言实现,在实现中使用到了numpy和matplotlib。如果你不熟悉这两个工具,请自行在网上搜索教程。
关于优化
大多数学习算法都涉及某种形式的优化。优化指的是改变x以最小化或者最大化某个函数的任务。
我们通常以最小化指代大多数最优化问题。最大化可经由最小化来实现。
我们把要最小化或最大化的函数成为目标函数(objective function)或准则(criterion)。
我们通常使用一个上标*表示最小化或最大化函数的x值,记做这样:
[x^* = arg; min; f(x)]
优化本身是一个非常大的话题。如果有兴趣,可以通过《数值优化》和《运筹学》的书籍进行学习。
模型与假设函数
所有的模型都是错误的,但其中有些是有用的。– George Edward Pelham Box
模型是我们对要分析的数据的一种假设,它是为解决某个具体问题从数据中学习到的,因此它是机器学习最核心的概念。
针对一个问题,通常有大量的模型可以选择。
本文不会深入讨论这方面的内容,关于各种模型请参阅机器学习的相关书籍。本文仅以最简单的线性模型为基础来讨论梯度下降算法。
这里我们先介绍一下在监督学习(supervised learning)中常见的三个符号:
m,描述训练样本的数量
x,描述输入变量或特征
y,描述输出变量或者叫目标值
训练集会包含很多的样本,我们用 表示其中第i个样本。
x是数据样本的特征,y是其目标值。例如,在预测房价的模型中,x是房子的各种信息,例如:面积,楼层,位置等等,y是房子的价格。在图像识别的任务中,x是图形的所有像素点数据,y是图像中包含的目标对象。
我们是希望寻找一个函数,将x映射到y,这个函数要足够的好,以至于能够预测对应的y。由于历史原因,这个函数叫做假设函数(hypothesis function)。
学习的过程如下图所示。即:首先根据已有的数据(称之为训练集)训练我们的算法模型,然后根据模型的假设函数来进行新数据的预测。
线性模型(linear model)正如其名称那样:是希望通过一个直线的形式来描述模式。线性模型的假设函数如下所示:
[h_{ heta}(x) = heta_{0} + heta_{1} * x]
这个公式对于大家来说应该都是非常简单的。如果把它绘制出来,其实就是一条直线。
下图是一个具体的例子,即: 的图形:
在实际的机器学习工程中,你会拥有大量的数据。这些数据会来自于某个数据源。它们存储在csv文件中,或者以其他的形式打包。
但是本文作为演示使用,我们通过一些简单的代码自动生成了需要的数据。为了便于计算,演示的数据量也很小。
import numpy as np
max_x = 10
data_size = 10
theta_0 = 5
theta_1 = 2
def get_data:
x = np.linspace(1, max_x, data_size)
noise = np.random.normal(0, 0.2, len(x))
y = theta_0 + theta_1 * x + noise
return x, y
这段代码很简单,我们生成了x范围是 [1, 10] 整数的10条数据。对应的y是以线性模型的形式计算得到,其函数是:。现实中的数据常常受到各种因素的干扰,所以对于y我们故意加上了一些高斯噪声。因此最终的y值为比原先会有轻微的偏离。
最后我们的数据如下所示:
x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
y = [6.66, 9.11, 11.08, 12.67, 15.12, 16.76, 18.75, 21.35, 22.77, 24.56]
我们可以把这10条数据绘制出来这样就有一个直观的了解了,如下图所示:
虽然演示用的数据是我们通过公式计算得到的。但在实际的工程中,模型的参数是需要我们通过数据学习到的。所以下文我们假设我们不知道这里线性模式的两个参数是什么,而是通过算法的形式求得。
最后再跟已知的参数进行对比以验证我们的算法是否正确。
有了上面的数据,我们可以尝试画一条直线来描述我们的模型。
例如,像下面这样画一条水平的直线:
很显然,这条水平线离数据太远了,非常的不匹配。
那我们可以再画一条斜线。
我们初次画的斜线可能也不贴切,它可能像下面这样:
最后我们通过不断尝试,找到了最终最合适的那条,如下所示:
梯度下降算法的计算过程,就和这种本能式的试探是类似的,它就是不停的迭代,一步步的接近最终的结果。
代价函数
上面我们尝试了几次通过一条直线来拟合(fitting)已有的数据。
二维平面上的一条直线可以通过两个参数唯一的确定,两个参数的确定也即模型的确定。那如何描述模型与数据的拟合程度呢?答案就是代价函数。
代价函数(cost function)描述了学习到的模型与实际结果的偏差程度。以上面的三幅图为例,最后一幅图中的红线相比第一条水平的绿线,其偏离程度(代价)应该是更小的。
很显然,我们希望我们的假设函数与数据尽可能的贴近,也就是说:希望代价函数的结果尽可能的小。这就涉及到结果的优化,而梯度下降就是寻找最小值的方法之一。
对于每一个样本,假设函数会依据计算出一个估算值,我们常常用来表示。即 。
很自然的,我们会想到,通过下面这个公式来描述我们的模型与实际值的偏差程度:
[(h_ heta(x^i) - y^i)^2 = (widehat{y}^{i} - y^i)^2 = ( heta_{0} + heta_{1} * x^{i} - y^{i})^2]
请注意, 是实际数据的值, 是我们的模型的估算值。前者对应了上图中的离散点的y坐标,后者对应了离散点在直线上投影点的y坐标。
每一条数据都会存在一个偏差值,而代价函数就是对所有样本的偏差求平均值,其计算公式如下所示:
[L( heta) = frac {1}{m} sum_{i=1}^{m}(h_ heta(x^i) - y^i)^2 = frac {1}{m} sum_{i=1}^{m}( heta_{0} + heta_{1} * x^{i} - y^{i})^2]
当损失函数的结果越小,则意味着通过我们的假设函数估算出的结果与真实值越接近。这也就是为什么我们要最小化损失函数的原因。
借助上面这个公式,我们可以写一个函数来实现代价函数:
def cost_function(x, y, t0, t1):
cost_sum = 0
for i in range(len(x)):
cost_item = np.power(t0 + t1 * x[i] - y[i], 2)
cost_sum += cost_item
return cost_sum / len(x)
这个函数的代码应该不用多做解释,它就是根据上面的完成计算。
我们可以尝试选取不同的 和 组合来计算代价函数的值,然后将结果绘制出来:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
theta_0 = 5
theta_1 = 2
def draw_cost(x, y):
fig = plt.figure(figsize=(10, 8))
ax = fig.gca(projection='3d')
scatter_count = 100
radius = 1
t0_range = np.linspace(theta_0 - radius, theta_0 + radius, scatter_count)
t1_range = np.linspace(theta_1 - radius, theta_1 + radius, scatter_count)
cost = np.zeros((len(t0_range), len(t1_range)))
for a in range(len(t0_range)):
for b in range(len(t1_range)):
cost[a][b] = cost_function(x, y, t0_range[a], t1_range[b])
t0, t1 = np.meshgrid(t0_range, t1_range)
ax.set_xlabel('theta_0')
ax.set_ylabel('theta_1')
ax.plot_surface(t0, t1, cost, cmap=cm.hsv)
在这段代码中,我们对 和 各自指定了一个范围进行100次的采样,然后以不同的 组合对来计算代价函数的值。
如果我们将所有点的代价函数值绘制出来,其结果如下图所示:
从这个图形中我们可以看出,当 越接近 [5, 2]时其结果(偏差)越小。相反,离得越远,结果越大。
直观解释
从上面这幅图中我们可以看出,代价函数在不同的位置结果大小不同。
从三维的角度来看,这就和地面的高低起伏一样。最高的地方就好像是山顶。
而我们的目标就是:从任意一点作为起点,能够快速寻找到一条路径并以此到达图形最低点(代价值最小)的位置。
而梯度下降的算法过程就和我们从山顶想要快速下山的做法是一样的。
在生活中,我们很自然会想到沿着最陡峭的路往下行是下山速度最快的。如下面这幅图所示:
针对这幅图,细心的读者可能很快就会有很多的疑问,例如:
对于一个函数,怎么确定下行的方向?
每一步该往前走多远?
有没有可能停留在半山腰的平台上?
这些问题也就是本文接下来要讨论的内容。
算法描述
梯度下降算法最开始的一点就是需要确定下降的方向,即:梯度。
我们常常用 来表示梯度。
对于一个二维空间的曲线来说,梯度就是其切线的方向。如下图所示:
而对于更高维空间的函数来说,梯度由所有变量的偏导数决定。
其表达式如下所示:
[ abla f({ heta}) = ( frac{partial f({ heta})}{partial heta_1} , frac{partial f({ heta})}{partial heta_2} , ... , frac{partial f({ heta})}{partial heta_n} )]
在机器学习中,我们主要是用梯度下降算法来最小化代价函数,记做:
[ heta ^* = arg min L( heta)]
其中,L是代价函数,是参数。
梯度下降算法的主体逻辑很简单,就是沿着梯度的方向一直下降,直到参数收敛为止。
记做:
[ heta ^{k + 1}_i = heta^{k}_i - lambda abla f( heta^{k})]
这里有几点需要说明:
收敛是指函数的变化率很小。具体选择多少合适需要根据具体的项目来确定。在演示项目中我们可以选择0.01或者0.001这样的值。不同的值将影响算法的迭代次数,因为在梯度下降的最后,我们会越来越接近平坦的地方,这个时候函数的变化率也越来越小。如果选择一个很小的值,将可能导致算法迭代次数暴增。
公式中的 称作步长,也称作学习率(learning rate)。它决定了每一步往前走多远,关于这个值我们会在下文中详细讲解。你可以暂时人为它是一个类似0.01或0.001的固定值。
在具体的项目,我们不会让算法无休止的运行下去,所以通常会设置一个迭代次数的最大上限。
线性回归的梯度下降
有了上面的知识,我们可以回到线性模型代价函数的梯度下降算法实现了。
首先,根据代价函数我们可以得到梯度向量如下:
[ abla f({ heta}) = (frac{partial L( heta)}{ partial heta_{0}}, frac{ partial L( heta)}{ partial heta_{1}}) = (frac {2}{m} sum_{i=1}^{m}( heta_{0} + heta_{1} * x^{i} - y^{i}) , frac {2}{m} sum_{i=1}^{m}( heta_{0} + heta_{1} * x^{i} - y^{i}) x^{i})]
接着,将每个偏导数带入迭代的公式中,得到:
[ heta_{0} := heta_{0} - lambda frac{partial L( heta_{0})}{ partial heta_{0}} = heta_{0} - frac {2 lambda }{m} sum_{i=1}^{m}( heta_{0} + heta_{1} * x^{i} - y^{i}) heta_{1} := heta_{1} - lambda frac{partial L( heta_{1})}{ partial heta_{1}} = heta_{1} - frac {2 lambda }{m} sum_{i=1}^{m}( heta_{0} + heta_{1} * x^{i} - y^{i}) x^{i}]
由此就可以通过代码实现我们的梯度下降算法了,算法逻辑并不复杂:
learning_rate = 0.01
def gradient_descent(x, y):
t0 = 10
t1 = 10
delta = 0.001
for times in range(1000):
sum1 = 0
sum2 = 0
for i in range(len(x)):
sum1 += (t0 + t1 * x[i] - y[i])
sum2 += (t0 + t1 * x[i] - y[i]) * x[i]
t0_ = t0 - 2 * learning_rate * sum1 / len(x)
t1_ = t1 - 2 * learning_rate * sum2 / len(x)
print('Times: {}, gradient: [{}, {}]'.format(times, t0_, t1_))
if (abs(t0 - t0_) < delta and abs(t1 - t1_) < delta):
print('Gradient descent finish')
return t0_, t1_
t0 = t0_
t1 = t1_
print('Gradient descent too many times')
return t0, t1
这段代码说明如下:
我们随机选择了 都为10作为起点
设置最多迭代1000次
收敛的范围设为0.001
学习步长设为0.01
如果我们将算法迭代过程中求得的线性模式绘制出来,可以得到下面这幅动态图:
最后算法得到的结果如下:
Times: 657, gradient: [5.196562662718697, 1.952931052920264]
Times: 658, gradient: [5.195558390180733, 1.9530753071808193]
Times: 659, gradient: [5.194558335124868, 1.9532189556399233]
Times: 660, gradient: [5.193562479839619, 1.9533620008416623]
Gradient descent finish
从输出中可以看出,算法迭代了660次就收敛了。这时的结果[5.193562479839619, 1.9533620008416623],这已经比较接近目标值 [5, 2]了。如果需要更高的精度,可以将delta的值调的更小,当然,此时会需要更多的迭代次数。
高维扩展
虽然我们举的例子是二维的,但是对于更高维的情况也是类似的。同样是根据迭代的公式进行运算即可:
[ heta_{i} = heta_{i} - lambda frac {partial L( heta)}{partial heta_i} = heta_{i} - frac{2lambda}{m} sum_{i=1}^{m}(h_ heta(x^{k})-y^k)x_i^k]
这里的下标i表示第i个参数,上标k表示第k个数据。
梯度下降家族BGD
在上面的内容中我们看到,算法的每一次迭代都需要把所有样本进行遍历处理。这种做法称为之Batch Gradient Descent,简称BGD。作为演示示例只有10条数据,这是没有问题的。
但在实际的项目中,数据集的数量可能是几百万几千万条,这时候每一步迭代的计算量就会非常的大了。
于是就有了下面两个变种。
SGD
Stochastic Gradient Descent,简称SGD,这种算法是每次从样本集中仅仅选择一个样本来进行计算。很显然,这样做算法在每一步的计算量一下就少了很多。
其算法公式如下:
[ heta_{i} = heta_{i} - lambda frac {partial L( heta)}{partial heta_i} = heta_{i} - lambda(h_ heta(x^k)-y^k)x_i^k]
当然,减少算法计算量也是有代价的,那就是:算法结果会强依赖于随机取到的数据情况,这可能会导致算法的最终结果不太令人满意。
MBGD
以上两种做法其实是两个极端,一个是每次用到了所有数据,另一个是每次只用一个数据。
我们自然就会想到两者取其中的方法:每次选择一小部分数据进行迭代。这样既避免了数据集过大导致每次迭代计算量过大的问题,也避免了单个数据对算法的影响。
这种算法称之为Mini-batch Gradient Descent,简称MBGD。
其算法公式如下:
[ heta_{i} = heta_{i} - lambda frac {partial L( heta)}{partial heta_i} = heta_{i} - frac{2lambda}{m} sum_{i=a}^{a + b}(h_ heta(x^k)-y^k)x_i^k]
当然,我们可以认为SGD是Mini-batch为1的特例。
针对上面提到的算法变种,该如何选择呢?
下面是Andrew Ng给出的建议:
如果样本数量较小(例如小于等于2000),选择BGD即可。
如果样本数量很大,选择 来进行MBGD,例如:64,128,256,512。
下表是 Optimization for Deep Learning 中对三种算法的对比
方法准确性更新速度内存占用在线学习BGD好慢高否SGD好(with annealing)快低是MBGD好中等中等是
算法优化
式7是算法的基本形式,在这个基础上有很多人进行了更多的研究。接下来我们介绍几种梯度下降算法的优化方法。
Momentum
Momentum是动量的意思。这个算法的思想就是借助了动力学的模型:每次算法的迭代会使用到上一次的速度作为依据。
算法的公式如下:
[v^t = gamma v^{t - 1} + lambda abla f( heta) heta = heta - v_t]
对比式7可以看出,这个算法的主要区别就是引入了,并且,每个时刻的受前一个时刻的影响。
从形式上看,动量算法引入了变量 v 充当速度角色——它代表参数在参数空间移动的方向和速率。速度被设为负梯度的指数衰减平均。名称动量来自物理类比,根据牛顿运动定律,负梯度是移动参数空间中粒子的力。动量在物理学上定义为质量乘以速度。在动量学习算法中,我们假设是单位质量,因此速度向量 v 也可以看作是粒子的动量。
对于可以取值0,而是一个常量,设为0.9是一个比较好的选择。
下图是momentum算法的效果对比:
对原来的算法稍加修改就可以增加动量效果:
def gradient_descent_with_momentum(x, y):
t0 = 10
t1 = 10
delta = 0.001
v0 = 0
v1 = 0
gamma = 0.9
for times in range(1000):
sum1 = 0
sum2 = 0
for i in range(len(x)):
sum1 += (t0 + t1 * x[i] - y[i])
sum2 += (t0 + t1 * x[i] - y[i]) * x[i]
v0 = gamma * v0 + 2 * learning_rate * sum1 / len(x)
v1 = gamma * v1 + 2 * learning_rate * sum2 / len(x)
t0_ = t0 - v0
t1_ = t1 - v1
print('Times: {}, gradient: [{}, {}]'.format(times, t0_, t1_))
if (abs(t0 - t0_) < delta and abs(t1 - t1_) < delta):
print('Gradient descent finish')
return t0_, t1_
t0 = t0_
t1 = t1_
print('Gradient descent too many times')
return t0, t1
以下是该算法的输出:
Times: 125, gradient: [4.955453758569991, 2.000005017897775]
Times: 126, gradient: [4.955309381126545, 1.9956928964532015]
Times: 127, gradient: [4.9542964317327005, 1.9855674828684156]
Times: 128, gradient: [4.9536358220657, 1.9781180992510465]
Times: 129, gradient: [4.95412496254411, 1.9788858350530971]
Gradient descent finish
从结果可以看出,改进的算法只用了129次迭代就收敛了。速度比原来660次快了很多。
同样的,我们可以把算法计算的过程做成动态图:
对比原始的算法过程可以看出,改进算法最大的区别是:在寻找目标值时会在最终结果上下跳动,但是越往后跳动的幅度越小,这也就是动量所产生的效果。
Learning Rate 优化
至此,你可能还是好奇该如何设定学习率的值。
事实上,这个值的选取需要一定的经验或者反复尝试才能确定。
关键在于,这个值的选取不能过大也不能过小。
如果这个值过小,会导致每一次迭代的步长很小,其结果就是算法需要迭代非常多的次数。
那么,如果这个值过大会怎么样呢?其结果就是:算法可能在结果的周围来回震荡,却落不到目标的点上。下面这幅图描述了这个现象:
事实上,学习率的取值未必一定要是一个常数,关于这个值的设定有很多的研究。
下面是比较常见的一些改进算法。
AdaGrad
AdaGrad是Adaptive Gradient的简写,该算法会为每个参数设定不同的学习率。它使用历史梯度的平方和作为基础来进行计算。
其算法公式如下:
[ heta_i = heta_i - frac{lambda}{sqrt{G_t + epsilon}} abla f( heta)]
对比式7,这里的改动就在于分号下面的根号。
根号中有两个符号,第二个符号比较好理解,它就是为了避免除0而人为引入的一个很小的常数,例如可以设为:0.001。
第一个符号的表达式展开如下:
[G_t = sum_{i = 1}^{t} abla f( heta){i} abla f( heta){i}^{T}]
这个值其实是历史中每次梯度的平方的累加和。
AdaGrad算法能够在训练中自动的对learning rate进行调整,对于出现频率较低参数采用较大的学习率;相反,对于出现频率较高的参数采用较小的学习率。因此,Adagrad非常适合处理稀疏数据。
但该算法的缺点是它可能导致学习率非常小以至于算法收敛非常的慢。
关于这个算法的直观解释可以看李宏毅教授的视频课程:ML Lecture 3-1: Gradient Descent。
RMSProp
RMS是Root Mean Square的简写。RMSProp是AI教父Geoff Hinton提出的一种自适应学习率方法。AdaGrad会累加之前所有的梯度平方,而RMSProp仅仅是计算对应的平均值,因此可缓解Adagrad算法学习率下降较快的问题。
该算法的公式如下:
[E[ abla f( heta_{i})^2]^{t} = gamma E[ abla f( heta_{i})^2]^{t - 1} + (1-gamma)( abla f( heta_{i})^{t})^{2} heta_i = heta_i - frac{lambda}{sqrt{E[g^2]^{t+1} + epsilon}} abla f( heta_{i})]
类似的,是为了避免除0而引入。 是衰退参数,通常设为0.9。
这里的 是t时刻梯度平方的平均值。
Adam
Adam是Adaptive Moment Estimation的简写。它利用梯度的一阶矩估计和二阶矩估计动态调整每个参数的学习率。
Adam的优点主要在于经过偏置校正后,每一次迭代学习率都有个确定范围,使得参数比较平稳。
该算法公式如下:
[m^{t} = eta_{1} m^{t-1} + (1-eta_{1}) abla f( heta) v^{t} = eta_{2} v^{t-1} + (1-eta_{2}) abla f( heta)^2 widehat{m}^{t} = frac{m^{t}}{1 - eta^{t}_1} widehat{v}^{t} = frac{v^{t}}{1 - eta^{t}_2} heta = heta - frac{lambda}{sqrt{widehat{v}^{t}} + epsilon}widehat{m}^{t}]
,分别是对梯度的一阶矩估计和二阶矩估计。, 是对,的校正,这样可以近似为对期望的无偏估计。
Adam算法的提出者建议 默认值为0.9,默认值为0.999,默认值为 。
在实际应用中 ,Adam较为常用,它可以比较快地得到一个预估结果。
优化小结
这里我们列举了几种优化算法。它们很难说哪种最好,不同的算法适合于不同的场景。在实际的工程中,可能需要逐个尝试一下才能确定选择哪一个,这个过程也是目前现阶段AI项目要经历的工序之一。
实际上,该方面的研究远不止于此,如果有兴趣,可以继续阅读 《Sebastian Ruder: An overview of gradient descent optimization algorithms》 这篇论文或者 Optimization for Deep Learning 这个Slides进行更多的研究。
由于篇幅所限,这里不再继续展开了。
算法限制
梯度下降算法存在一定的限制。首先,它要求函数必须是可微分的,对于不可微的函数,无法使用这种方法。
除此之外,在某些情况下,使用梯度下降算法在接近极值点的时候可能收敛速度很慢,或者产生Z字形的震荡。这一点需要通过调整学习率来回避。
另外,梯度下降还会遇到下面两类问题。
局部最小值
局部最小值(Local Minima)指的是,我们找到的最小值仅仅是一个区域内的最小值,而并非全局的。由于算法的起点是随意取的,以下面这个图形为例,我们很容易落到局部最小值的点里面。
这就是好像你从上顶往下走,你第一次走到的平台未必是山脚,它有可能只是半山腰的一个平台的而已。
算法的起点决定了算法收敛的速度以及是否会落到局部最小值上。
坏消息是,目前似乎没有特别好的方法来确定选取那个点作为起点是比较好的,这就有一点看运气的成分了。多次尝试不同的随机点或许是一个比较好的方法,这也就是为什么做算法的优化这项工作是特别消耗时间的了。
但好消息是:
对于凸函数或者凹函数来说,不存在局部极值的问题。其局部极值一定是全局极值。
最近的一些研究表明,某些局部极值并没有想象中的那么糟糕,它们已经非常的接近全局极值所带来的结果了。
鞍点
除了Local Minima,在梯度下降的过程中,还有可能遇到另外一种情况,即:鞍点(Saddle Point)。鞍点指的是我们找到点某个点确实是梯度为0,但它却不是函数的极值,它的周围既有比它小的值,也有比它大的值。这就好像马鞍一样。
如下图所示:
多类随机函数表现出以下性质:在低维空间中,局部极值很普遍。但在高维空间中,局部极值比较少见,而鞍点则很常见。
不过对于鞍点,可以通过数学方法Hessian矩阵来确定。关于这点,这里就不再展开了,有兴趣的读者可以以这里提供的几个链接继续探索。
参考资料与推荐读物
Wikipeida: Gradient descent
Sebastian Ruder: An overview of gradient descent optimization algorithms
吴恩达:机器学习
吴恩达:深度学习
Peter Flach:机器学习
李宏毅 - ML Lecture 3-1: Gradient Descent
PDF: 李宏毅 - Gradient Descent
Intro to optimization in deep learning: Gradient Descent
Intro to optimization in deep learning: Momentum, RMSProp and Adam
Stochastic Gradient Descent – Mini-batch and more
刘建平Pinard - 梯度下降(Gradient Descent)小结
多元函数的偏导数、方向导数、梯度以及微分之间的关系思考
[Machine Learning] 梯度下降法的三种形式BGD、SGD以及MBGD
H. 如何用Python进行线性回归以及误差分析
数据挖掘中的预测问题通常分为2类:回归与分类。
简单的说回归就是预测数值,而分类是给数据打上标签归类。
本文讲述如何用Python进行基本的数据拟合,以及如何对拟合结果的误差进行分析。
本例中使用一个2次函数加上随机的扰动来生成500个点,然后尝试用1、2、100次方的多项式对该数据进行拟合。
拟合的目的是使得根据训练数据能够拟合出一个多项式函数,这个函数能够很好的拟合现有数据,并且能对未知的数据进行预测。
代码如下:
importmatplotlib.pyplot as plt
importnumpy as np
importscipy as sp
fromscipy.statsimportnorm
fromsklearn.pipelineimportPipeline
fromsklearn.linear_modelimportLinearRegression
fromsklearn.
fromsklearnimportlinear_model
''''' 数据生成 '''
x = np.arange(0,1,0.002)
y = norm.rvs(0, size=500, scale=0.1)
y = y + x**2
''''' 均方误差根 '''
defrmse(y_test, y):
returnsp.sqrt(sp.mean((y_test - y) **2))
''''' 与均值相比的优秀程度,介于[0~1]。0表示不如均值。1表示完美预测.这个版本的实现是参考scikit-learn官网文档 '''
defR2(y_test, y_true):
return1- ((y_test - y_true)**2).sum() / ((y_true - y_true.mean())**2).sum()
''''' 这是Conway&White《机器学习使用案例解析》里的版本 '''
defR22(y_test, y_true):
y_mean = np.array(y_true)
y_mean[:] = y_mean.mean()
return1- rmse(y_test, y_true) / rmse(y_mean, y_true)
plt.scatter(x, y, s=5)
degree = [1,2,100]
y_test = []
y_test = np.array(y_test)
fordindegree:
clf = Pipeline([('poly', PolynomialFeatures(degree=d)),
('linear', LinearRegression(fit_intercept=False))])
clf.fit(x[:, np.newaxis], y)
y_test = clf.predict(x[:, np.newaxis])
print(clf.named_steps['linear'].coef_)
print('rmse=%.2f, R2=%.2f, R22=%.2f, clf.score=%.2f'%
(rmse(y_test, y),
R2(y_test, y),
R22(y_test, y),
clf.score(x[:, np.newaxis], y)))
plt.plot(x, y_test, linewidth=2)
plt.grid()
plt.legend(['1','2','100'], loc='upper left')
plt.show()
该程序运行的显示结果如下:
[ 0. 0.75873781]
rmse=0.15, R2=0.78, R22=0.53, clf.score=0.78
[ 0. 0.35936882 0.52392172]
rmse=0.11, R2=0.87, R22=0.64, clf.score=0.87
[ 0.00000000e+00 2.63903249e-01 3.14973328e-01 2.43389461e-01
1.67075328e-01 1.10674280e-01 7.30672237e-02 4.88605804e-02
......
3.70018540e-11 2.93631291e-11 2.32992690e-11 1.84860002e-11
1.46657377e-11]
rmse=0.10, R2=0.90, R22=0.68, clf.score=0.90
I. 马科维茨有效前沿python求出每个点的配置比例
1.马科维茨有效前沿中每个点的配置比例可以通过求解其凸组合来确定。这需要解决一个线性规划问题,目标是最大化有效前沿上的点到要优化的点的距离,约束条件是各点的权重和为1,且每个权重大于等于0。通过求解该线性规划问题,可以得到每个点在有效前沿上的配置比例。
2.深入分析
2.1 根据马科维茨有效前沿的定义,其上每个点可以由多个极点通过凸组合得到。配置比例就是各极点在该凸组合中所占的权重。这些权重满足二次型约束:w1+w2+...+wn=1,wi≥0,i=1,2,...,n。
2.2 求解配置比例的关键在于构建一个线性规划模型。目标函数设为maximizeρ,其中ρ代表有效前沿上点到要优化的点的欧几里得距离。约束条件为wi≥0,w1+w2+...+wn=1。通过求解该线性规划问题,可以得到最优的权重配置,这些权重值即为各极点在有效前沿点上的配置比例。
2.3 上述线性规划问题可以通过python中的凸优化库cvxopt来求解。要先构建线性规划问题的矩阵形式,再使用cvxopt.solvers.lp这个函数进行求解。函数输入为目标函数矩阵、约束矩阵和变量下界上界,输出为最优化权重向量,这即为所求的配置比例。
2.4 求解配置比例需要先确定马科维茨有效前沿,这需要使用极小化方法来寻找要优化的目标函数的极小点。常用方法有梯度下降法、Newton法以及interior point method等。通过这些方腊禅法可以找到目标函数的所有极小点,构建出有效前沿,这为后续的配置比例计算提供了必要的条件。
3.建议
3.1 在马科维茨有效前沿的计算中,应采用既定的优化方法,如牛顿法,来确保找到全局最优解。这有助于构建出完备的有效前沿,为后续配置比例计算提供准确的计算基础。
3.2 线性规划建模时,目标函数和约束条件应表达清晰准确。各矩阵应事先规范化,以避免由于数据量级差异导致的计算误差。
3.3 凸优化库的选择上,推荐使用经过验证的优化库,如cvxopt。这类库运算速度较快,且可以直接求解various 类型的凸规划问题,避免由于算法实现带来的误差。
3.4 配置比例的计算结果还需要进行正确性验证。可以通过计算有效前沿上各点的凸组合,与原有效前沿点的坐标进行比较,看其误差是否在可接受范围内。这一验证过程是保山毕证最终计算结逗局芹果正确的必要步骤。
J. python 线性模型
因变量是你自己确定的,一般主成分得分是作为自变量的,叫主成分回归分析