深度学习完全初学者指南:人工神经网络
深度学习入门！在 15 分钟内攻克人工神经网络的基础知识
本文是 深度学习 系列完整初学者指南的一部分。
Image by Lolame via Pixabay
什么是深度学习？
这是从例子中学习。差不多就是这么回事！
在非常基础的层面上，深度学习是一种机器学习技术。它教会计算机通过层层过滤输入，以学习如何预测和分类信息。观察可以是图像、文本或声音的形式。
深度学习的灵感来自人脑过滤信息的方式。它的目的是模仿人类大脑如何工作来创造一些真正的魔法。
GIF via GIPHY
深度学习试图模仿新大脑皮层中神经元层的活动。
这是一个名副其实的人工神经网络。
在人脑中，大约有 1000 亿个神经元。每个神经元都与大约 10 万个相邻的神经元相连。这就是我们试图创造的，但是在某种程度上是为机器服务的。
Image by geralt on Pixabay
这对神经元、轴突、树突等意味着什么？神经元有一个体，树突和一个轴突。来自一个神经元的信号沿着轴突传递到下一个神经元的树突。信号传递的连接被称为突触。
Image by mohamed_hassan on Pixabay
神经元本身是没用的。但是当你有很多这样的人时，他们会一起创造出一些真正的奇迹。这就是深度学习算法背后的想法！你从观察中获得输入，然后把你的输入放到一个层中。该层创建一个输出，该输出又成为下一层的输入，依此类推。这种情况反复发生，直到你的最终输出信号！
因此神经元(或节点)获得一个或多个信号(输入值)，这些信号通过神经元。该神经元传递输出信号。把输入层想象成你的感官:例如，你看到的、闻到的和感觉到的东西。这些是一次观察的独立变量。这些信息被分解成数字和计算机可以使用的二进制数据位。(您需要标准化或规范化这些变量，以便它们在相同的范围内。)
那么突触呢？每个突触都被分配了权重，这对人工神经网络至关重要。重量是人工神经网络学习的方式。通过调整权重，人工神经网络决定信号传递的程度。当你训练你的网络时，你决定如何调整权重。
人工神经网络是如何学习的？
有两种不同的方法可以让程序做你想做的事情。首先，有特别指导和硬编程的方法。在这种方法中，您确切地告诉程序您想要它做什么。然后还有神经网络。在神经网络中，你告诉你的网络输入和你想要的输出，让它自己学习。通过允许网络自己学习，我们可以避免输入所有规则的必要性。对于神经网络，你可以创建架构，然后让它去学习。一旦它训练好了，你可以给它一个新的图像，它将能够区分输出。
Photo by Annie Spratt on Unsplash
有不同种类的神经网络。一般分为前馈和反馈网络。
前馈网络是一个包含输入、输出和隐藏层的网络。信号只能单向传播(向前)。输入数据传递到执行计算的层。每个处理元件基于其输入的加权和进行计算。新值成为新的输入值，用于下一层(前馈)。这将贯穿所有层，并决定输出。前馈网络通常用于例如数据挖掘。
一个反馈网络(例如，一个递归神经网络)有反馈路径。这意味着它们可以使用环路双向传输信号。神经元之间所有可能的连接都是允许的。由于这种类型的网络中存在环路，因此它成为一个非线性动态系统，不断变化直到达到平衡状态。反馈网络通常用于优化问题，其中网络寻找相互关联的因素的最佳安排。
大多数现代深度学习架构都是基于人工神经网络(ann)的。它们使用多层非线性处理单元进行特征提取和转换。每个后续层使用前一层的输出作为其输入。他们所学的形成了概念的层次结构。在这个层次结构中，每一级都学会将其输入数据转换成越来越抽象和复合的表示。
Image by ahmedgad on Pixabay
这意味着，例如，对于一幅图像，输入可能是像素矩阵。第一层可以对边缘进行编码并组成像素。下一层可能构成边缘的排列。下一层可能编码鼻子和眼睛。下一层可能会识别出图像包含人脸，等等。
神经元内部发生了什么？输入节点接收数字形式的信息。信息以激活值的形式呈现，每个节点都有一个编号。数字越高，激活程度越大。
基于连接强度(权重)和传递函数，激活值传递到下一个节点。每个节点对其接收的激活值求和(它计算加权和)，并基于其传递函数修改该和。接下来，它应用一个激活函数。激活函数是应用于这个特定神经元的函数。由此，神经元知道是否需要传递信号。激活通过网络运行，直到它到达输出节点。然后输出节点以我们可以理解的方式给我们信息。您的网络将使用成本函数来比较输出和实际预期输出。模型性能由代价函数来评价。它表示为实际值和预测值之间的差值。您可以使用许多不同的成本函数，查看网络中的误差。你在努力减少功能损失。(本质上，损失函数越低，就越接近你想要的输出)。信息返回，神经网络开始学习，目标是通过调整权重来最小化成本函数。这个过程被称为反向传播。
有兴趣了解更多关于成本函数的知识吗？查看 神经网络中使用的成本函数列表，以及堆栈交换上的应用
在前向传播中，信息被输入到输入层，并通过网络向前传播，以获得我们的输出值。我们将这些值与我们的预期结果进行比较。接下来，我们计算误差并反向传播信息。这允许我们训练网络并更新权重。反向传播允许我们同时调整所有的权重。在这个过程中，由于算法的结构，你可以同时调整所有的权重。这使您可以看到神经网络中的每个权重对错误的哪一部分负责。
渴望更多？你可能想读一下 Yann LeCun 等人的 高效反向投影 ，以及 Michael Nielsen 的 神经网络和深度学习 。
当您将重量调整到最佳水平后，您就可以进入测试阶段了！
什么是加权和？
神经元的输入可以是来自训练集的特征，也可以是来自前一层神经元的输出。两个神经元之间的每个连接都有一个独特的突触，并附有独特的权重。如果你想从一个神经元到下一个神经元，你必须沿着突触行进，并支付“通行费”(重量)。然后，神经元将激活函数应用于来自每个传入突触的加权输入的总和。它将结果传递给下一层的所有神经元。当我们谈论更新网络中的权重时，我们谈论的是调整这些突触上的权重。
一个神经元的输入是前一层所有神经元的加权输出之和。每个输入都乘以与将输入连接到当前神经元的突触相关联的权重。如果在前一层中有 3 个输入或神经元，则当前层中的每个神经元将具有 3 个不同的权重:每个突触一个。
那么什么是激活函数呢？
简而言之，一个节点的激活函数定义了该节点的输出。
激活函数(或传递函数)将输入信号转换成输出信号。它在 0 到 1 或-1 到 1 的范围内映射输出值。这是一个抽象概念，代表细胞内动作电位发放的速率。这是一个数字，代表细胞被激发的可能性。最简单的，函数是二元的:是(神经元触发)或否(神经元不触发)。输出可以是 0 或 1(开/关或是/否)，也可以是范围内的任何值。例如，如果您使用映射范围在 0 和 1 之间的函数来确定图像是一只猫的可能性，输出 0.9 将显示您的图像实际上是一只猫的概率为 90%。
Photo by minanafotos on Pixabay
我们有什么选择？有许多激活函数，但以下是四个非常常见的函数:
- 阈值函数这是一个阶跃函数。如果输入的合计值达到某个阈值，该函数将传递 0。如果它等于或大于零，那么它将传递 1。这是一个非常严格，简单，是或不是的函数。
Example threshold function
- Sigmoid 函数:该函数用于逻辑回归。与阈值函数不同，它是一个从 0 到 1 的平滑渐进过程。它在输出层非常有用，并且大量用于线性回归。(线性回归是统计学和机器学习中最知名的算法之一)。
Example sigmoid function
- 双曲正切函数这个函数非常类似于 sigmoid 函数。与从 0 到 1 的 sigmoid 函数不同，该值从-1 到 1 变到 0 以下。虽然这不是生物学中发生的事情，但这个函数在训练神经网络时会给出更好的结果。神经网络有时会在使用 sigmoid 函数进行训练时“卡住”。当有很多强烈的负面输入使输出接近于零时，就会发生这种情况，这会扰乱学习过程。
Example hyperbolic tangent function (tanh)
- 整流函数这可能是神经网络宇宙中最流行的激活函数。这是最有效的，也是生物学上最合理的。尽管它有一个扭结，但在 0°扭结后，它是平滑和渐变的。举例来说，这意味着你的输出要么是“否”，要么是“是”的百分比这个函数不需要标准化或其他复杂的计算。
Example rectifier function
想潜得更深吗？查看 Xavier Glorot 等人的 深度稀疏整流器神经网络
比如说，你想要的值是二进制的。你在寻找一个“是”或“否”。你想使用哪个激活功能？从上面的例子中，你可以使用阈值函数，或者你可以使用 sigmoid 激活函数。sigmoid 函数可以给你一个肯定的概率。
那么，权重究竟是如何调整的呢？
你可以使用蛮力方法来调整权重，并测试数千种不同的组合。即使是最简单的只有五个输入值和一个隐藏层的神经网络，你也会得到 10⁷⁵可能的组合。在世界上最快的超级计算机上运行这个程序需要的时间比迄今为止宇宙存在的时间还要长。
Photo by skorchanov on Pixabay
然而，如果您使用梯度下降**，您可以查看权重的倾斜角度，并找出它是正还是负，以便继续向下倾斜，在您寻求达到全局最小值的过程中找到最佳权重。**
如果你使用梯度下降**，你可以查看重物的倾斜角度，并找出它是正还是负。这允许你继续下坡，在你寻求达到全局最小值的过程中找到最佳权重。**
****梯度下降是一种求函数最小值的算法。你会一遍又一遍地看到类似的情况，有人被困在山顶，试图下来(找到最小值)。大雾使她看不清路，所以她使用梯度下降法到达山脚。她看了看她所在的山的陡峭程度，然后向最陡的下坡方向走去。你应该假设陡度不会立即变得明显。幸运的是，她有一个工具可以测量陡度。不幸的是，这个工具需要很长时间。她想尽可能少地使用它，以便在天黑前下山。真正的困难是选择她使用工具的频率，这样她就不会偏离轨道。在这个类比中，人就是算法。山的陡度是该点误差面的斜率。她走的方向是该点误差曲面的梯度。她使用的工具是微分(误差表面的斜率可以通过对该点的平方误差函数求导来计算)。在进行另一次测量之前，她行进的速率是算法的学习速率。这不是一个完美的类比，但它让你很好地理解了梯度下降是怎么回事。机器正在学习模型应该采用的梯度或方向，以减少误差。
随机梯度下降
梯度下降要求代价函数是凸的，但如果不是呢？
正常梯度下降将卡在局部最小值而不是全局最小值，从而导致网络不合格。在正常梯度下降中，我们将所有行插入同一个神经网络，看一看权重，然后调整它们。这被称为批量梯度下降。在随机梯度下降中，我们一行接一行，运行神经网络，查看成本函数，调整权重，然后移动到下一行。实际上，您正在调整每一行的权重。
随机梯度下降有更高的波动，这允许你找到全局最小值。它被称为“随机的”,因为样本是随机混合的，而不是作为一个单独的组或当它们出现在训练集中时。它看起来可能会慢一些，但实际上更快，因为它不必将所有数据加载到内存中，并在数据一起运行时等待。批量梯度下降的主要优势在于它是一种确定性算法。这意味着，如果你有相同的起始重量，每次运行网络，你会得到相同的结果。随机梯度下降总是随机工作。(您也可以运行小批量梯度下降，您可以设置行数，一次运行那么多行，然后更新您的权重。)
已经提出并使用了对基本随机梯度下降算法的许多改进，包括隐式更新(ISGD)、动量法、平均随机梯度下降、自适应梯度算法(AdaGrad)、均方根传播(RMSProp)、自适应矩估计(Adam)等等。
喜欢这个吗？你可能想看看安德鲁·特拉斯克 和迈克尔·尼尔森的 神经网络和深度学习
因此，这里有一个用随机梯度下降训练人工神经网络的快速演练:
- 1:将权重随机初始化为接近 0 的小数字
- 2:将数据集的第一个观测值输入到输入图层中，每个要素位于一个输入结点中。
- 3: 前向传播——从左到右，神经元的激活方式是每个神经元的激活受到权重的限制。你传播激活，直到你得到预期的结果。
- 4:比较预测结果和实际结果，测量产生的误差。
- 5: 反向传播 —从右到左，误差反向传播。根据权重对误差的影响程度来更新权重。(学习率决定了我们更新权重的程度。)
- 6: 强化学习(重复步骤 1-5，每次观察后更新权重)或 批量学习(重复步骤 1-5，但只在一批观察后更新权重)。
- 7:当整个训练集已经通过 ANN 时，这是一个时期。重复更多的纪元。
就是这样！你现在知道了在人工神经网络中发生的事情背后的基本思想！
Photo by Sam Mathews on Unsplash
还跟我在一起？ 过来看第三部 ！
和往常一样，如果你对这些信息做了什么酷的事情，请在下面的回复中留下评论，或者随时联系 LinkedIn@ annebonnerdata！
辛普森悖论导论
辛普森悖论是一个迷人的现象，说明了因果关系在推理中的重要性。如果你还没有读过这本书，你应该去看看为什么的书——这是这十年来最重要的书之一。
我认为辛普森悖论很好地介绍了因果推理的价值。在这篇文章中，我将通过几个例子来阐明它是如何产生的，以及我们如何开始形式化我们对它的思考。
击球率
肯·罗斯是辛普森悖论的一个常见例子。他指出，当观察德瑞克·基特和大卫·贾斯蒂斯的击球率时，一个奇怪的矛盾出现了:
David Justice has a better batting average every year, but Derek Jeter has a better batting average overall!
德瑞克·基特每个赛季的平均打击率都更差，但总的来说却有更好的平均打击率！
如果我们盯着这一点看一会儿，随着一些事情的突出，原因就变得很清楚了。第一年，德里克和大卫的分数都很低；第二年他们都得了高分。最重要的是，德里克的低分仅仅来自于几次击球，因此对他的总体平均得分影响不大。另一方面，大卫恰恰相反——他的较低分数占了他击球尝试的绝大部分。即使大卫在没有加权的两年中做得稍好，当我们把蝙蝠的这种不平衡与它们在分数上类似的巨大差异结合起来时，我们得到了这个结果。辛普森悖论告诉我们，总体统计数据的适当加权平均值可能与我们在每个子组中看到的趋势相矛盾。
下面是另一个简单的例子——这次是图形化的。比方说，非常年轻的人(往往比成年人矮)经常打篮球，而且打得很好。比方说，年纪大的人(往往更高)停止练习，通常实际上很不擅长投篮。
现在，假设我们测量他们的投篮表现，希望最终证明高个子更擅长篮球。我们得到的是这样的:
A graphical example of Simpson’s paradox. Taller kids and taller adults are both better at basketball, but in general being tall makes you worse!
又一个难题！对于每个年龄段的人来说，个子高意味着你是一个更好的篮球运动员。总的来说，个子高会让你变得更糟！
我们都知道个子高会让你变得更优秀，解释这个结果的正确方法是用年龄来分隔，而不是用总体人口。但是为什么以及我们是如何知道的呢？我们如何将经验法则推广到更复杂的情况？
因果推理
这里有一个更难的例子，它反映了现实生活中的一项医学研究。假设你生病了，给你开了一种药，对以前的病人有如下效果:
Uh oh. The drug is good for men, good for women, and bad for people.
你没看错。如果你是男性，这种药物会增加你康复的机会；如果你是女性，这也会增加你康复的几率。但如果你是一个人，它更有可能伤害你。你吃这种药吗？
值得庆幸的是，上述说法在逻辑上是不连贯的，是错误的。但是很容易阅读数据并以这种方式看待它。我们如何解决这个问题？
答案
我们知道，这里有三个因素给了我们辛普森悖论。有(1)患者的性别，(2)药物的有效性，以及(3) 患者有多大可能属于治疗或非治疗组。请记住，这是组的不平衡，使整体平均加权可能不同于每个子组的趋势。
那么问题就变成了，因果关系在哪里？可以肯定地说，你的性别不会因为你的治疗或康复而改变。另一方面，性别似乎会影响你选择治疗的可能性——女性在 75%的情况下选择治疗，而男性只有 12%。
鉴于我们以这种方式检查数据，我们清楚地相信性别可能会影响恢复。最后(也很明显)，治疗会影响你康复的机会，而不是相反。我们可以画一个漂亮的因果图来总结我们所有的因果知识。
A simple causal diagram of what’s going on. Gender alters the likelihood of accepting treatment. We also expect the likelihood of recovery to be a function of both receiving treatment AND gender.
哪里出了问题？希望很清楚，图表中“受性别影响的治疗选择”部分——将性别与治疗联系起来的箭头——是导致我们悖论的原因。
我们不希望那个信息改变我们的决定。一个人选择的可能性有多大？这种治疗并不能告诉我们这种治疗最终对他们有多有效。为了准确了解这种药物对整个人类的疗效，我们必须消除这种联系。换句话说，我们应该平均每个组的百分比，而不是根据大小来加权。这消除了两性治疗选择不均衡的影响。
Simpson’s paradox is resolved!
它概括了吗？
在上面的例子中，信任子组而不是整体平均值是有意义的。其他情况下，可以反过来。
举一个例子(这个例子来自 Judea Pearl 的书)，我们研究了一种用于预防心脏病发作的药物，并评估了它对低血压和高血压患者的影响。我们看到了类似的问题——这种药物减少了低血压和高血压患者的心脏病发作，但总体上增加了心脏病发作！我们该怎么办？在这种情况下，当我们绘制因果图时，它可能看起来像这样:
Causal Diagram for our second example. In this case, the blood pressure is a mediating variable in the path from treatment to heart attack, and thus we should trust the overall statistic.
这看起来非常类似于我们的第一个例子，但它不是。仔细看。我们小组的位置和治疗方法都变了。
在第一个病例中，我们的亚组(性别)影响了接受治疗的。在第二个例子中，我们的亚组(血压)受到接受治疗的的影响。
在这种情况下，心脏病发作的风险和患者的血压都是治疗的因果下游，所以我们应该而不是按血压细分。血压是治疗发挥作用的一种机制——它是从治疗到心脏病发作风险之间因果关系的中介变量。
与之前的案例形成对比，在之前的案例中，治疗和恢复之间的关系受到性别的影响。性别是治疗选择和恢复的上游；为了了解治疗和康复之间的真正关系，我们不得不以此为条件。当处理没有图形双亲时，我们不需要辅助变量。当它发生时(并且该母公司有一个“后门路径”来影响利息)，我们会这样做。
另一种方法
思考辛普森悖论的另一种方式是通过纯统计的语言。首先，让我们稍微绕一圈，回顾一下基础知识。
联合概率可以很容易地分解成条件的乘积。这里有一个有趣的例子——你有了一只新狗，你对它的个性很感兴趣！有些狗很年轻，有些很聪明(有些不聪明)，有些则很有活力。
下面是代表狗的性格概率分布的方框。我们把正方形的面积设为 1。放下一个别针，它所代表的狗就会有下面的性格特征。
在所有可能的狗中，你可以看到一些是年轻的，一些是聪明的，一些是精力充沛的，还有一些是三者兼而有之。一只随机的狗具有一些特定特征组合的可能性可以很容易地被分解——下面是一个简单的例子来说明当我们寻找年轻、聪明和精力充沛的狗时(Y=1，I=1，E=1)。
Factoring a joint probability into the product of conditionals.
为什么你会想把一个关节放入这样的条件句中呢？在辛普森悖论的背景下，我们将通过重温上面的第一个例子来回答这个问题。
在那个例子中，我们在看(T)治疗、(G)性别和®康复的联合概率。我们想询问在给定其他性别和治疗组合的情况下，R=1(成功康复)的可能性。让我们将关节分解，这样我们就可以使用它:
All possible ways to rewrite that joint as a product of conditioned variables.
有六种方法可以分解三个变量的结合。然而，在这种情况下，我们试图用观察关系(在这种情况下，可以说‘假设某人成功康复，他是一个男人的概率’之类的话)，而不是用询问他们之间的因果关系(‘如果我们强迫治疗的概率等于 1，一个男人的康复会怎么样？’).T5 这里我们介绍珍珠的 do-operator 。
代替 p(R|T)，一个关于治疗和恢复之间相关性的观察，我们实际上关心 p(R |do(T))——给定我们因果干预和强制治疗(通过意志力设定 T=1)的恢复概率。do 运算符限制了我们如何分解关节，因为我们必须将 do(T)保持在关于恢复的条件线的右侧，以回答关于治疗如何影响恢复的因果问题。
所以我们对 p(R|G，do(T)) 感兴趣。这给了我们两个选择:
在第二个选择中，我们发现 p(G|do(T))是一个无意义的关系。我们现在认为是一种因果关系，这个术语回避了一个问题“如果你接受治疗，你的性别改变的可能性有多大？”如果当我们改变(T)时性别不变，我们就不能在因果关系中连贯地使用这样一个术语。
另一方面，第一个等式没有逻辑缺陷。我们注意到 p(do(T)|G)正好是 1——我们在任何理性条件下的 do(T)概率将是 1，因为这是 do 运算符的定义。我们可以去掉这个术语，简化我们的表述。
另一种说法是，基于性别接受治疗的可能性是而不是我们想要的影响我们对治疗成功概率的评估的东西，所以这个术语必须删除。我们知道这个术语是坏消息，因为我们有因果直觉。我们的条件必须遵循我们的因果本能，匹配我们的因果图，使得子节点的概率总是(并且仅仅)取决于它们的父节点。
我们只剩下 p(R | do(T)，G) * p(G) 。它向我们展示了康复的可能性与该性别人群的比例有关，而不是与总人口有关。这通过一个更正式的过程给了我们和我们的因果图一样的答案。
现实生活，伙计
辛普森悖论可能出现在现实生活中——就像加州大学伯克利分校众所周知的有争议的性别偏见问题。然而，它真正的美妙之处在于阐释和介绍了统计学中因果推理背后的基本原理。
我们所做的一切在标准统计学的背景下都没有意义，但确实有意义。我希望这篇关于辛普森悖论的介绍让你对如何开始形式化你的因果逻辑有所体会。
长期以来，统计学缺乏讨论因果推理和关系的语言。谢天谢地，统计世界似乎正在积极转变。随之而来的将是科学的一大福音。
辛普森悖论
理解统计错觉
辛普森悖论是一个封闭的术语，指的是在全球层面分析时，群体层面的趋势被隐藏或逆转的情况。这些都是需要思考的 fum 例子，是混合效果模型的主要候选对象，这也是我上一篇文章的主题。
说明
下面是一个来自维基百科文章的例子:蓝色和红色组有明显的积极趋势，然而用一条线回归这两个组会产生消极趋势。
在分析医疗保健索赔趋势时，通常会出现类似的情况(当然也可以扩展到其他行业)。下图显示了 4 种医疗服务——两种住院服务，由于需要多日住院，因此费用更高，还有两种门诊服务。一项常见的健康计划计划旨在将低强度/低成本的住院服务转移到门诊服务中:这样做可以为会员节省大量的自付额/共同保险费用。这产生了一个有趣的趋势，因为相对于其他低强度的门诊服务，低强度的住院服务通常强度更高。
这两个部分分别代表每年每千人的利用率(分析利用率时调整成员增长/衰退的常用指标)和单位成本(通常称为“成本/混合”,因为字段的聚合可以隐藏基础服务利用率的变化，而单位成本不会发生任何变化。
随着图表从第 1 期(P1)移动到第 2 期(P2)，50 个单位的廉价住院治疗转移到门诊，在那里它们相对更贵。每项服务的单位成本不变。总的来说，这 50 个单位现在的价格是 9000 美元，而以前的价格是 12000 美元，但是潜在的趋势是骗人的。由于每个时段中高成本服务的比例较高，这两种单位成本趋势都有所增加。住院病人利用趋势下降，被门诊病人抵消。在 IP/OP 级别，该计划看起来像一场灾难！总体水平与现实相符:利用率没有变化，但一些服务现在成本更低了。
技术
如果你发现辛普森悖论的存在，就很容易解决它。在我的混合效应模型文章中，我们通过测试数据的同质性和分析残差图来确定这一点。在上面的维基百科图片中，你可以注意到，随着 x 轴从左向右移动，每组的残差从负值变为正值——残差中的这种相关性是一个好迹象，表明你错过了一些东西。
这强调了数据可视化和解释性数据分析的重要性。我们可以在这篇文章的结尾指出，虽然一张图片胜过千言万语，但有些图片需要千言万语的解释——这就是数据科学家和分析专业的价值。
在 PyBullet 中模拟 ML 的图像——快速简单的方法
当将深度强化学习(RL)应用于机器人时，我们面临一个难题:当深度学习需要数十万甚至数百万个样本时，我们如何训练机器人来完成一项任务？为了在从未见过的物体上实现 96%的抓取成功率，谷歌和伯克利的研究人员通过 580，000 次真实世界的抓取尝试训练了一个机器人代理。这项壮举花费了七个机器人和几个星期来完成。如果没有谷歌资源，对于一般的 ML 从业者来说，这似乎是没有希望的。我们不能指望使用物理机器人轻松地运行数十万次迭代训练，这容易磨损，需要人工监督，两者都不便宜。如果我们能够预先训练这种 RL 算法，以大幅减少现实世界中所需的尝试次数，这将更加可行。
随着深度学习的出现，RL 技术已经成熟，但对数据的需求也越来越大。为了弥合这一差距，许多研究人员正在探索训练数据的合成生成，利用 3D 渲染技术来制作任务环境的模型。虽然这种技术在模拟环境中产生了奇迹，但它并不能很好地推广到现实世界。任何与深度学习密切合作的人都知道，利用训练数据中的细微差别来“欺骗”任务是多么有效。在实践中，对于这些经过模拟训练的模型来说，真实世界是样本外的，因此不出所料，它们会失败。
最近的研究致力于减少“现实差距”，即真实世界和 3D 渲染传真之间的显著差异，以在模拟中预先训练机器人代理，从而大幅减少所需的真实世界训练量。通过积极地随机化模拟的外观和动态，模型学习理论上应该推广到现实世界的特征。根据最近的研究，这甚至可以将训练机器人所需的训练次数减少 99% 。对于产生数十万甚至数百万训练图像不可行的计算机视觉任务，这种方法非常有吸引力。
这个概念
比方说，我们希望生成合成训练数据来预训练一个分类器，该分类器可以区分一组小对象。我们希望选择一个对象，将其放入一个虚拟的“垃圾箱”，并呈现一个屏幕截图。理想情况下，我们希望随机重复这个过程数千次，以创建一个丰富的对象图像数据集。为此(受到 OpenAI 和 T2 谷歌所做工作的启发)，我最近开始了为自己的应用程序生成训练数据的旅程。这让我想到了 PyBullet ，这是一个基于子弹物理 SDK 为机器人和机器学习应用程序设计的 Python 模块。
与 Maya 或 Blender 等其他 3D 渲染解决方案不同，PyBullet 专注于机器人技术，并具有关节、动力学模拟、正向和反向运动学等概念的原生实现。此外，它可以很容易地与包管理器一起安装，允许您将其他 Python 包，如 NumPy 和 TensorFlow，无缝集成到您的模拟中。虽然类似于专有的物理模拟软件 MuJoCo ，PyBullet 既免费又易于安装，因此对于希望进行模拟和机器人实验的人来说，它是一个很好的选择。本指南旨在为那些希望生成训练图像的人提供介绍，但是官方 PyBullet 快速入门指南可以在这里找到。请注意，本教程假设您使用的是 Python 版本 3.6 以上(如果您使用的是旧版本，请尝试为本教程创建一个虚拟环境)。
装置
首先，我们将在 Python 环境中安装 PyBullet。
pip install numpy
pip install pybullet
请注意，首先安装 NumPy 是可选的，但在渲染图像时建议这样做，因为在 C/C++和 Python 之间复制图像缓冲区的开销很大。
入门指南
现在您已经成功安装了所有的东西，让我们直接进入 PyBullet 会话。PyBullet 依赖于客户端-服务器模型，其中您的 Python 会话从模拟服务器发送和接收数据。主要有两种类型的服务器:直接和图形用户界面。顾名思义，GUI 允许您在 GUI 中查看物理模拟。“直接”对于无标题渲染非常有用，可用于在没有 GPU 的情况下进行高效渲染。当我们准备渲染成千上万的图像作为训练数据时，直接将是选择的模式，但现在我们将选择 GUI。让我们从初始化 PyBullet 服务器开始。首先启动一个 python 交互式会话，然后输入以下内容:
import pybullet as pb
physicsClient = pb.connect(pb.GUI)
这一步是必要的，因为它实例化了您将使用的世界模拟。physicsClient 变量包含该服务器的唯一 ID。这使得 PyBullet 可以运行多个服务器，甚至跨多台机器，并通过跟踪这些 id 从单个脚本管理它们。PyBullet 中的许多实体都用它们的 id 来指代，包括对象、碰撞器和纹理。
装载飞机
现在，您应该会看到一个模拟窗口弹出。太好了！这意味着你的模拟在起作用。您应该会看到一个主视口，其中包含用于 RGB、深度和分段数据的附加窗口。世界应该是完全空的。让我们首先为我们的模拟创建一个平面。这架飞机带有 pybullet_data 包，可以很容易地加载到任何模拟中。
import pybullet_data
pb.setAdditionalSearchPath(pybullet_data.getDataPath())
planeId = pb.loadURDF(‘plane.urdf’)
这有几个作用。首先，它加载另一个名为 pybullet_data 的模块，其中包含许多有用东西的占位符 URDFs。这被添加到路径中，这样当我们最终调用 loadURDF 时，它知道在哪里可以找到‘plane . urdf’并将其加载到我们的场景中。现在我们应该有一个平面横跨我们的场景。
加载要渲染的对象
PyBullet 中最受支持的文件格式是 URDF、SDF 和 MJCF。这些格式很容易加载和配置，并有自己的专门加载功能。然而，我们在网上找到的许多 3D 模型都是简单而流行的波前 OBJ 文件格式。加载这些是可能的，但是需要知道更多关于 PyBullet 中模型的具体细节，所以我们将咬紧牙关学习用 OBJ 文件加载我们自己的模型。
我们将使用程序生成的对象数据集进行模拟任务。这套包括 1000 个 OBJ 文件，其中包括相应的对撞机，材料和 URDF 版本。在这里下载那些，并把它们解压到你的项目文件夹中的一个名为“random_urdfs”的文件夹中。
物理模拟中的每个典型 3D 对象都由一个描述对象出现方式的视觉网格和一个描述对象物理尺寸以及它将如何与世界交互的碰撞器网格组成，例如它如何与其他物理对象碰撞。碰撞器网格，或简称为碰撞器，通常是我们视觉网格的简化表示，用它来计算碰撞要便宜得多。同时使用这两种类型的网格使我们的物理模拟能够高效运行，同时仍然具有视觉吸引力。
视觉形状
让我们从使用 createVisualShape 函数加载一个对象的可视形状开始。我们可以从程序生成的对象数据集中的一个对象开始。
visualShapeId = pb.createVisualShape(
shapeType=pb.GEOM_MESH,
fileName='random_urdfs/000/000.obj',
rgbaColor=None,
meshScale=[0.1, 0.1, 0.1])
第一个参数 shapeType 告诉 PyBullet 我们正在加载哪种形状。这些包括球体、长方体、圆柱体、平面、胶囊和网格。在我们的例子中，我们希望加载 OBJ 文件中描述的自定义网格。您还可以描述 RGB 颜色(作为描述红色、绿色、蓝色和 alpha 通道的长度为 4 的列表)，我们在这里将其设置为 None 。我们稍后将替换对象的纹理。最后，OBJ 文件本质上没有单位，所以我们可能需要指定一个网格比例，以确保我们的对象不会太大或太小。
碰撞形状
现在我们有了一个视觉形状来参考我们的对象，但是没有物理形状。让我们使用 createCollisionShape 函数来加载它。
collisionShapeId = pb.createCollisionShape(
shapeType=pb.GEOM_MESH,
fileName='random_urdfs/000/000_coll.obj',
meshScale=[0.1, 0.1, 0.1])
类似于 createVisualShape 函数，我们需要使用 shapeType 指定我们想要的网格，以及在哪里找到带有文件名的网格数据。我们可以也应该指定网格比例，这样碰撞器就有了和视觉形状一样的形状。
多体
在这一点上，我们有一个视觉形状和一个碰撞器的参考，但这两个没有连接，我们还没有对象。这就是创建多体的用武之地。这个函数将这两者结合起来，并创建一个对象的实例。
multiBodyId = pb.createMultiBody(
baseMass=1.0,
baseCollisionShapeIndex=collisionShapeId,
baseVisualShapeIndex=visualShapeId,
basePosition=[0, 0, 1],
baseOrientation=pb.getQuaternionFromEuler([0, 0, 0]))
我们将碰撞器的 id 和可视形状传递给 createMultiBody 来创建一个新的物体。这个物体需要一个位置和方向来知道它在物理世界中的位置，还需要一个质量来知道它应该如何与其他物体交互。用长度为 3 的 x、y、z 坐标列表设定位置。为了简单起见，我们现在将质量设置为 1.0。
方向有点复杂，需要我们在一个长度为 4 的四元数列表中指定。四元数是三维空间中旋转的一种有用的四维表示，在 Blender、Maya 和 Unity 等计算机图形软件中非常重要，用于避免旋转伪影，如万向节锁。PyBullet 有一个方便的函数，叫做getquaternonfromeuler，可以让我们用绕 x、y、z 轴的旋转来描述旋转，人类更容易可视化和理解。这种旋转的表示被称为欧拉角。
注意:如果收到类似“不支持的纹理图像格式[random_urdfs/000/000.png]”的消息，那么忽略它。无论如何，我们将在下一步覆盖默认纹理。
将纹理应用到对象
成功！我们已经生成了一个物理对象，并将其放置在我们的世界中。回想一下之前我们想要随机化 3D 模型的外观和动态。让我们从外观开始，在这里下载并提取可描述纹理数据集。这包含了 47 个纹理类别的成千上万的图像。将这些内容提取到一个名为“dtd”的文件夹中。让我们随机选择一个纹理，并将其应用到我们的模型中。
import os, glob, randomtexture_paths = glob.glob(os.path.join('dtd', '**', '*.jpg'), recursive=True)
random_texture_path = texture_paths[random.randint(0, len(texture_paths) - 1)]
textureId = pb.loadTexture(random_texture_path)
pb.changeVisualShape(multiBodyId, -1, textureUniqueId=textureId)
此代码片段递归地获取文件夹“dtd”中所有纹理的路径列表，扩展名为“”。jpg”并随机抽取一个样本。然后将它传递给函数 loadTexture 来加载纹理并创建一个 ID。然后，它使用 changeVisualShape 函数将带有 textureId 的纹理应用于带有 multiBodyId 的对象。现在你的物体应该已经应用了纹理。如果愿意，也可以应用相同的方法将纹理附加到平面上。
注意-1 对于必需的参数 jointIndex 是必需的，该参数指定了您想要应用纹理的链接。链接对于创建具有运动部件层次的对象非常有用，例如具有许多关节的机器人。这是 PyBullet 中的一个重要概念，但超出了本教程的范围，所以我们提供了一个值“-1”来告诉它将它应用于基本对象而不是任何链接。
开始我们的模拟
到目前为止，我们已经创建了我们的模拟，但没有什么是移动的。这有几个原因。首先，我们没有对场景应用任何重力，所以我们的对象是漂浮在空间中的。让我们尽快做到这一点。
pb.setGravity(0, 0, -9.8)
这给场景中的物体施加了-9.8 米/秒的向下负加速度，就像我们在地球表面经历的一样。然而，你的物体仍然没有像你预期的那样下落。这是因为模拟没有步进。
有两种方法可以逐步完成你的模拟:手动调用 stepSimulation 进行一个时间步长的模拟，或者调用 setRealTimeSimulation 进行实时渲染。出于演示的目的，我们将使用后者，但是一旦我们开始在一个循环中渲染一批对象，你可能会希望使用 stepSimulation 。
pb.setRealTimeSimulation(1)
设置为“1”可启用，设置为“0”可禁用。现在你的物体应该下落并落地。顺便说一下，默认的时间步长是 1/240 秒，但是如果需要，您可以使用 setTimeStep 函数来覆盖这个默认值。
如何渲染你的图像
我们现在有一个方法来加载一个平面，在我们的平面上放置一个 OBJ 文件，对对象应用纹理，并渲染场景。拼图的最后一块是渲染我们的对象的图像。本质上，我们希望在我们的对象上放置一个虚拟相机，将其向下指向我们的对象，并拍摄照片。
无论我们是要解决分类还是分割任务，PyBullet 都通过 getCameraImage 函数使我们的模拟拍照变得相当简单，该函数同时渲染一个 RGB 图像、一个分割蒙版和一个深度缓冲区。然而，我们需要首先指定相机的外部和内部属性，这个过程比看起来要复杂一些。
相机外在属性本质上描述了我们的相机在世界中的位置，比如它的方向和位置。相机固有属性描述相机本身的属性，例如视野(FOV)和其传感器的纵横比(例如 4:3、16:9 等。).在 PyBullet 中，我们分别用视图矩阵和投影矩阵来描述这些属性。PyBullet 快速入门指南链接到一个很棒的资源，可以了解更多关于相机的外在和内在属性及其对应矩阵的信息。
视图矩阵
摄像机视图矩阵是一个复杂的 4x4 矩阵，但用最简单的术语来说，它描述了摄像机的位置和指向。有一些有用的辅助函数可以通过更直接地指定位置和旋转来创建这个矩阵。函数 computeViewMatrix 可以创建这个矩阵来交换三个向量。
viewMatrix = pb.computeViewMatrix(
cameraEyePosition=[0, 0, 3],
cameraTargetPosition=[0, 0, 0],
cameraUpVector=[0, 1, 0])
顾名思义， cameraEyePosition 描述了摄像机在 x 、 y 、 z 坐标中的物理位置。我们将把它放在平面正上方 3 个单位处，以我们的对象为中心。接下来我们描述我们希望相机面对的点 cameraTargetPosition 。我们想让它直视我们的物体，所以我们让它指向原点。最后，我们需要描述摄像机的方向。我们用camera uvector来做这个，它是一个指向我们相机顶部的向量。我们沿着 y 轴指出这一点。我们的相机现在应该指向正下方，这样它渲染的图像平面就与我们生成的平面完全平行了。
投影矩阵
投影矩阵很像视图矩阵，可以使用一些辅助函数来创建。在这种情况下，computeProjectionMatrixFOV函数以最简单、最贴切的方式描述了我们相机的固有属性。
projectionMatrix = pb.computeProjectionMatrixFOV(
fov=45.0,
aspect=1.0,
nearVal=0.1,
farVal=3.1)
fov 参数以度为单位指定摄像机 fov。它描述了相机的视野有多“宽”,原则上类似于真实相机的焦距。纵横比参数描述了摄像机的纵横比，为了简单起见，我们在这里将其设置为 1.0。 nearVal 和 farVal 参数分别描述了摄像机渲染对象的最小和最大距离。由于我们的相机在平面上方 3 个单位，我们将 farVal 设置为 3.1，以包含该平面。这些参数在 OpenGL 文档中有定义，你可以在这里阅读关于它们的更多信息。
getCameraImage
我们已经描述了相机的属性，现在准备渲染一些图像。让我们用 getCameraImage 来做这件事。
该函数返回三个图像缓冲区:rgbImg、depthImg 和 segImg。 rgbImg 是一个 uint8 图像，带有相机视觉效果的红色、绿色、蓝色和 alpha 通道。depthImg 是一种浮点灰度图像，描述单个渲染像素与相机的距离。它可以用于模拟真实世界深度传感器的视野，如 Microsoft Kinect。最后， segImg 是图像的分割掩模，其中每个像素包含具有对象 id 的唯一整数。这些对于训练机器人智能体的分割算法来说是非常宝贵的，例如需要识别物体以分类到各自的箱子中的机械臂，或者需要识别行人、街道标志和道路的无人驾驶汽车。
width, height, rgbImg, depthImg, segImg = pb.getCameraImage(
width=224,
height=224,
viewMatrix=viewMatrix,
projectionMatrix=projectionMatrix)
一旦调用了这个函数，您应该会在查看器的三个窗口中看到您渲染的图像。如果是这样的话，那么恭喜你！您已经渲染了您的第一个对象。您可能希望将这些图像缓冲区转换为 NumPy 矢量，并用图像库(如 PIL 或 imageio)保存它们。
结论
PyBullet 使我们创建虚拟场景以生成训练图像的任务变得相当简单。利用本教程中讨论的概念，您应该能够生成自己的数据集用于分类或分割。额外的增强可能值得探索，例如向相机和渲染对象的位置和方向添加噪声。人们也可能希望用你希望训练的真实世界环境的模型来代替飞机。希望这篇教程是你开发自己的机器学习物理模拟的坚实起点。
模拟故障
一百万次泵故障为什么不再多一次？如何模拟制造领域的故障？
Photo by Elijah O’Donnell on Unsplash
在制造业领域，我最常被问到的一个问题是，这台机器什么时候会出故障。在化学工程领域，机器常常是泵。许多过程和可靠性工程师根据行业图表和多年的泵工作经验来估算。但是很少有人能告诉我，在我们的工作区域内，泵的正常寿命是多少。因为这些泵不会经常出故障。
这就引出了一个问题，“我怎样才能更好地理解这种泵模式？”
在这种情况下，我经常使用的算法是蒙特卡罗模拟。正如您将看到的，这种风格将帮助我们建立一个基线，让我们更好地了解什么时候可能会发生失败。
我们将要谈论的:
- 如何建立蒙特卡洛模拟？
- 如何模拟泵故障(或任何有明确数据的过程)？
- 如何获得工程团队的信任？
让我们开始吧。下面我们设置了一些虚拟数据来输入我们的算法。
import numpy as np# Generate positive portions of the distribution. positive_num = np.random.randint(low=0, high=100, size = 200)#Generate negative portions of the distribution.negative_num = np.random.randint(low=-100, high=0, size = 300)#Joining the positive and negative datasets together.pump_pressure_changes = np.concatenate([negative_num, positive_num], axis = 0)
在上面的代码中，我们使用 numpy 生成了一些随机数据。这个数据集的负值比正值多。您可以将这些点中的每一点视为泵压的正变化和负变化。事实上，在现实世界的问题中，你应该只取泵压读数的第一个差值，就可以得到与你的问题相似的数据。
现在让我们来构建工作马函数。我们将一步步地了解这段代码的每一部分。
####################################################################
#Monte Carlo Simulation
####################################################################
def Base_Monte_Carlo(current_pressure, final_pressure, max_threshold, iterations, data):
hour_distribution = np.empty(shape=(1,0))
for i in np.arange(iterations):
print("Iteration number: {}".format(i))
hours_til_failure = 0
psi_left = current_pressure
not_done = True
while(not_done):
sample_chunk = np.random.choice(data)
print('Sample drawn: {}'.format(sample_chunk))
hours_til_failure += 1
psi_left += sample_chunk
print('psi left {}'.format(psi_left))
# Check if the predefined cutoff has been hit
if psi_left <= final_pressure:
# This is where the Late-Stopping or hope would be added.
hour_distribution = np.append(hour_distribution, hours_til_failure)
not_done = False
if psi_left > max_threshold:
print("print had a spike reseting back to normal")
psi_left = max_threshold
print(psi_left)
return(hour_distribution)
在此功能中，参数的含义如下
current_pressure =你的泵现在所处的压力，或者你想要模拟的压力。
最终压力=泵出现故障或在产品中产生不良结果的压力。例如，在化学领域，这可能意味着你的产品开始燃烧，因为它在管道中移动变慢了。在真正的解决方案中，这是你需要向熟悉你的流程的人询问的事情。
max_threshold =机械有物理可能的极限。例如，如果你将一杯水倒入水槽，排水管很容易跟上，但如果你将消防水管放入同一个水槽，水就会溢出。我们必须尽可能模拟这些物理障碍，这就是上限的原因。
迭代次数=这是我们想要运行模拟的次数。
现在让我们来看看这个函数在做什么。
hour_distribution = np.empty(shape=(1,0))
这一行代码只是给我们一个地方来存储每次模拟达到“失败”状态所花费的时间。在这个例子中，我们假设样本中的每个数据点都是以一个小时为间隔采集的。
for i in np.arange(iterations):
print("Iteration number: {}".format(i))
hours_til_failure = 0
psi_left = current_pressure
not_done = True
这段代码是模拟开始的地方。它会为您想要的每次迭代生成一个模拟。然后，它让我们知道我们在打印语句的模拟中进行了多远。
hours_til_failure 和 psi_left 是我们在每次迭代开始时重置的两个变量。hours_til_failure 实际上记录了我们的泵在每次迭代中运行了多长时间而没有出现故障。psi_left 是一个变量，我们稍后将根据 final_pressure 对其进行验证，以查看我们是否满足故障标准。
while(not_done):
sample_chunk = np.random.choice(data)
print('Sample drawn: {}'.format(sample_chunk))
hours_til_failure += 1
psi_left += sample_chunk
print('psi left {}'.format(psi_left))
这里，当我们的“not_done”变量为真时，这段代码将从我们的数据集中随机抽取一个样本，给泵已经运行的时间加上一个小时，并从泵的当前压力中加上或减去适当的量。
if psi_left <= final_pressure:
hour_distribution = np.append(hour_distribution, hours_til_failure)
not_done = False
if psi_left > max_threshold:
print("print had a spike resetting back to normal")
psi_left = max_threshold
第一个 if 语句检查我们是否低于我们的 final_pressure。如果是这种情况，那么模拟运行的小时数被附加到 hour_distribution，模拟重新开始。
第二个 if 语句检查我们的压力是否超过了允许的最大值。如果发生这种情况，它会将压力调整回可能的范围。
让我们试着用一些虚构的压力和限制来运行 1000 次。
distribution = Base_Monte_Carlo(current_pressure=1000,
final_pressure=200,
max_threshold=1500,
iterations=1000,
data=pump_pressure_changes)
最后，该函数返回我们的泵在所有迭代中持续的小时分布。我们可以使用 matplotlib.pyplot hist 函数对此进行图形化。
from matplotlib import pyplot as plt
# Plotting the histogram.
plt.hist(x = distribution, bins = 100)
plt.xlabel('Survival time in hours')
plt.ylabel('Number of occurences')
在本例中，我们可以看到泵在运行 90 小时左右最频繁地达到其故障状态。然而，我们也看到，假设它曾经运行了将近 500 个小时。
最后一点要注意的是，一定要从你的数据集中去掉任何异常事件，例如，机器的启动，或者像自然灾害这样的外部事件引起的意外停机。
现在，我发现当试图为过程工程师或可靠性工程师运行这些模型时，通常最好从更广泛的变量开始作为模型的输入。例如，如果他们告诉你，“我们知道在每小时 200 磅的压力下，我们的化学物质开始燃烧”。将您的模型设置为 150 失败，这样，如果您在不同的操作条件下遇到挑战，您可以说您对机器的运行时间能力持乐观态度。一旦建立了信任，您就可以将模型磨练成精确的运行参数。
我们在本帖中讨论了以下内容:
- 如何建立蒙特卡洛模拟？
- 如何模拟泵故障(或任何有明确数据的过程)？
- 如何获得工程团队的信任？
有很多有趣的方法，你可以调整这个模型，你可以真实地模拟任何改变状态的过程。另一个使用它的例子是一个进程的能量消耗。同样，这是一个基线模型，有助于理解该过程。
我希望听到你的反馈和建设性的想法。
用张量网络模拟量子电路
奇点研究
用物理学改进机器学习模型
在本文中，我们将讨论用于模拟量子电路的开源软件 Quimb 。Quimb 如此有趣和有用的原因之一是它能够用张量网络进行计算。如果你不熟悉张量网络，我们会在这里谈一点，但你也可以看看我在 Medium " 黑洞机器学习"上的另一篇文章，在那里我讨论了张量网络在物理学中的历史，以及它们在机器学习和人工智能中的用途。如果你仍然好奇并喜欢钻研一点代码，你也可以看看我的 Quimb 教程的 Github repo。
量子电路
Quantum Circuit using IBM’s Quantum Circuit Composer
量子电路是与音乐家使用的五线谱非常相似的图表。你从左到右读它们，就像音乐一样，五线谱上的每一个方框都代表了在被称为门的量子位上执行的一个矩阵操作。这些线代表量子位在时间中运行，通过这些门。在电路的最右侧，我们通常会看到小的测量符号，类似于汽车中的速度计，或者自行车打气筒上显示气压的仪表。它们表示量子位的测量值，并及时返回该量子位的状态向量。由于量子物理的概率性质，这些测量通常要进行多次。每次测量都将返回不同的值，有可能处于两种状态之一。这是量子计算机与经典计算机真正不同的一点。它们天生具有概率性，量子位不像笔记本电脑或手机中的经典二进制位那样简单地以 0 或 1 的形式存在。在上面由 IBM 的量子计算电路合成器制造的量子电路中，我们看到一个哈达玛门，由一个带有“**H”**的小盒子表示，它对单个量子位 q[0]进行操作。接下来是一个’ CNOT’ 门，它对前两个量子位 q[0]和 q[1]进行操作。电路中的第三个是三个量子位 q[1]、q[2]和 q[3]上的算子，由一个 8x8 矩阵表示。最后，对量子位 q[1]进行测量。你可以在 IBM 注册一个免费账户，在云中的实际 IBM 量子计算机上创建和运行你自己的量子电路。
张量网络
假设你有一个可以在经典计算机上有效模拟的量子电路(这并不总是可能的，因为在许多情况下经典计算机很难模拟计算)。有效模拟一些量子电路的一种方法是使用张量网络。张量网络是张量的图形表示，可以认为是多维数组。一个秩为 0 的张量仅仅是一个标量。秩为 1 的张量是一个矢量。一个秩为 2 的张量是一个矩阵，以此类推到秩为 3 后变得难以可视化的更高维数组。
谷歌的 TensorNetwork 库
Google’s TensorNetwork Library Logo
谷歌最近发布了一个以其著名的 TensorFlow 作为后端运行的库。由于 TensorFlow 的流行和广泛使用，在 TensorFlow 上建立一个库是一个巨大的优势，但 TensorNetwork 库本身目前并不是最用户友好的，对新手来说可能有点令人生畏。如果您还不熟悉 TensorFlow，并且还没有扎实的编程背景，那么在尝试执行基本张量网络构造和操作之外的任何操作时，对抛出的一些错误进行故障排除可能会令人头疼。据说这个库在 Github 页面上是稳定的，但根据代码作者的说法，它在很大程度上仍处于 alpha 版本。取决于在这个项目上投入了多少努力，这个库可能还需要一段时间才能像 TensorFlow 一样在研究之外被广泛采用。
一个用户友好的选择…Quimb
Quimb Logo
Quimb 是一个开源软件，用于执行量子电路仿真和构建张量网络，在张量网络上执行操作，使用张量网络开发算法，以及为机器学习训练张量网络。它为那些已经熟悉量子力学和量子计算的人提供了出色的文档，对于那些不熟悉的人，我正在创建一个 Github 教程库，解释代码背后的数学和物理知识，并提供许多额外的示例和解释。第一个例子，查看第一个教程的笔记本。Quimb 更出色的地方在于，它使用了 NetworkX ，一个用于创建、操作和研究复杂网络的结构、动态和功能的 Python 包，来创建漂亮的张量网络图可视化:
Tensor Network using NetworkX in Quimb
张量网络和机器学习
那么，除了量子物理学之外，张量网络还有什么用途呢？它们可以在机器学习中用作深度神经网络中的层。这弥合了经典机器学习和量子机器学习新领域之间的差距，量子机器学习使用量子变分电路作为量子神经网络。在许多情况下，使用张量网络模拟这种量子神经网络可以显著提高机器学习算法的性能。因此，即使你不打算很快钻研量子计算，如果你正在使用神经网络和机器学习，你仍然可以通过将它们集成到你当前的机器学习模型中来受益于量子计算技术的使用。
使用 sklearn 在 Python 中模拟(复制)R 回归图
当谈到数据科学和机器学习工作负载时，R 和 Python 是最流行和最强大的语言。Python 通常被视为具有易于理解的语法的通用语言，而 R 则用于统计分析，有大约 12000 个可用的包。
有几十篇文章从主观角度比较了 Python 和 R。在这里，我不会厚此薄彼。
在本帖中，我们将讨论使用 sklearn 在 Python 中对 R 回归图的复制。
大多数 R 的功能可以很容易地直接在 Python 中转换，但是如果不使用自定义函数，有些功能很难找到对等的功能。例如，R 中用于回归模型的 plot () 函数在 Python 中并不具有所有绘图的直接等价物。
让我们用 r 中可用的忠实数据集来讨论一下，该数据集包含两个变量 喷发 (喷发时间分钟)和*(等待下一次喷发的时间)的 272 个观测值。该数据集揭示了美国怀俄明州黄石国家公园老忠实间歇泉的喷发间隔时间和喷发持续时间。*
现在，让我们将 R 中的简单线性回归模型拟合到该数据集，以根据喷发时间预测等待时间。
对于验证线性假设，上述每一个图都有其自身的重要性。我们不打算在这里深究。让我们关注 Python 代码来拟合相同的线性回归模型。
导入所有必要的库并加载所需的数据。
现在让我们使用 sklearn.linear_model 在 faithful 数据集上拟合一个线性回归模型。
现在让我们使用 sklearn 逐一关注所有的回归图。
- 剩余剧情
这是由 R 中的 plot() 函数生成的第一个图，有时也被称为残差与拟合图。通过绘制拟合值和残差之间的散点图，它有助于验证线性假设。
如果该图描述了任何特定或规则的模式，则假设目标变量和预测值之间的关系在本质上是非线性的，即存在非线性。并且曲线中没有图案是所选特征和目标变量之间线性的标志。
使用 Seaborn 中的 residplot()函数可以获得 Python 中的相同绘图。这里，第一个和第二个参数分别指向拟合(预测)值和目标变量。lowess=True 确保绘制 lowess(平滑)回归线，并使用 line_kws 参数自定义该线的属性。
2。QQ 剧情
该图描述了残差(误差)是否正态分布。如果这些点靠近法线，那么残差被假设为正态分布。
在 Python 中，使用 seaborn 中可用的 probplot()函数可以实现相同的绘图。这里，残差作为参数传递给函数。
3。比例位置图
通常，它用于猜测残差的同方差。它是标准化残差与拟合值的平方根图。如果它没有描述特定的模式，那么拟合的回归模型支持同质性假设。
使用 Seaborn 中的 regplot()函数可以获得 Python 中的相同绘图。这里，第一个和第二个参数分别指向拟合值和平方根标准化残差。
4。杠杆图
通常，它用于猜测异常值对回归拟合的影响。目前，我不知道如何在 Python 中为基于 sklearn 的合身模型绘制相同的内容。一旦我能够弄清楚，我会更新相同的。使用 statsmodels 很容易使用内置的杠杆图绘制，但我不打算在这里讨论它。
如果你已经知道如何为基于 sklearn 的拟合模型绘制相同的图。在评论里告诉我，我会加进去的！
使用几何布朗运动在 Python 中模拟股票价格
对几何布朗运动的动力学和主要驱动因素的高级描述，并附有示例 Python 代码
The picture was taken from @StockMarketMeme
介绍
如果我们能准确预测股票价格在近期或远期的变化，那就太好了。我们会很富有，但几乎不可能做出准确的预测。股票价格的变动涉及太多难以建模的因素。人类心理学就是其中之一。当然，投资者是根据经验证据和股市指标做出决策的。然而，他们还是人类。不同的人会对数据做出不同的解读，这是因为他们的风险偏好和当时的情绪。当一个国家的突发新闻导致一个投资者买入一只股票时，它会导致另一个投资者卖出同样的股票。因此，预测股票价格是一项困难的工作，但我们仍然有有价值的工具，可以帮助我们了解股票价格在某种程度上的运动。
在本文中，我们将讨论如何使用 Python 构建几何布朗运动(GBM)模拟。在构建脚本时，我们还探索了 GBM 模型背后的直觉。我将不会进入它的推导的理论背景。这超出了本文的范围。我更关心对 GBM 需要什么参数、它的组件是什么以及它如何创建预测的高层次理解。我将尝试一种自下而上的方法，从组件开始构建 GBM 的逻辑。我们在这里开发的仿真模型是一个离散时间模型。因此，这里讨论的所有数学都是连续随机过程的几何布朗运动的离散时间类比。在本文的最后，我们将学习如何使用 GBM 创建模拟，您将拥有完整的代码。
内容
- GBM 做什么
- GBM 模拟需要什么输入参数
- GBM 的组成部分:漂移和扩散
- 建立离散时间 GBM 模型
- 做预测
1。GBM 做什么
在解释相关概念时，我在整篇文章中使用 E.ON 的股票价格作为例子。E.ON 是一家位于德国的电力公司，也是欧洲最大的电力公司之一。我通过 Quandl 的 Python 包从 Xetra Exchange 检索它的股票价格(以欧元计)。这里有一个可以显示股票价格的链接:【Investing.com 链接
为了能够使用 Quandl，你需要注册并从其网站获得授权令牌，还需要安装***【Quandl】***Python 包。假设您已经完成了这些步骤，您可以使用下面的代码来提取股票价格数据。
几何布朗运动被广泛用于金融中的股票价格建模，人们选择它是有原因的。在下面的线图中，x 轴表示 2019 年 1 月 1 日至 2019 年 7 月 31 日之间的天数，y 轴表示以欧元为单位的股价。我希望你只关注图中主要的、持续时间较长的趋势，忽略小的波动。你会意识到股票价格是波动的。它从 1 月到 4 月增加，然后到 5 月中旬减少，接着是另一个增加系列，直到 6 月中旬，最后到 7 月底减少。现在，我要你关注短期波动。我们来调查一下七月。你可以说大势下行，但股价不只是平稳下跌。在它到七月底的旅程中，它不断上升和下降，这表明在这个运动中没有一个可解释的模式。
E.ON stock prices in the 1st half of 2019
因此，我们已经确定了股票价格运动的两个组成部分。七月，
- E.ON 的股价往往会下跌，而
- 每天都有随机的冲击导致股价下跌，形成不规则的线图。
因此，我们需要一个合适的模型来考虑股票价格的两种类型的运动。这就是几何布朗运动发挥作用的地方。GBM 有两个组件来完成这项工作。一个组件包含长期趋势，而另一个组件应用随机冲击。我们将在后面的章节中讨论这些。
2.GBM 模拟需要什么输入参数
以下是我们的 GBM 模拟模型将采用的输入参数。我们将逐一详细讨论它们。为了这篇文章，我将使用 E . ON 7 月份的股票价格来预测 8 月份的情况。另外，你应该注意到，当我用日期解释事情时，我说的是交易日。我们假设一周只由工作日组成，这意味着周五之后是周一。因此，我希望大家在阅读输入参数时记住这些内容。
我们检索从 开始日期 到 结束日期 之间的历史股票价格。然后使用我们的 GBM 模型，我们将得到我们的模拟直到 pred_end_date 。我们可以玩这些变量，创建不同的设置。
Date Close
0 2019-07-01 9.612
1 2019-07-02 9.761
2 2019-07-03 9.856
3 2019-07-04 9.800
4 2019-07-05 9.664Date Close
18 2019-07-25 9.517
19 2019-07-26 9.460
20 2019-07-29 9.488
21 2019-07-30 9.226
22 2019-07-31 9.054
注意，只检索交易日的股票价格，从上面的数据可以看出这一点。
1。所以
这是初始股价。远期股票价格模拟的根源就是这个初始股票价格值。我们将这个初始值乘以一些表达式，得到八月份每个交易日的预测值。在我们这里， 所以 就是 2019 年 7 月 31 日的收盘价。
9.054
2。dt
这是我们模型中的时间增量。它意味着我们假设的时间单位。请记住，我们有每个交易日股票的收盘价，那么我们模型中的时间增量是 1 天。请注意 dt 遵循历史数据中的时间步长。如果股票价格数据每 7.2 小时发布一次，那么 dt 将是 0.3，因为它对应于 7.2 小时(= 0.3 天)。
1
3。T
T 表示我们预测时间范围的长度。在指定 T 和 dt 的值时，我们应该小心并保持一致。例如，在我们的例子中， T 应该是 22 天，因为我们希望预测 8 月份的 22 个交易日，当给 dt 赋值时，根据我们对 T 的声明，我们应该记住 dt 必须用天数来表示。这两个参数的时间单位必须相同。现在，看一个不同的例子，假设我们的数据中每个交易日有两个股票价格值，我们知道我们将对八月份的 22 个交易日进行预测。在这个设置下，我们的时间增量将是 0.5 天，我们将 0.5 分配给 dt 。因此，在应用 GBM 后，我们将有 44 个预测，8 月的每个交易日有两个预测。
还记得上一部分中的例子吗？我们假设股票价格数据每 7.2 小时发布一次。如果我们希望我们的预测时间范围是 72 小时( T = 72)，那么合适的 dt 本身就是 7.2 小时，而不需要将其转换为 0.3 天。所以，最后，我们会有 72 / 7.2 = 10 个预测，每个预测之间间隔 7.2 小时。
如果你在构建和使用 GBM 模型时遵循这一思路，那么在不同环境下对不同股票使用你的模型就会变得容易得多。这让许多人感到困惑，这就是为什么我在这里试图将它标准化:)我们试图达到的最终目的是正确计算*(在下一部分解释)。它是在我们的预测时间范围内的时间点的数量，在时间增量大小方面应该与我们的历史数据一致。我知道这是一件简单的事情，但建立一条推理路线总是一个好主意，以防止将来我们代码的不同应用中的潜在混淆。*
我们可以使用本节开始时声明的变量**【pred _ end _ date】**来推断八月份的交易日数。使用下面的代码，我们可以通过计算( end_date + 1 天)和 pred_end_date 之间的工作日，提取我们的模型将预测股票价格的交易日数。在我们的案例中，我们需要的是 2019 年 8 月 1 日至 2019 年 8 月 31 日之间的交易天数。
*22*
**4。n
该参数在 dt 和 T 赋值后自动出现。它是预测时间范围内的时间点数量。在我们的例子中，我们的时间增量是 1 天，我们将得到 22 个交易日的预测。这意味着，我们有 22 个不同的时间点(天),最终我们将有 22 个预测。
*22.0*
5。t
这是一个数组，我们在其中显示了模型中的时间进程。它就像一个计时器，我们通过计算经过的时间点来测量时间。建筑数组 t 由 N 计算而来。记住我们有 22 个时间点(天)来做预测。这就是为什么数组 t 从 1 开始，一直到 22。当我们使用数组 t 中的一个元素时，它的意思是，模型中过去了那么多时间点。随着我们的预测距离 的日期越来越远，因此 (起始日)，随着时间点的推移，它们会因这个数组 t 而波动得更大。在我们的例子中， t 将是我们模型下面的数组。
*[ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22]*
其他一些例子:
一、 如果 dt 是 0.5 天(每天两个股票价格) T 是 22 天，那么 t :
*[1 2 3 4 5 6 7 8 9 10 ... 44]*
二世。 如果 dt 是 7.2 小时而 T 是 72 小时，那么 t :
*[1 2 3 4 5 6 7 8 9 10]*
我们之所以可以使用数组 t 作为一个正整数序列，直到 N 的原因是，我们在时间尺度上以同样的方式对待历史股价数据和我们的预测。我的意思是，历史数据中的 1 天和我们预测的 1 天是一样的。在网上的一些其他文章中，你可能会遇到这样的情况，数组 t 被缩小，上升到 1。在这些情况下，预测时间范围的 1 天变化的意义和我们的输入参数 mu 、 sigma、 和数组必须被缩放以说明时间-尺度差异。在本文中，我们对历史数据和预测使用相同的时间尺度，以使 GBM 模型更具解释性。
6。穆
这是所选历史日期范围内股票价格的平均回报。在我们的例子中，我们选择 7 月的交易日作为历史数据，我们应该使用 7 月的股票价格计算。在计算 mu 之前，我们要计算每个交易日的收益。计算如下。**
Return calculation for time point k
****[0.015501456512692377, 0.00973260936379476, -0.005681818181818097, -0.013877551020408265, -0.0032077814569537026, 0.00384096335513348, -0.013340227507755901, -0.005554973273241884, -0.008326306913996506, -0.0004251248804337693, -0.0005316321105793732, -0.0013829787234042447, 0.0030893789283050936, 0.0036108751062020944, 0.0010582010582012237, 0.0052854122621563355, 0.002944269190326022, -0.0022017194380374075, -0.0059892823368707165, 0.0029598308668074727, -0.0276138279932545, -0.018642965532191694]****
Mean return of the historical stock prices
在计算了七月份所有日子的收益之后，我们需要使用这些收益计算算术平均值，以获得 mu 。然后我们将使用进行我们的 漂移 分量的计算。它会对股价的长期走势产生影响。此外，正如您在下面看到的， mu 为~ -0.3 %，这表明 7 月份的平均回报率为负，我们应该在计算 8 月份的远期预测时考虑到这一点。**
**-0.0026705996693330543**
需要注意的一点是，我们不需要用一个数字乘以平均回报率来进行调整。我们忠于历史数据的时间增量，并以历史数据存在的方式创建预测序列。我们已经在上面讨论了这个问题，但是为了更详细地说明，假设我们每 7.2 小时检索一次股票价格数据，并且我们有过去 72 小时的数据。此外，我们还想预测接下来的 36 个小时。当我们从历史数据中计算平均回报时，该平均回报占 7.2 小时的相关时间步长(两个数据点之间有 7.2 小时)。在进行预测时，遵循我们对历史和未来时间尺度的一致性规则，我们将 7.2 小时分配给 dt。 这导致 5 个时间点(36/7.2)，它们之间有 7.2 小时，我们的数组 t 变成了【1 2 3 4 5】。这意味着我们可以在我们的模型中使用 mean return 作为 mu 而无需对其进行调整。**
7。西格玛
这是 7 月份股票价格收益的标准差。西格玛很重要，因为，它将有助于把随机冲击纳入我们的预测。还记得上面 E .关于股票价格的线图吗，股票价格一天比一天持续上涨和下跌。还有，这个小运动的幅度和方向似乎是随机的。在这里，∑将帮助我们确定运动的幅度。 sigma ，本身并没有给模型增加我们需要的随机性。当选取随机值时，我们将使用标准的正态随机变量。西格玛将通过缩放随机冲击的幅度做出贡献，以便根据股票价格的历史波动性出现小波动。我们不希望任何不相关的随机值来自标准正态分布。下面，你可以看到 sigma 是如何计算的。
The volatility of the historical stock prices
**0.009362443851856311**
记得我们说过，我们不会将 mu 与任何数字相乘进行调整，因为我们在进行预测时会忠于历史数据中的时间步长。同样的逻辑也适用于此，我们也不需要调整。**
8。b**
这个数组是我们在模型中添加随机性的数组。这里重要的是 scen_size 变量。由于这是一个模拟模型，我们希望在预测时间范围内有许多股票价格的假设路径。将 2 赋给 scen_size 意味着，最终我们会有 2 个不同的股价序列。每个场景的逻辑都是一样的。因此，我将讨论我们是否只创建一个预测系列。还记得在 t 的讨论中，我们声明了一个计算经过时间点的时间推进数组。这里，数组 b ，对于每个相应的预测时间点，存储来自标准正态分布的随机数。这些随机数会给模型增加随机冲击，这一点我们从文章开始就一直在讨论。这就是为什么我们将使用 numpy.random.normal() 从标准正态分布中生成随机值。同样，我们不将这个随机值与任何数字相乘进行调整，遵循与**和相同的推理。****
Creating random shocks
在上面的例子中，我们假设我们的历史股票价格数据是每日的，这就是为什么我们将 1 分配给 dt 。同样，我们想要 5 天的预测，这意味着 T 的值是 5。接下来， N 被计算为 5 个时间点(5 个 1 天周期)。因此，对于每个时间点，我们使用 numpy.random.normal() 方法生成一个来自标准正态分布(mean=0，std=1)的随机数。这些随机值最终组成了我们的数组 b 。
在我们的例子中， T 应该是 22 天。下面的代码中，我们为两种不同的场景创建了两个长度为 22 的随机值数组。
**{'1': array([ 1.37139893, -1.60468185, 0.84629302, 0.84599139, 0.1392499 , -0.23658454, 0.66158754, 1.44276449, 0.25582332, -0.00432083, 1.39636253, 0.21584903, 0.15232219, 1.41902316, 0.15672926, -0.76728718, 0.26303217, 1.61346848, 0.16319213, 1.08515534, -1.43681955, -1.89643601]),'2': array([-1.19967538, -0.66905682, -1.31963702, -0.04374916, 1.36622513, 0.39133558, -0.05713623, -1.16944076, 0.51635002, -0.53486539, 1.864583 , -1.2722681 , 0.06091857, -0.86510476, 0.1338046 , 0.76766205, -0.00538846, 0.92632565, -0.13265282, -0.12242566, 0.4031829 , 0.74874924])}**
9。w
【W】是布朗路径，它决定了股票价格从开始时间点( So )到另一个时间点 t )如何波动。你要区分 b 和 W 。在下一节中，它们之间的区别将会更加清楚，但我还是想在这里简单提一下。 b 是在预测下一个时间点的股价时，被施加到一个时间点股价上的随机冲击。所以，假设，在时间点 3，股价是 S_3。当预测时间点 4 时，b【4】作为随机冲击应用于 S3。 W 则相反，是路径。这意味着它包括了自预测时间范围开始以来所有随机冲击的影响。是随机性并入 所以 (初始股价)直到我们关注的特定时间点的总效应。假设我们预测时间点 4。我们需要将所有随机冲击应用到时间点 4 到 ，因此 。所以我们用索引小于等于 4 的数组 b 元素的累积和【W】【4】代替b**【4】。下面是上一部分中我们讨论数组 b 的例子的延续。****
Creating the Brownian paths
**{'1': array([ 1.37139893, -0.23328292, 0.6130101 , 1.45900149, 1.59825139, 1.36166685, 2.02325439, 3.46601888, 3.72184221, 3.71752138, 5.11388391, 5.32973294, 5.48205514, 6.9010783 , 7.05780756, 6.29052038, 6.55355255, 8.16702103, 8.33021316, 9.4153685 , 7.97854894, 6.08211294]),'2': array([-1.19967538, -1.8687322 , -3.18836922, -3.23211838, -1.86589325, -1.47455767, -1.5316939 , -2.70113466, -2.18478464, -2.71965004, -0.85506704, -2.12733514, -2.06641657, -2.93152133, -2.79771673, -2.03005468, -2.03544314, -1.10911749, -1.24177031, -1.36419597, -0.96101307, -0.21226383])}**
以上是我们案例的代码。继上一部分的数组 b 计算之后，我们根据上面的 W (k)表达式取累加和，创建数组 W 。
我们对 GBM 模型输入参数的讨论到此结束。我们通过例子理解了它们背后的原因，在接下来的部分中，我们将从它的组件中构建 GBM 模型。这将使本节中的输入参数对您更有意义。
3.GBM 的组成部分:漂移和扩散
记得在第一节，我们已经确定了几何布朗运动的两个组成部分。一个是股票价格的长期趋势，另一个是短期的随机波动。现在，我们将为它们命名。我们称长期趋势为漂移，称短期波动为扩散。在本节中，我希望你假设我们在一个特定的时间点(k-1)，我们将预测时间点(k)。然而，我们在上一节中谈到的相同假设仍然有效。提醒一下:
假设 1:(k-1)和(k)之间的时间段长度，也就是 dt，符合历史数据频率。
假设 2: 我们模拟中的时间是通过计算时间段来进行的。
这些假设使我们忽略了将 mu ，【sigma】， z (k)与一个时间长度值相乘进行调整，我们可以组成我们的数组 t 作为一个时间点计数器【1 2 3…】。****
我们继续吧。我们在一个时间点(k-1)。在这个时间点的股票价格需要服从股票价格的长期趋势，同时当我们使用它来创建时间点(k)的预测时，暴露于随机冲击。因此，我们需要将 GBM 的两个组成部分应用到股票价格中。
让我们用数学术语来看看这些组件是什么:
**1。漂移
Drift for time point k
漂移 反映股票价格的长期趋势。大家可以回忆一下上一节，我们已经知道了 mu 和 sigma 。它们是从 7 月份的历史股价数据中计算出来的常量值。因此，我们的 漂移 是一个常数值，如果我们假设不存在随机冲击，我们可以将其应用于时间点(k-1)的股价，如上面的表达式所示。如果我们继续应用 漂移 而没有任何随机冲击，那么如果漂移为正，股价就会平稳上升；如果 漂移 漂移为负，股价就会下跌。你应该注意到股票价格永远不会变成 0，因为 漂移 是通过指数函数应用的。这就是 漂移 的故事。现在，让我们用**扩散移动扩散*。*****
**2。扩散
Diffusion for time point k
扩散反映短期波动。从第 2 节可以回忆起，数组 b 存储了我们需要的随机冲击信息，它从标准正态随机变量 z 中检索随机冲击信息。在计算扩散分量时，我们将随机值 z (k)乘以∑*。在这里，您可以注意到随机性是如何融入 GBM 模型的。另外， 扩散 组件使得创建不同的股票价格预测场景成为可能。 漂移 分量不变。我们不能从中创建不同的场景，但 扩散 组件帮助我们创建尽可能多的场景，因为它涉及维纳过程(它创建独立、平稳和正态分布的随机冲击)。*****
因此，当考虑到时间点(k-1)的股票价格对时间点(k)进行预测时，我们需要将这两个因素的影响相加。总的组合效应为我们提供了对时间点(k)的预测。如果只有 漂移 成分，我们总是知道股票价格在下一个时间点会变成什么，但是 扩散 成分通过引入随机冲击来破坏这种平滑性。因此，我们没有观察到股票价格的单调上升或下降趋势。下面是总组合效应如何应用于时间点(k-1)的股价。
Adding drift and diffusion for 1 time point ahead prediction
4.建立离散时间 GBM 模型
到目前为止，我们学习了为什么我们需要 GBM，离散时间 GBM 模型的参数是什么，以及如何提前 1 个时间点进行预测。现在，我们建立一般的闭合形式的方程的几何布朗运动调整为离散时间的背景。在上一节中，我们展示了如何通过将 S (k-1)乘以 exp( 漂移 + 扩散 ) 来预测 S (k)。作为此的延续，我们将尝试预测 S (k)使用 So (已知的初始股价)。
为此，我们首先将 S (k-1)替换为
**S(k-2)exp(漂移(k-1)+扩散 (k-1) )。
然后我们将S***(k-2)替换为*****
S(k-3)exp(漂移 (k-2) + 扩散 (k-2) ) 。
如果我们继续这样做，最后，我们将有 所以 乘以许多指数项。您可以在下面的步骤结束时看到我们获得的结果。
之后，我们扩展上面的等式，这次分别表示和 扩散 项。**
以上是倒数第二个方程式为 S (k)。这意味着，我们可以从 预测 S (k)所以 通过添加 k 个 漂移 和累积直到 k 的组合效应记住从上一节开始， 漂移 不会从一个时间点变化到另一个时间点。这就是为什么，要预测前方 k 个时间点，我们只需多次添加 漂移 k。这与 扩散 不同。在第 2 节讨论 W 时，我们了解到它是布朗路径。你可以回忆一下那部分。基于该讨论，为了预测 k 个时间点，我们应该考虑在路径上遇到的所有随机冲击的总效应。在建立等式的最后一步，我们通过替换下面的等式来修改上面的等式。**
然后，它给出了方程的最终形式，其中的参数我们都很熟悉，我希望我们在第 2 节中进行的所有讨论现在能帮助您理解我们为什么以这种方式创建这些参数。
我们知道所以我们知道 mu 我们知道 sigma 我们知道数组 t 我们知道数组 W 。我们可以一步轻松地为预测时间范围内的所有时间点创建预测。我们将在下一节中完成这项工作。
Geometric Brownian Motion equation
5.做预测
这是我们文章的最后一部分，也是最有趣的部分。我们现在可以创建预测。首先要做的是计算预测时间范围内所有时间点的 漂移 。你已经记住了**阵列t。我们只需将它与 漂移 相乘，就可以得到一个 漂移 的数组。该 漂移 数组包含预测时间范围内所有时间点的总 漂移 。****
关于 扩散 ，我们将有一个 扩散 数组用于每一个场景。记住，我们使用变量 scen_size 来控制我们想要多少场景。现在，对于这个问题，我们只有两种情况。你想涨多少就涨多少。这就是为什么你会看到 2 个不同的阵列，分别标有‘1’和‘2’。但是，对他们每个人来说，这是同样的逻辑。我们将布朗路径(数组 W )乘以 【西格玛】 。****
**drift:
[-0.00271443 -0.00542885 -0.00814328 -0.01085771 -0.01357214 -0.01628656 -0.01900099 -0.02171542 -0.02442985 -0.02714427 -0.0298587 -0.03257313 -0.03528756 -0.03800198 -0.04071641 -0.04343084 -0.04614526 -0.04885969 -0.05157412 -0.05428855 -0.05700297 -0.0597174 ]diffusion:
{'1': array([ 0.01283965, -0.0021841 , 0.00573927, 0.01365982, 0.01496354, 0.01274853, 0.01894261, 0.03245041, 0.03484554, 0.03480509, 0.04787845, 0.04989933, 0.05132543, 0.06461096, 0.06607833, 0.05889464, 0.06135727, 0.07646328, 0.07799115, 0.08815086, 0.07469872, 0.05694344]),'2': array([-0.01123189, -0.0174959 , -0.02985093, -0.03026053, -0.01746932, -0.01380546, -0.0143404 , -0.02528922, -0.02045492, -0.02546257, -0.00800552, -0.01991706, -0.01934671, -0.0274462 , -0.02619347, -0.01900627, -0.01905672, -0.01038405, -0.011626 , -0.01277221, -0.00899743, -0.00198731])}**
以下是预测。因为我们有两个场景，所以作为下面操作的结果，我们将有两个股票价格预测序列。第一行代码就是我们在上一节中推导出的最终 GBM 方程。第二行只是将 So 添加到预测序列中，因为它是起点，我们可能希望在我们的图中看到它。
**[[9.054 9.1461394 8.98533403 9.03226024 9.07940588 9.06660652 9.02202324 9.05347139 9.15171815 9.14879651 9.12362735 9.21862966 9.21223827 9.2003776 9.29815157 9.28656346 9.19509744 9.19278238 9.30740419 9.29636703 9.36583938 9.21564323 9.02891936][9.054 8.92860643 8.84880033 8.71645344 8.6892655 8.77726795 8.78560531 8.75710402 8.63826721 8.65659865 8.59000918 8.71758658 8.59101135 8.57261126 8.48040716 8.46802049 8.50598079 8.48249529 8.53318693 8.49949308 8.46674294 8.47572541 8.51221286]]**
作为最后要做的事情，我们现在创建我们的预测的情节。正如你在下面看到的，我们有两个股票价格预测场景。一个上升，另一个下降。我们问题中的 漂移 是负的，你可以从 漂移 数组中观察到。这意味着股票价格的长期趋势是下降的(记住我们只考虑 7 月份的历史数据。我是在这个基础上说的！).然而，尽管有一个负的 漂移 ，我们有一个向上的预测。这是由于随机冲击。记住，我们使用标准的正态随机变量来创建它们。在价格上涨的情况下，很明显，这些随机冲击超过了 漂移 效应，导致了这样一个系列。
Two different scenarios for August stock prices
结论
在本文中，我们学习了如何在离散时间环境中使用几何布朗运动来构建股票价格的模拟模型。下面是完整的代码。当您在注册后从 Quandl 获取授权令牌并安装所需的 Python 包时，您可以立即使用代码。在每次运行中，你会得到不同的股价情景。是因为我们使用 np.random.normal() 没有设置种子。每次我们运行模型，我们将有一个不同的数组 W ，这将导致不同的预测。
几何布朗运动帮助我们看到股票价格可能遵循的路径，并让我们为即将到来的事情做好准备。当然，预测准确的未来是不可能的，但是这些统计方法给了我们创造可靠的交易和对冲策略的机会。
新浪 _ 剧情娱乐与 ggplot
我最近被要求重新创建一个 ggplot 图表。
knitr::include_graphics("images/image002.jpg")
通过查看它，我能够确定一些虚拟数据，让自己开始…
*# dummy data*
small_group_size <- 20
large_group_size <- 2000
set.seed(123)
small_sample_1 <-
tibble(value1 = rnorm(small_group_size, 70, 20),
value2 = rnorm(small_group_size, 150, 20),
group = "Group A")
small_sample_2 <-
tibble(value1 = rnorm(small_group_size, 30, 20),
value2 = rnorm(small_group_size, 180, 20),
group = "Group B")
control <-
tibble(value1 = rnorm(large_group_size, 50, 20),
value2 = rnorm(large_group_size, 220, 20),
group = "Group C")
*# let's put it all together and add a gender column, as well*
data_collected <- bind_rows(
small_sample_1,
small_sample_2,
control) %>%
mutate(person_id_number = row_number()) %>%
mutate_at(c("value1", "value2"), round, 2) %>%
mutate(gender = sample(c("M", "F"), n(), replace = T)) %>%
mutate(group = factor(group, levels = c("Group A", "Group B", "Group C")))
*# and tidy things up (if you want... object nuking, everying but)*
rm(list = setdiff(ls(), "data_collected"))
我发现在 ggforce 包里也有类似的东西，叫做 sina_plot。相似但不完全…来看看开箱的新浪 _ 剧情吧
data_collected %>%
ggplot() +
aes(x = group, y = value1) +
ggforce::geom_sina()
我看到的图表中有一些有趣的特征，是直线 sina_plot 没有给我的。首先，y 轴值被分类，刻度是某种形式的对数。加上 x 轴抖动是一个不规则的。让我们看看我们能否实现这一点…
given_levels <- data_collected %>%
mutate(value_measured = value1,
level_group = (value1 %/% 5) * 5,
level_group = pmax(level_group, 0)) %>%
group_by(group, level_group) %>%
mutate(count_of_entries = n(),
is_group_count_even = count_of_entries %% 2 == 0)
现在，让我们用颜色来看看。
given_levels %>%
ggplot() +
aes(x = group, y = level_group) +
ggforce::geom_sina()
因此，看起来我必须确定数据准备中的 x 轴和 y 轴位置，并为不同的组使用一个面…并使用 geom_point 图表…稍后再担心颜色。
given_levels_and_width_values <- given_levels %>%
mutate(count_of_entries = n(),
width_of_group = count_of_entries %/% 2,
width_position = row_number() - width_of_group - 1) %>%
mutate(width_position = if_else(is_group_count_even & width_position == 0, width_of_group, width_position)) %>% *# this is tricky... Basically, if there's an even number, we want to relocate that x=0 point somehow...by putting it at the end... this is actually going to cause us problems later :-)*
mutate(width_position = case_when(is_group_count_even & width_position > 0 ~ width_position - 0.5,
is_group_count_even & width_position < 0 ~ width_position + 0.5,
TRUE ~ width_position))
让我们现在画出来
given_levels_and_width_values %>%
ggplot() +
aes(width_position, level_group) +
geom_point() +
facet_grid(. ~ group,
scales = "free_x")
太棒了。稍微整理一下…
distance <- 5
unique_grouping_levels <- unique(given_levels_and_width_values$group)
dummy_data <- tibble(column_to_group = rep(unique_grouping_levels, 2),
width_position = c(rep(-distance, length(unique_grouping_levels)),
rep(distance, length(unique_grouping_levels))),
level_group = rep(NA_integer_, length(unique_grouping_levels)*2))
median_lines <- given_levels_and_width_values %>%
group_by(group) %>%
summarise(level_group = median(value_measured, na.rm = T),
to = max(width_of_group, na.rm = T) + distance/2,
from = -to) %>%
gather(key, width_position, -1, -2) %>%
select(-key)
break_steps <- seq(min(given_levels_and_width_values$level_group, na.rm = T), max(given_levels_and_width_values$level_group, na.rm = T), 25)
nearly_finished_plot <- given_levels_and_width_values %>%
ggplot() +
aes(x = width_position, y = level_group) +
geom_point(shape = 1) +
geom_blank(data = dummy_data) +
labs(x = NULL,
y = NULL) +
scale_y_continuous(breaks = break_steps) +
facet_grid(. ~ group,
scales = "free_x",
**switch** = "both") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
strip.background = element_rect(fill = "white"),
strip.placement = "outside",
panel.spacing.x = unit(0, "lines")) +
geom_line(data = median_lines, colour = "red", size = 1)
nearly_finished_plot
我们可能还想看到由另一列着色的点…
only_colour_points <- given_levels_and_width_values %>%
filter(!is.na(gender)) %>%
arrange(gender) %>%
mutate(count_of_entries = n(),
width_of_group = count_of_entries %/% 2,
width_position = row_number() - width_of_group - 1) %>%
mutate(width_position = if_else(is_group_count_even & width_position == 0, width_of_group, width_position)) %>% *# this is tricky... Basically, if there's an even number, we want to relocate that x=0 point somehow...by putting it at the end... this is actually going to cause us problems later :-)*
mutate(width_position = case_when(is_group_count_even & width_position > 0 ~ width_position - 0.5,
is_group_count_even & width_position < 0 ~ width_position + 0.5,
TRUE ~ width_position))
nearly_finished_plot +
geom_point(data = only_colour_points, aes(colour = gender)) *# these colour points are actually being overlayed over the original points, and loaded in from the left... this is beneficial because it allows you to see proportions (roughly) of each group (and even NA values, if there are any)*
好吧，这一切都很好。我们已经制定了数据准备步骤和所有的小 ggplot 添加…从这里，让我们把这些步骤变成函数。(1)分配宽度组(箱)(2)分配宽度位置(3)如果适用，通过颜色变量分配宽度位置(4)绘制！
现在，让我们尽可能的动态化…这样 y 轴和 faceting(分组)变量可以在函数调用中指定。
assign_width_groups <- **function**(data, column_to_group, column_to_bin, min_value_to_plot, bin_size) {
column_to_bin <- rlang::enquo(column_to_bin)
column_to_group <- rlang::enquo(column_to_group)
data %>%
mutate(value_measured = !! column_to_bin,
column_to_group = !! column_to_group,
level_group = ((!! column_to_bin) %/% bin_size) * bin_size,
level_group = pmax(level_group, min_value_to_plot)) %>%
group_by(column_to_group, level_group) %>%
mutate(count_of_entries = n(),
is_group_count_even = count_of_entries %% 2 == 0)
}assign_width_positions <- **function**(data) {
data %>%
mutate(count_of_entries = n(),
width_of_group = count_of_entries %/% 2,
width_position = row_number() - width_of_group - 1) %>%
mutate(width_position = if_else(is_group_count_even & width_position == 0, width_of_group, width_position)) %>%
mutate(width_position = case_when(is_group_count_even & width_position > 0 ~ width_position - 0.5,
is_group_count_even & width_position < 0 ~ width_position + 0.5,
TRUE ~ width_position))
}plotting_function <- **function**(data, distance = 5, colour_selection = NULL) {
colour_selection <- rlang::enquo(colour_selection)
unique_grouping_levels <- unique(data$column_to_group)
dummy_data <- tibble(column_to_group = rep(unique_grouping_levels, 2),
width_position = c(rep(-distance, length(unique_grouping_levels)),
rep(distance, length(unique_grouping_levels))),
level_group = rep(NA_integer_, length(unique_grouping_levels)*2))
median_lines <- data %>%
group_by(column_to_group) %>%
summarise(level_group = median(value_measured, na.rm = T),
to = max(width_of_group, na.rm = T) + distance/2,
from = -to) %>%
gather(key, width_position, -1, -2) %>%
select(-key)
break_steps <- seq(min(data$level_group, na.rm = T), max(data$level_group, na.rm = T), 25)
plot <- data %>%
ggplot() +
aes(x = width_position, y = level_group) +
geom_point(shape = 1) +
geom_blank(data = dummy_data) +
labs(x = NULL,
y = NULL) +
scale_y_continuous(breaks = break_steps) +
facet_grid(. ~ column_to_group,
scales = "free_x",
**switch** = "both") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
strip.background = element_rect(fill = "white"),
strip.placement = "outside",
panel.spacing.x = unit(0, "lines")) +
geom_line(data = median_lines, colour = "red", size = 1)
*# we also need to add a final step for when a colour component is selected, as well.*
**if** (!is.null(colour_selection)) {
only_colour_points <- data %>% *# this is the data after both steps of the data-prep...*
mutate(colour_column = !! colour_selection) %>%
filter(!is.na(colour_column)) %>%
arrange(colour_column) %>%
assign_width_positions()
plot <- plot +
geom_point(data = only_colour_points, aes(colour = colour_column))
}
plot
}
现在，只需几行代码就可以将它们组合在一起
data_collected %>%
assign_width_groups(gender, value1, 0, 5) %>%
assign_width_positions() %>%
plotting_function(colour_selection = group)
中线也很棘手。首先，它落在彩色点的后面(这很容易修复)，其次，它是从实际值而不是入库值计算出来的(这就是为什么它不完全在某个显示点上或之间)。
还有一个小问题，但我认为它不会被注意到，直到我能费心去解决它(这不是对数刻度…我认为没有它看起来更好)。
不管怎样，所有这些功能都准备好了，我们现在可以把它们放入一个漂亮的小 flexdashboard 中…
结论
天啊，这花了很长时间才解决。与 Tableau 相比，在 ggplot / R 中定制的(和可定制的)图表会花费很长时间…但是现在我知道怎么做了…小心世界:-)实现下一次将是一个 sinch。
新加坡统一价格预测器
用机器学习预测新加坡房价
Photo by Peter Nguyen on Unsplash
新加坡连续第五年被经济学人信息部选为最昂贵的国家之一。作为世界上最昂贵的城市之一，住宿占了很大一部分费用。
为了给新加坡普通民众提供负担得起的住房，新加坡政府一直在新加坡提供一些公共住房选择，更普遍的说法是 HDB 公寓。
对于这个附带项目，我已经采用了 HDB 转售价格数据，并训练了一个机器学习模型来预测新加坡的统一转售价格。
根据智能国家倡议，新加坡政府在网上发布了大量新加坡数据，公众可以访问。此外，有一些很好的 API 可供开发人员集成到他们的应用程序中。
概观
这篇文章被分成几个部分。以下是我所涉及的部分的概述:
- 将数据可视化
- 特征创建
- 训练/测试分割
- 实验
- 创建机器学习模型
- 全体
- 韵律学
- 结果
- 结论
事不宜迟，让我们深入了解一下如何自己创建一个统一价格预测器！
将数据可视化
在开始一个机器学习项目之前，我们执行一些标准的可视化程序来理解我们的数据是很重要的。只有了解我们的数据，我们才能知道清理和预处理数据的正确方法。
以下是该数据中可用的列:
'month', 'town', 'flat_type', 'block', 'street_name', 'storey_range', 'floor_area_sqm', 'flat_model', 'lease_commence_date', 'resale_price'
转售价格将是这个项目的目标标签，其余的列将形成我们的特征矩阵。让我们继续画一些图表。
Distribution Plot for Floor Area against Resale Price
我们可以从图表中推断出，建筑面积越高，公寓的转卖价格就越高。通过绘制分布图，我们可以发现是否有任何异常值，如果它们看起来太不寻常，就将其删除。
除此之外，我们还可以使用箱线图来了解我们的数据分布。下面是一个简单的箱线图，它有助于向我们展示给定特定特征的目标标签的分布。
Image from Flowingdata
在我们的项目中，boxplot 可以帮助我们了解不同公寓类型中的公寓价格分布。
另一个有助于我们理解特征的重要可视化工具是关联矩阵:
Correlation matrix of different flat’s features
那么，从这个矩阵中能得到什么呢？建筑面积和转售价格之间的颜色是其他组合中最浅的。这意味着建筑面积和转售价格高度正相关——更高的建筑面积，更高的转售价格。
相比之下，厕所数量和转售价格的相关性非常低，这意味着厕所数量在预测公寓的转售价格方面可能不是很有用。
特征创建
Photo by Kelly Sikkema on Unsplash
该数据集中可用的特征相当小。但是，没有什么能阻止我们创造更多！
拥有更多功能肯定会帮助机器学习模型学习数据集的模式。
一些我认为有用的功能:
- 方向。这套公寓所在的城镇在这个数据集中给出了。我们如何改进这一功能？我创建了一个新特征，将城镇映射到下面列出方向/区域之一:
"NORTH_EAST","EAST""WEST","NORTH","CENTRAL"
2.房间号&卫生间号。根据给定的公寓类型，我从 HDB 网站上查看，找出房间和厕所的数量，并创建了两个新功能。
3.层高。数据中提供的楼层范围要素是字符串格式:4 到 6，7 到 9，我认为这是次优的，因为如果您的要素是序数，您应该将要素转换为数字。
原因在于，楼层 7 低于楼层 9，将特征建模为数字将有助于模型了解这种关系。如果我们将特征保留为字符串而不进行任何处理，模型将会错过这些重要的信息。
所以我创建了一个随机生成器，它将在数据集中给定的最低楼层和最高楼层之间生成一个数字。请看下面的例子，以便更好地理解:
# Example, :
4 To 6 - The generator output: 4 or 5 or 6
7 TO 9 - The generator output: 7 or 8 or 9
4.**销售年度&月。**数据集中最初给出的月份是这样的格式:2018–4，因此通过将它们分解为“年”和“月”将允许模型选取序数关系。
新创建的功能在提高模型性能方面的效果如何？我将运行多个实验，并与基本模型进行一些比较。我将在结果的后面部分介绍这一点。
训练/测试分割
在将模型投入生产之前，我们需要将数据分为训练集和测试集，以验证模型的准确性。创建我们的训练测试集时需要注意的一些重要事项:
- 训练和测试集的分布应该类似。训练和测试集的一个不好的例子是，训练集中的单元都来自东部地区，而测试集中的单元都来自西部地区。显然，在这种条件下训练的模型不会非常有效。
- 测试集的分布和现实生活中的数据应该是相似的。测试集是一个安全的评估，让您在投入生产之前了解系统的性能。如果你的测试集没有捕捉到真实生活数据的分布，那么你从测试集中得到的结果不会给你多少信心。您的系统很有可能在生产中表现不佳。
为了产生一个好的训练和测试数据集，我们可以使用**分层抽样。**在我们进行分层抽样之前，我们需要确定一个重要的关键特征，并确保您的训练和测试集分布基于该特征是相似的。
例如，在我们的项目中，我选择位置作为分层抽样的关键特征。就位置而言，列车和测试设备中的单元具有相似的分布是一个好主意。也就是说，如果您的列车组在中部有 50%的平坦度，在西部有 30%的平坦度，在东部有 20%的平坦度，那么您应该在测试组中保持这个比例。
实验
数据科学的很大一部分是关于进行实验的。创建新功能、更改模型、调整超参数等等。所有这些活动的有效性需要通过运行大量的实验来衡量。
在这一部分，我将介绍一些技巧和工具，它们可以帮助你更有效地进行实验。
MLflow
Mlflow 是一个帮助您简化机器学习开发的平台。MLflow 下有多个子项目，例如:
- MLflow Tracking —帮助您跟踪实验的参数、结果和指标，并通过交互式用户界面进行比较
- MLflow 项目s——帮助您通过 Docker 和 Conda 将代码打包成可重复运行的代码
- MLflow 模型 —帮助您打包模型，以便您可以轻松地与他人共享
在这个项目中，我使用 MLflow 来跟踪我的实验，并在训练后记录模型，以便我可以很容易地与他人分享我的结果。下面的代码片段显示了 MLflow 提供的交互式 UI:
MLflow UI
对于您在实验中所做的任何更改或额外处理，留下一些描述是一个很好的做法。这将有助于您稍后分析结果。
使用较小的样本量进行实验
我从一位了不起的教练那里学到了这个建议——杰瑞米·霍华德。他是 fastai 的创始人，该公司的使命是让人工智能教育对任何人免费开放。
这是一门非常棒的关于机器学习的免费在线课程。如果你有兴趣亲自动手，并开始你的机器学习之旅，一定要去看看他们的课程。
那么这实际上意味着什么呢？在前面的部分中，我们已经解释了训练测试分割过程，这样我们就可以为我们的模型获得一些反馈。现在，我们可以更进一步，从我们的训练集中创建一个更小的样本，用于实验。
我们为什么要这么做？当我们设计模型或添加新功能时，我们希望从我们正在进行的所有活动中获得快速反馈。等待 5-10 分钟去做一件事是不会有成效的。因此，拥有更小的样本量，你可以更快地运行你的实验，获得反馈并做出改变。
样本量应该有多小？经验法则——它应该很小，但足够有代表性，可以帮助你在不到 20 秒的时间内完成实验，并获得对结果的粗略估计。
创建我们的机器学习模型
Sklearn 提供了许多方便的库，帮助我们轻松地创建模型。您可以创建多个模型，并在一个实验中运行所有模型，请参考下面的代码:
elasticnet = make_pipeline(RobustScaler(), ElasticNet(alpha =0.00001, random_state=1))svm = make_pipeline(RobustScaler(),SVR(C=140))lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))GBoost = GradientBoostingRegressor(loss='huber',min_samples_split=10,n_estimators=200)experiment_description = [
'Add in sale year and sales month feature',
'Add more data from 2015'
]run_experiment([elasticnet,svm,lasso,GBoost],experiment_description)
代码可能看起来非常简单，但是它需要大量的知识来理解这些模型背后的理论。例如，你需要知道模型背后的数学理论是什么，超参数是什么，何时何地我们应该使用这个模型。
每一个模型都值得用一整篇文章来彻底解释其背后的理论和机制。在网上查找资源，那里有大量有用的资料。
全体
Photo by Debby Hudson on Unsplash
这是机器学习中非常有用的技术，可以提高你的整体模型的准确性。它所做的是创建一堆不同的分类器(例如:random forest、elasticnet、xgboost)，并集成这些分类器的结果以改善结果。
我们可以使用的一种简单而常见的集成技术是简单地对每个分类器的预测进行平均。
你可能想知道它有什么帮助？背后的基本原理是，不同的分类器从不同的角度学习数据集的模式。通过收集和平均结果，我们可以看到一些性能提升。
堆垛
我在这里使用的整体技术是堆叠。与平均模型相比，这是一种稍微先进的技术。然而，它在引擎盖下做的事情非常简单！下面是堆叠的简单架构:
Image from Wikimedia commons
我们有多个分类器对数据进行学习，并输出许多不同的预测。我们没有对预测进行平均，而是将其发送给元学习器，元学习器是另一个机器学习模型，可以给我们最终的预测。
元学习者决定不同分类器所做预测的权重。一些分类器可能比其他分类器表现更好，因此元学习器将优化并找到最佳组合。
虽然整个过程听起来令人生畏，但这实际上只需要很少的代码，因为繁重的工作是由我们出色的 sklearn 库完成的！
韵律学
是时候看看结果了！但是在分析结果之前，我们首先需要了解指标是什么。
我们这里使用的指标是均方根对数误差(RMSLE)。我相信你们大多数人都听过和学习过均方根误差(RMSE)。这两个指标的区别在于，在计算均方根之前，我们对目标值和预测值应用对数。参见下面的等式:
Equation of root mean square logarithmic error
为什么我们用 RMSLE 而不是 RMSE？通过做对数，我们只查看实际值和预测值之间的相对差异/百分比差异。
Comparison between RMSE and RMSLE
如果目标值和预测值之间的差异在数值上很大，RMSE 会给予很大的惩罚。但 RMSLE 正在研究百分比差异，并对两种情况给出相同的损失分数。
在我们的项目中，不同公寓的价格在不同的水平上变化，我们希望所有的误差都以百分比为基础来处理，这是使用 RMSLE 的主要原因。
结果
Comparison of results between the effects of feature creation
损失分数越低，目标值和预测值之间的差异越小，因此意味着更好的结果。通过以上比较，我们可以知道，在所有的特征创建中，创建平面位置到区域的映射是最有帮助的。
通过组合我们创建的所有额外功能，总体改进是非常显著的。您可以比较第一列和最后一列的结果，以观察损失分数的改善情况。
我们用于执行所有这些实验的数据大小只是数据的一个非常小的子集。我们可以运行完整的实验(使用所有的训练数据)并比较“最基本的模型”——模型 A 和“具有所有额外特征的模型”——模型 B ，以验证我们的特征创建实际上是有帮助的。请参见下面的结果:
Comparison between most basic models and models with all extra features
在我们使用的所有型号中，型号 B 的表现要好得多。所以功能创建肯定很有帮助！
我想强调一下这两个模型中性能最好的模型，即堆叠分类器。模型 B 的堆积分类器的损失分数比模型 A 低 15%,这是一个惊人的改进。
结论
简而言之，我已经强调了从零开始构建统一价格预测器的整个过程。你可以从这个 GitHub 库下载源代码。
除了机器学习代码，我还创建了一个 web 应用程序来执行实时预测。
web 应用程序代码在同一个存储库中，如果您觉得对您的用例有用，可以参考它们。
使用 GraphX 跨设备为匿名用户匹配单个用户标识
Public domain image, Pxhere
有可能你正在阅读这篇文章，因为你有一个数字营销的用例。 Gartner 称之为跨设备识别(XDID) 数字营销人员需要关注的新兴技术。此外，当前对 cookie 管理的监管变化进一步增加了概率模型和数据科学在单个用户匹配中的重要性。
但是数字营销并不是唯一可以从单一用户 ID 中获益的用例。其他有价值的用例围绕着推荐系统和个性化，并且与行业无关。
识别用户的不同选项有:
账户登录:一种高度确定性的用户识别方式。这种形式的身份识别的质量是以注册和登录用户数量通常较少为代价的。因此，并非所有返回的用户都必须被识别为单个用户。对于拥有多个帐户的用户也是如此。
Cookie 标识:通常是一种高质量的标识形式，但不能跨设备工作。cookie 许可管理和 cookie 寿命的法规变化也带来了越来越多的问题。
设备指纹:一种识别用户身份的概率方法，但在设备指纹的使用方面受到越来越多的监管和技术限制。此方法也不能跨设备工作。
IP 匹配:识别单个用户的较低质量的概率选项。然而，IP 是跨设备可用的。与 3G/4G 网络或公共 Wifi 一样，动态共享 IP 也存在问题。无论如何，IP 是最可行的选择。(尽管根据 GDPR，知识产权是“个人身份信息”。)
Mozilla Lightbeam, public domain image, Wikimedia
在 Zoopla，消费者可以浏览我们的房源列表并发送线索。但是，我们在这个过程中对用户注册和登录只有相对较弱的激励。因此，许多 Zoopla 用户实际上是匿名的，并且通常使用多种设备。为了创建一个强大的个性化路线图，我们必须找到一种方法来创建单个用户 id，它可以跨设备工作，也可以处理匿名流量。在 Zoopla 每天超过 200 万次会话的规模上，这是一个有趣的数据科学问题，也是一个困难的数据工程挑战。
这篇博文展示了我们如何使用点击流数据创建一个概率性的单用户 ID。值得注意的是，用户在我们的过程中保持匿名。我们的目标不是用户的个人身份，而是更好地理解他们与我们产品的交互。
步骤 1:创建候选图
当匿名消费者使用 Zoopla 的酒店门户时，我们会像其他网站一样生成一个会话 ID。我们的目标是找到一个数据点，我们可以使用它将最有可能来自同一个用户的潜在会话 id 链接在一起。我们的目标是链接会话 IDs)跨设备以及(2)随着时间的推移。
我们用于链接会话 id 的数据是来自点击流的 IP 地址和时间戳:
我们使用时间戳等数据为 ML 模型创建附加功能。例如，我们在模型中使用时间戳的方式:
- 在数据准备中，过滤我们在特定时间窗口内观察到的 IP 和会话，以减少我们的数据总量，同时提高我们模型的质量
- 为了创建额外的功能，我们将图形嵌入与机器学习模型相结合。
第一步是创建我们的图形顶点。每个顶点是一个会话，对应的 IP 特征向量总结了我们在每个会话中看到的相关 IP，包括我们在会话中观察到的时间范围:
val vertices =
preparedClickstreamData
.groupBy("session_id", "ip")
.agg(
count("*").alias("ip_count"),
//Additional ts aggregations etc
)
.withColumn(“ip_count_struct”, struct("ip", "ip_count”))
.groupBy("session_id")
.agg(
//Spark SQL 2.4
map_from_entries(“ip_count_struct”).alias("ip_vector"),
//Additional ts aggregations etc
)
the vertices data model
我们还需要共享至少一个 IP 地址的候选会话图的边:
val udfFunc = udf((r: Seq[String]) => r.combinations(2).toSeq)val edges =
preparedClickstreamData
.groupBy("ip")
.agg(
collect_set("session_id").alias("session_id_set"),
//Additional ts aggregations
)
.withColumn("session_edges", explode(udfFunc($"session_id_set")))
.withColumn("src_session_id", col("session_edges")(0))
.withColumn("dst_session_id", col("session_edges")(1))
the edge data model
这已经为我们提供了一个相当有用的用户图表，其中包含了有见地的数据。为了保持简单并保持在 Spark 环境中，我们将使用 GraphX 作为我们的图形分析引擎。
虽然 Spark SQL 和 Dataframes 非常适合我们的数据准备，但 GraphX 使用 Spark RDDs。以下代码为我们的顶点和边属性创建数据模型，并为 GraphX 创建相应的 rdd:
import org.apache.spark.graphx._
import org.apache.spark.rdd.RDDcase class VertexAttributes(...)
case class EdgeAttributes(...)val vertexIdEncoder = (s: String) => java.util.UUID.fromString(s).getMostSignificantBits()
val udfVertexIdEncoder = udf(vertexIdEncoder)val vertexRdd: RDD[(VertexId, NodeAttributes)] = vertices
// handle collisions in production
.withColumn("vertexid", udfVertexIdEncoder(col("session_id")))
.rdd
.map(row => (row.getAs("vertexid"), VertexAttributes(
session_id=row.getAs("session_id"),
...
)))val edgeRdd: RDD[Edge[EdgeAttributes]] = edges
// handle collisions production
.withColumn("srcVertexId", udfVertexIdEncoder(col("src_session_id")))
.withColumn("dstVertexId", udfVertexIdEncoder(col("dst_session_id")))
.rdd
.map(row => Edge(
row.getAs("srcVertexId"), row.getAs("dstVertexId"),
EdgeAttributes(
ip=row.getAs("ip"),
...
)
))val graph = Graph(vertexRdd, edgeRdd)
graph.triplets.take(1)
我们可以遍历该图，并遍历具有相应顶点和边属性的所有图三元组:
这实现了我们的候选图。
步骤 2:过滤候选图
下一步，我们过滤候选图。根据您的使用情形和数据质量，这可能需要:
- 处理 vertexIdEncoder()函数中由长类型顶点 ID 转换引起的顶点 ID 冲突。
- 根据时间戳或会话/ ip 计数属性过滤顶点和边，例如潜在的公共 wifi 热点，以减少误报等。
- 从机器人和爬虫中移除顶点，即通过查看顶点度数。此外，顶点度数是机器学习模型的一个重要特征。使用 GraphX outerJoinVertices()用度数信息更新顶点属性，如下所示:
val vertexDegrees: VertexRDD[Int] = graph.degrees
val degreeGraph = graph.outerJoinVertices(vertexDegrees) {
(id, oldVertexAttr, vertexDegOpt) =>
vertexDegOpt match {
//Return a new vertex attribute case class with additional degree data
case Some(vertexDegOpt) => newVertexAttr(degree=vertexDegOpt, ...)
//No degree means zero edges
case None => newVertexAttr(degree=0, ...)
}
}
步骤 3:计算概率会话相似性
在最后一步中，我们计算候选图中会话顶点之间的相似性度量。
会话相似性的直接度量是相关会话 IP 向量的余弦相似性。假设 id1 和 id2 是两个给定的会话 id，具有以下观察到的 IP:
id1 -> [ip1, ip1, ip1, ip2, ip2]
id2 -> [ip1, ip1, ip2, ip2, ip2]
考虑到空间，我们可以将 IP 向量及其相似性表示为:
id1 -> v1 = (3, 2)
id2 -> v2 = (2, 3)
cosine_sim(v1, v2) = cos(ϴ) = v1.v2 / (|v1||v2|)
cosine_sim(v1, v2) = 0.92
捕获顶点度数后，我们可以选择折叠共享多个 IP 的会话中的重复边。这些只会创建具有相同余弦相似性的三元组。这移除了冗余的图形三元组，并加速了顶点相似性计算:
val dedupedGraph = graph.groupEdges(
(edgeAttr1, edgeAttr2) => newEdgeAttr(...)
)
为了通过余弦会话相似性过滤我们的候选图，我们使用 GraphX 子图()来过滤边三元组。首先，我们为组合的 IP 向量空间实现余弦相似度:
object CosineSim {
def dotProduct[A](l: Map[A, Long], r: Map[A, Long]): Long = {
val k = (l.keySet ++ r.keySet).toSeq
val v1 = k.map(l.getOrElse(_, 0L))
val v2 = k.map(r.getOrElse(_, 0L))
(v1 zip v2).map { case (x, y) => x * y }.sum
}
def vectorNorm[A](v : Map[A, Long]): Double = Math.sqrt(dotProduct(v, v))
def cosineSim[A](l: Map[A, Long], r: Map[A, Long]): Double = {
val d = dotProduct(l, r)
val n1 = vectorNorm(l)
val n2 = vectorNorm(r)
if (n1 == 0 || n2 == 0)
0.0
else
d.toDouble / n1 / n2
}
}
并且简单地过滤余弦相似度高于期望阈值的图形三元组:
val cosineFilteredGraph = graph.subgraph(
epred = triplet => CosineSim.cosineSim(triplet.srcAttr.ip_vector, triplet.dstAttr.ip_vector) > threshold
)
通过这些简单的步骤，你就有了第一个版本的概率单用户 ID 图。为了访问单个用户连接的会话，我们从图中提取连接的组件:
val singleUserGraph = cosineFilteredGraph.connectedComponents()
val singleUserComponents = singleUserGraph.vertices
我们可以查找单个会话 ID:
// Lookup a single user ID in the graph for a session
val sessionId = “fa199794-0c60-4088-bccd-858c67a954ed”
val vertexId = vertexIdEncoder(sessionId)
val singleUserId: VertexId = singleUserComponents.filter {
case (id, component) => id == vertexId
}.map(_._2).collect().head
// Get all other sessions connected to this single user ID
singleUserComponents.filter {
case (id, component) => component == singleUserId
}.map(_._1).collect
或者收集整个数据集用于进一步处理:
// Collect the data
singleUserComponents.aggregateByKey(Seq[VertexId]())(
(accum:Seq[VertexId], v:VertexId) => accum :+ v,
(v1:Seq[VertexId], v2:Seq[VertexId]) => v1 ++ v2
)
GraphX 是一个非常高效的图形分析框架，适用于单跳类型的问题，比如计算我们的余弦相似度。但是对于更复杂的查询和图形，例如从图形中提取连接的组件，性能会迅速下降。下面的 DAG 是 Spark 程序，用于创建本文中讨论的单个用户 id:
用新数据逐步更新图表也不是一个简单的问题。但这为我们提供了一个批量解决方案来改进我们的推荐系统，并提供了一个丰富的数据模型来为各种批量 ML 模型创建功能。
信用
Zoopla 的数据科学是一个伟大的团队努力，我站在巨人的肩膀上。特别感谢 Zoopla 的高级数据科学家 Behzad Behzadan ，他对我们的用例进行了最初的研究工作，并且是这个项目的主要个人贡献者。
Jan 是公司数据转型方面的成功思想领袖和顾问，拥有将数据科学大规模应用于商业生产的记录。他最近被 dataIQ 评为英国 100 位最具影响力的数据和分析从业者之一。
在领英上连接:https://www.linkedin.com/in/janteichmann/
阅读其他文章:https://medium.com/@jan.teichmann
Python 中的奇异值分解示例
https://www.pexels.com/photo/woman-writing-on-a-whiteboard-3862130/
奇异值分解(SVD)有着广泛的应用。这些包括降维、图像压缩和数据去噪。本质上，SVD 表明一个矩阵可以表示为其他三个矩阵的乘积。用数学术语来说，SVD 可以写成如下形式:
其中 n 为行数(即样本数) p 代表维数。
假设我们有一个矩阵 一个 。
为了确定关联的矩阵，我们必须先求出矩阵的特征向量乘以矩阵的转置。**
因为
回想一下特征向量和特征值的定义:
后者也可以表示为:
因此，回到我们的例子:
我们得到了一个四次多项式。
求解后，我们用其中一个特征值代替λ。
在将矩阵乘以向量 x 之后，我们得到如下:
在解方程时，我们得到:
我们把特征向量代入的列。****
然后，我们对矩阵乘以矩阵的转置重复同样的过程。**
求解后，我们得到了 V 的表达式。
最后， S 是一个对角矩阵，其值是任一个的特征值的平方根
这就是事情变得有趣的地方。我们知道所有三个矩阵的乘积等价于左边的矩阵。
我们可以排除特征，但仍然保持原始矩阵的近似值。假设矩阵是组成图像的列和行或像素的数据集，我们可以在理论上使用新形成的矩阵来训练模型，并达到相当的(如果不是更好的话)(由于维数灾难)精确度。
密码
让我们看看如何在 Python 中应用奇异值分解。首先，导入以下库。
***import numpy as np
from sklearn.datasets import load_digits
from matplotlib import pyplot as plt
from sklearn.decomposition import TruncatedSVD
float_formatter = lambda x: "%.2f" % x
np.set_printoptions(formatter={'float_kind':float_formatter})
from sklearn.ensemble import RandomForestClassifier***
在接下来的教程中，我们将尝试对手写数字进行分类。幸运的是，scikit-learn库提供了一个包装器函数，用于将数据集导入我们的程序。
***X, y = load_digits(return_X_y=True)***
该数据集包含 1797 幅 8×8 的图像。如果你指定了return_X_y=True,这个函数将把像素作为一维数组返回。
***X.shape***
y 包含每个数字的标签。
**y**
让我们来看看第一个数字。正如我们所看到的，它只是一个长度为 64 的数组，包含像素亮度。
**image = X[0]**
如果我们想使用matplotlib查看图像，我们必须首先重塑数组。
**image = image.reshape((8, 8))plt.matshow(image, cmap = 'gray')**
接下来，我们将使用奇异值分解来看看我们是否能够仅使用每行的两个特征来重建图像。函数返回的 s 矩阵必须使用diag方法转换成对角矩阵。默认情况下，diag将创建一个相对于原始矩阵为 n x n 的矩阵。这导致了一个问题，因为矩阵的大小不再遵循矩阵乘法的规则，其中一个矩阵中的列数必须与另一个矩阵中的行数相匹配。因此，我们创建一个新的 m x n 矩阵，并用对角矩阵填充它的第一个 n x n 部分。
**U, s, V = np.linalg.svd(image)S = np.zeros((image.shape[0], image.shape[1]))S[:image.shape[0], :image.shape[0]] = np.diag(s)n_component = 2S = S[:, :n_component]
VT = VT[:n_component, :]A = U.dot(Sigma.dot(VT))print(A)**
**plt.matshow(A, cmap = 'gray')**
我们可以通过取 U 和 S 矩阵的点积得到缩减的特征空间。
**U.dot(S)**
原始与缩减的特征空间
让我们比较随机森林模型在使用原始手写数字训练时和使用从奇异值分解获得的缩减特征空间训练时的准确性。
我们可以通过查看公开得分来衡量模型的准确性。如果你对 OOB 的概念不熟悉，我鼓励你看看兰登森林的这篇文章。
**rf_original = RandomForestClassifier(oob_score=True)rf_original.fit(X, y)rf_original.oob_score_**
接下来，我们用 2 个组件创建并装配一个TruncatedSVD类的实例。值得一提的是，与前面的例子不同，我们使用的是 2/64 特性。
**svd = TruncatedSVD(n_components=2)X_reduced = svd.fit_transform(X)**
精简数据集中的每个图像(即行)包含 2 个特征。
**X_reduced[0]**
看一看图像，很难区分图像由什么数字组成，它很可能是 5 而不是 0。
**image_reduced = svd.inverse_transform(X_reduced[0].reshape(1,-1))image_reduced = image_reduced.reshape((8,8))plt.matshow(image_reduced, cmap = 'gray')**
在精简的数据集上训练随机森林分类器后，我们获得了 36.7%的微弱准确率
**rf_reduced = RandomForestClassifier(oob_score=True)rf_reduced.fit(X_reduced, y)rf_reduced.oob_score_**
我们可以通过取explained_variance_ratio_属性的和得到总方差解释。我们通常希望达到 80%到 90%的目标。
**svd.explained_variance_ratio_.sum()**
让我们再试一次，只是这一次，我们使用 16 个组件。我们查看包含在 16 个特征中的信息量。
**svd = TruncatedSVD(n_components=16)X_reduced = svd.fit_transform(X)svd.explained_variance_ratio_.sum()**
我们获得了与使用原始图像训练的模型相当的准确度，并且我们使用了 16/64=0.25 的数据量。
**rf_reduced = RandomForestClassifier(oob_score=True)rf_reduced.fit(X_reduced, y)rf_reduced.oob_score_**
用于数字分类的正弦神经网络
一种正弦基函数的隐层神经网络
我在之前的帖子里讲过正弦基函数的神经网络(正弦基函数的神经网络)。这是具有正弦基函数的神经网络的基本实现。此后，我将把具有正弦基函数的神经网络称为“正弦神经网络(SNN)”。这篇文章假设你了解神经网络及其工作原理。如果你不知道，有惊人的故事来学习什么是媒体中的神经网络。
A fancy image by Gertrūda Valasevičiūtė on Unsplash
今天，我将介绍为数字分类而实现和训练的 SNNs。数字分类是神经网络应用的一个众所周知的领域。但是，对新结构进行首次评估是一个很好的主题。
另一个不同的主题，我也将谈到的是实现 SNN 与隐藏层。当我从零开始实现神经网络架构时，我严重地遇到了很多矩阵维数的问题，尤其是哪些矩阵应该在哪些维数上才能进行正确的计算。每当我试图实现带有隐藏层的神经网络时，我都会对此感到困惑。既然我终于自己说清楚了，那我也帮你理解一下 NNs 中的维数，矩阵及其乘法。您还将认识到，snn 不需要任何激活函数来创建线性，因为它们本身就是非线性函数。
Python 代码中的架构、维度和公式
Sinusoidal Neural Networks — Structure and Formulas
前馈操作非常简单。我们只是将矩阵相乘并求和(点积。)
反向传播是令人头疼的问题。你可能也猜到了，我们将在反向传播中使用链式法则。问题是矩阵的维数。每个输出应该有哪个维度？我一时无法回答这个问题。我终于明白了，恭喜我，😄
矩阵的维数
在这一节中，dy_dx 表示 y 相对于 x 的导数。因此，要注意变量的命名，以理解正在发生的事情。
derivative of y with respect to x
derror_doutput 具有我们的输出变量的维数，(1xK)。
doutput_dlayer 通常需要与我们的隐藏层具有相同的维度。但是，我们不会让它发生。我们将把它保留为 KxN，因为这个变量将进入一个带有 derror_doutput 的操作。此操作的输出应该与隐藏层具有相同的维度。
derror_dlayer 具有我们隐藏层的维度。(详见上文。)
dlayer _ dweights_in_1 将与 weights _ in _ 1 具有相同的尺寸。因为 weights_in_1 矩阵中的每个变量都应该变化。
derror_dweights_in_1 将具有 weights_in_1 矩阵的转置维数。当我们更新权重时，我们把它固定回来(转置它)。
doutput _ dweights_in_2 将与 weights _ in _ 2 具有相同的维度。因为 weights_in_2 矩阵中的每个变量都应该改变。
derror_dweights_in_2 将具有 weights_in_2 矩阵的转置维数。当我们更新权重时，我们把它固定回来(转置它)。
我没有写 weights_out_1 和 weights_out_2，因为它们的维数与 weights_in_1 和 weights_in_2 相同。
哦，终于实现了…
构建神经网络类:
import numpy as npclass NeuralNetwork:
#DEFINE INITIALS
def __init__(self, x, y):
self.input = x
self.y = y
self.layer = np.zeros((1,32))
#weights connecting to hidden layer
self.weights_in_1 = np.random.uniform(-1, 1,(self.layer.shape[1],self.input.shape[1]))
self.weights_out_1 = np.random.uniform(-1, 1,(self.layer.shape[1],self.input.shape[1]))
#weights connecting to output.
self.weights_in_2 = np.random.uniform(-1, 1,(self.y.shape[1],self.layer.shape[1]))
self.weights_out_2 = np.random.uniform(-1, 1,(self.y.shape[1],self.layer.shape[1]))
self.output = np.random.uniform(-1, 1,self.y.shape)
print('Output:',self.output)
print('Y:',self.y)
#DEFINE FEEDFORWARD
def feedforward(self,cnt):
self.layer = np.sum(self.weights_out_1*self.input[cnt]*np.sin(np.pi*self.input[cnt]*self.weights_in_1),axis = 1) #forwarding to hidden layer
self.output[cnt] = np.sum(self.weights_out_2*self.layer*np.sin(np.pi*self.layer*self.weights_in_2),axis = 1) #forwarding to output
#function for derivative of output with respect to hidden layer.
def doutput_dlayer_func(self):
return self.weights_out_2*np.sin(np.pi*self.layer*self.weights_in_2)+np.pi*self.layer*self.weights_in_2*self.weights_out_2*np.cos(np.pi*self.layer*self.weights_in_2)
#function for derivative of hidden layer with respect to weights on first level.
def dlayer_d_weights_func(self,cnt):
doutput_dweights_in = self.weights_out_1 * np.square(self.input[cnt]) * np.pi * np.cos(np.pi*self.input[cnt]*self.weights_in_1)
doutput_dweights_out = self.input[cnt]*np.sin(np.pi*self.input[cnt]*self.weights_in_1)
return doutput_dweights_in,doutput_dweights_out #function for derivative of output with respect to weights on second level.
def doutput_d_weights_func(self):
doutput_dweights_in = self.weights_out_2 * np.square(self.layer) * np.pi * np.cos(np.pi*self.layer*self.weights_in_2)
doutput_dweights_out = self.layer*np.sin(np.pi*self.layer*self.weights_in_2)
return doutput_dweights_in,doutput_dweights_out#DEFINE BACKPROPAGATION
def backprop(self,cnt):
error = np.square(self.y[cnt]-self.output[cnt])
print(cnt,'___',np.sum(error))
derror_doutput = self.y[cnt]-self.output[cnt]
#print('___',cnt,'___',derror_doutput)
#calculate update amount for weights_in_1 and weights_out_1
#application of chain rule for weights on first level.
doutput_dlayer = self.doutput_dlayer_func()
derror_dlayer = np.dot(derror_doutput,doutput_dlayer)
dlayer_dweights_in_1,dlayer_dweights_out_1 = self.dlayer_d_weights_func(cnt)
derror_dweights_in_1 = derror_dlayer*dlayer_dweights_in_1.T
derror_dweights_out_1 = derror_dlayer*dlayer_dweights_out_1.T
#application of chain rule for weights on first level.
doutput_dweights_in_2, doutput_dweights_out_2 = self.doutput_d_weights_func()
derror_dweights_in_2 = derror_doutput*doutput_dweights_in_2.T
derror_dweights_out_2 = derror_doutput*doutput_dweights_out_2.T
self.weights_in_1 += derror_dweights_in_1.T*0.0001
self.weights_out_1 += derror_dweights_out_1.T*0.0001
self.weights_in_2 += derror_dweights_in_2.T*0.0001
self.weights_out_2 += derror_dweights_out_2.T*0.0001
#PREDICT THE TEST DATA
def feedforward_test(self,test_data):
self.layer = np.sum(self.weights_out_1*test_data*np.sin(np.pi*test_data*self.weights_in_1),axis = 1)
test_output = np.sum(self.weights_out_2*self.layer*np.sin(np.pi*self.layer*self.weights_in_2),axis = 1)
return test_output
def predict(self,input_):
predictions = []
for elm in input_:
#print('___',elm)
predictions.append(self.feedforward_test(elm).tolist())
return np.array(predictions)
#SAVE WEIGHTS
def save_weights(self,dir_in_1 = './weights_in_1.npy',dir_out_1 = './weights_out_1.npy',
dir_in_2 = './weights_in_2.npy',dir_out_2 = './weights_out_2.npy'):
np.save(dir_in_1,self.weights_in_1)
np.save(dir_out_1,self.weights_out_1)
np.save(dir_in_2,self.weights_in_2)
np.save(dir_out_2,self.weights_out_2)
#IMPORT WEIGHTS
def import_weights(self,dir_in_1 = './weights_in_1.npy',dir_out_1 = './weights_out_1.npy',
dir_in_2 = './weights_in_2.npy',dir_out_2 = './weights_out_2.npy'):
self.weights_in_1 = np.load(dir_in_1)
self.weights_out_1 = np.load(dir_out_1)
self.weights_in_2 = np.load(dir_in_2)
self.weights_out_2 = np.load(dir_out_2)
测试:
from sklearn.datasets import load_digits
from sklearn import preprocessinglb = preprocessing.LabelBinarizer() #create label binarizer
digits_object = load_digits()
images = digits_object.data
images = (images/np.max(images))+0.01
labels = digits_object.target.reshape((1797,1))
lb.fit(labels) # fit your labels to a binarizing map.
labels = lb.transform(labels) #binarize your labels (example: from [0,3,2] to [[1,0,0,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0]])#split your data-set
x = images[0:1500]
test_x = images[1500:-1]
y = labels[0:1500]
test_y = labels[1500:-1]#create neural network instance.
nn = NeuralNetwork(x,y)#TRAIN NETWORK FOR 1000 TIMES (EPOCH = 1000).
for gen_cnt in range(1000):
for cnt in range(1500):
print('Epoch: {}'.format(gen_cnt))
nn.feedforward(cnt)
nn.backprop(cnt)#PREDICT THE TEST DATA
predictions = nn.predict(test_x)
error = np.around(predictions,decimals = 2)-test_y#CHECK BELOW
结果和结论:
我不知道为什么在神经网络中不常用正弦基函数。因为结果其实是令人满意的。我刚刚使用了 1 个隐藏层，得到了相当好的结果。我们不需要为选择不同类型的激活函数而感到困惑。非线性已经位于正弦基函数之下。我们的神经网络应该做的唯一一件事就是选择要求和的正弦函数。
编辑:我修改了预测函数，以便能够看到准确性。您可以在下面找到新功能:
def predict(self,input_):
predictions = []
for elm in input_:
#print('___',elm)
result_local = self.feedforward_test(elm)
result_temp = np.zeros(result_local.shape)
result_temp[np.argmax(result_local)] = 1
predictions.append(result_temp.tolist())
return np.array(predictions)
在运行预测函数后也添加这个脚本。
accuracy = 1-((np.count_nonzero(predictions-test_y)/2)/predictions.shape[0])
print('Accuracy is {}'.format(accuracy*100))
我在 1000 次迭代后检查了准确度。准确率为 86.82%。这可以通过添加更多的隐藏层或通过用更多的纪元训练网络来改善。(训练集:[0，1500]，测试集:[1500，1797])
每位数据科学家将面临的六大挑战以及如何克服这些挑战
Supplements for thought? Image made by author.
信息时代给人类带来了有史以来最大的科技类工作爆炸之一。虽然优步、脸书、AirBnB 和亚马逊等成功大公司背后的强大力量是他们的独创性和对消费者的便利，但他们的成功也可以归功于对数据的利用——而且是大量数据。不可否认，数据现在是大家趋之若鹜的新资源。甚至像可口可乐和百事可乐这样的公司也在利用消费者数据来制定更好的营销活动和产品战略决策。这导致数据被称为“新石油”，这句话最初是由数学家克莱夫·亨比说的，后来被经济学家轰动一时。虽然许多人在字面上和比喻上都不同意这种说法，但重点是显而易见的:数据具有非凡的价值，并且在未来只会变得更有价值。理解如何在充满数据的世界中导航对于维护人类和地球其他部分的社会、经济、环境和身体健康至关重要。数据科学家可以通过将大量与数据相关的模糊方面分解成清晰的术语来简化这一过程。因此，我想了一些方法来提高所有人的数据科学意识，同时坚持真正的科学。
1。)数据是补充，不是饭。
数据是有价值的，可以提供大量的信息。这些信息可以以多种方式用于澄清或混淆当今世界的许多问题。例如，在可再生能源的日常运营中，独立服务运营商(ISOs)进行预测，以预测在假设的天气条件下某一天将使用多少能源。这有助于运营商更好地管理电力资源，并提供一个清晰的估计，即一天需要多少能源，以及它将来自哪里。然而，天气条件可能会快速变化，在许多情况下，这些预测成为最佳估计，更多的能源最终被使用或取自原本不会像它们那样被利用的来源。
Image captured from the California ISO webpage (http://www.caiso.com/TodaysOutlook/Pages/default.aspx)
这种情况下的数据被用作理解感知条件的手段，但它不是实际情况。因此，为什么数据应该用于做出更好的决策，但不应该完全依赖于整个故事。在复杂系统中尤其如此，在复杂系统中，元素可以被迅速改变，结果也可以被彻底改变。对于数据科学家来说，培养对手头问题的良好直觉和本能仍然是很好的实践。重要的是要知道复杂系统的模型只有在 95–99%的置信度下才是好的。有必要采取适当的措施来应对不确定性。
2。)了解生态系统
许多数据科学家喜欢一头扎进问题中，认为因为他们对模型有很好的理解，他们可以在任何领域产生有意义的结果。这一点与前者有关。数据科学家实际了解他们研究领域的理论和概念知识是非常重要的。如果你是一名训练有素的数据科学家，对空气动力学一无所知，但你接受了一份在飞机上创建模型检查机翼应力数据的工作，你将很多人置于危险之中。围绕航空学和航空学有太多的知识。您可能认为您正在构建一个不错的模型，但实际上您只是通过一个基于感知器的算法来输入变量并产生随机结果。你不会真的知道如何有效地传达你从模型中得到的信息。
Perceptron does not equate to instant weather forecasting. Image created by the author.
关键要点是，数据科学家是充满活力的个体，但我们需要彻底了解我们将要工作的生态系统。如果你要建立一个雷暴预测的模型，与你要建立一个股票市场萧条和繁荣的模型相比，你需要一个完全不同的概念工具集。这迫使有抱负和漫游的数据科学家成为世界上最大的自学者。
3。)仅仅因为数据可用就使用它不是一个好主意。
可以理解。你已经在办公桌前用拳头敲打木头台面好几天了，因为你找不到任何数据来完善你的模型。当你感到希望的机会慢慢溜走，就像刚炸好的鸡肉上的油脂一样，你在谷歌上找到了一个看起来最理想的数据库。你下载数据集，并做你通常的探索技术。散点图、箱线图、相关矩阵和一个微小的线性判别分析测试之后，您决定这是您唯一的机会，您将在您的模型中使用它。问题是，你后来发现数据集是由大卫·杜伊(David Duey)制作的，他是一位虚构的私人顾问，他对另一个不完整的数据集应用了一种“转换算法”，从而产生了你一直在搜索的数据集。显然，这其中的问题是众多而明显的。这个数据集很可能是无用的，如果它曾经服务于任何目的，它就不应该在产生后被传播。
这里的一个关键提示是，不管你的数据来自哪里，总是要问数据本身无法回答的问题。数据是从哪里来的？数据是谁做的？数据是怎么做出来的？该数据库公开多久了？还有谁在自己的研究中使用了这些数据？问这些问题不仅会让你的模型更好，而且从长远来看，会让你成为一名更好的数据科学家。
4。)一种模式不适合所有人。
*虽然数据确实可以以多种方式使用，但没有一种模式适合所有情况。是的，只要假设得到满足，数据得到正确处理，方法就可以应用于不同的数据集。然而，使用一种先验算法完全是浪费时间，而这种算法对于预测杂货店的销售模式是非常准确的，例如降雨模式。首先，杂货店的“*定律”与热力学定律的运作方式不同。仅仅因为你的模型显示买牛奶的人倾向于买黄油，并不意味着仅仅因为天气热，一天的晚些时候就会下雨。还有许多其他变量可以用来预测降雨事件。
远离毁灭的建议，总是在数据集上测试不同模型的功效。可以制作基线模型来测量它们与旧模型相比的准确性。如果您发现精度和准确性有任何妥协，或者如果您怀疑大量的准确性，这可能是一个好迹象，表明您的模型不是最适合您的数据。
5。)要知道你的大脑是可以忽悠你的。
Do you see a pattern here? I don’t. Image taken from google.
很容易犯这个错误。你得到数据并对其进行测试，你似乎已经发现了一些关于你所研究的东西的关键信息。在某些情况下，这可能是一个实际的结果，但在大多数情况下，你可能会欺骗自己相信你已经找到了一个不存在的模式。人类的大脑想要找到模式，并以此为生。整天看着曲线、线条、点和像素是欺骗自己相信你发现了一个实际上并不存在的模式的快速方法。在庆祝之前，一定要反复运行你的模型，彻底检查你的数据是否有任何问题。你可能会意识到你错过了一条重要的信息。
6。)总是质疑结果。
不管 p 值或模型预测的准确性如何，您都应该对最终结果持怀疑态度。这也应该在你的报告或对结果的讨论中表现出来。即使您的结果似乎始终如一地为您提供信息，改进您正在尝试改进的系统的某个方面，但重要的是要认识到随机性仍然存在，并可能导致您的输出出现严重偏差。你应该小心信息的内容，以及它如何在更广泛的受众手中被歪曲。采取这一步骤可确保您对您的数据做出符合道德的决定，并防止以后对整个社会造成灾难性事件。
Question data and the results they provide. Image taken from google.
记住这些提示，我认为数据科学领域内外的个人都可以获得关于数据的新观点。即使你不直接处理数据，问一些很少有人会问的关于给定数据集的结果的问题也是很重要的。如果你是一名数据科学家，重要的是要问清楚你的模型及其结果将如何影响那些可能错误解读或出于个人原因歪曲信息的人的行动。数据是有价值的，在未来几年，它将继续在我们的生活中发挥重要作用。让每个人更多地意识到它的局限性、优点和缺点，对于创造一个更实际的社会是很重要的。
给有抱负的数据科学家的六条建议
Source: https://www.maxpixel.net/Art-Colourful-Gears-Creativity-Cogs-Colorful-1866468
找到工作前积累经验
数据科学是一个需求巨大的领域，部分原因是它似乎需要作为数据科学家的经验才能被聘为数据科学家。但我合作过的许多最优秀的数据科学家都有从人文科学到神经科学的不同背景，要脱颖而出需要经验。作为一名即将进入数据科学职业生涯的新毕业生或分析专业人士，构建一个工作组合来展示该领域的专业知识可能是一项挑战。我在数据科学职位的招聘过程中经历过两种情况，我想列举一些能够帮助我找到数据科学家工作的关键经历:
- 亲身体验云计算
- 创建新的数据集
- 把东西粘在一起
- 支持一项服务
- 创造惊人的视觉效果
- 写白皮书
我将详细阐述这些主题，但数据科学的关键主题是能够构建为公司增加价值的数据产品。能够构建这些端到端数据产品的数据科学家是一笔宝贵的财富，在从事数据科学职业时展示这些技能是有用的。
亲身体验云计算
许多公司都在寻找过去在云计算环境中有经验的数据科学家，因为这些平台提供了支持数据工作流和预测模型扩展到海量数据的工具。你还可能在日常工作中使用云平台，比如亚马逊网络服务(AWS)或谷歌云平台(GCP)。
好消息是，这些平台中的许多都提供了免费的层来熟悉平台。例如，AWS 有免费的 EC2 实例，并免费使用 Lambda 等服务来满足少量请求，GCP 提供 300 美元的免费积分来试用该平台的大部分内容，Databricks 提供了一个社区版，您可以使用它来接触该平台。有了这些免费选项，你将无法处理大量数据集，但你可以在这些平台上积累经验。
我的一个建议是在这些平台上尝试不同的特性，看看是否可以使用一些工具来训练和部署模型。例如，在我的模型服务文章中，我利用了我已经熟悉的工具 SKLearn，并研究了如何将模型包装成 Lambda 函数。
我的创业数据科学系列的第二部分主要关注 Python。
towardsdatascience.com](/data-science-for-startups-model-services-2facf2dde81d)
创建新的数据集
在学术课程和数据科学竞赛中，您通常会得到一个干净的数据集，其中项目的重点是探索性数据分析或建模。然而，对于大多数现实世界的项目，您需要执行一些数据整理，以便将原始数据集整理成对分析或建模任务更有用的转换数据集。通常，数据管理需要收集额外的数据集来转换数据。例如，为了更好地理解美国富裕家庭的资产配置，我过去曾在美联储数据公司工作过。
美国最近通过的税改法案引发了许多关于该国财富分配的问题…
medium.freecodecamp.org](https://medium.freecodecamp.org/clustering-the-top-1-asset-analysis-in-r-6c529b382b42)
这是一个有趣的项目，因为我使用第三方数据来衡量第一方数据的准确性。我的第二个建议是更进一步，建立一个数据集。这可以包括抓取网站、从端点采样数据(例如 steamspy )或将不同的数据源聚合到新的数据集。例如，我在研究生学习期间创建了一个自定义的星际争霸回放数据集，它展示了我在一个新的数据集上执行数据管理的能力。
2009 年，我发表了一篇关于预测《星际争霸:育雏战争》建造顺序的论文，使用了不同的分类…
towardsdatascience.com](/reproducible-research-starcraft-mining-ea140d6789b9)
把东西粘在一起
我喜欢看到数据科学家展示的技能之一是让不同的组件或系统协同工作以完成任务的能力。在数据科学的角色中，可能没有一个清晰的模型产品化的路径，您可能需要构建一些独特的东西来启动和运行系统。理想情况下，数据科学团队将获得工程支持来启动和运行系统，但原型制作是数据科学家快速行动的一项重要技能。
我的建议是尝试将不同的系统或组件集成到数据科学工作流中。这可能涉及到动手使用工具，如气流，以原型数据管道。它可以包括在不同系统之间建立一座桥梁，比如我开始用 Java 连接星际争霸育雏战争 API 库的 JNI-BWAPI 项目。或者它可以涉及在一个平台内将不同的组件粘合在一起，例如使用 GCP 数据流从 BigQuery 中提取数据，应用预测模型，并将结果存储到云数据存储中。
我正在进行的关于在创业公司建立数据科学学科系列的第七部分。您可以找到所有……
towardsdatascience.com](/data-science-for-startups-model-production-b14a29b2f920)
支持一项服务
作为一名数据科学家，您经常需要提供公司内其他团队可以使用的服务。例如，这可能是一个提供深度学习模型结果的 Flask 应用程序。能够原型化服务意味着其他团队将能够更快地使用您的数据产品。
这篇文章演示了如何使用 Keras 构建的深度学习模型来设置端点以服务于预测。它…
towardsdatascience.com](/deploying-keras-deep-learning-models-with-flask-5da4181436a2)
我的建议是获得使用诸如 Flask 或 Gunicorn 等工具的实践经验，以便设置 web 端点，以及 Dash 以便用 Python 创建交互式 web 应用程序。尝试在一个 Docker 实例中设置这些服务之一也是一个有用的实践。
创造惊人的视觉效果
虽然伟大的作品应该是独立的，但在解释为什么一个分析或模型是重要的之前，通常有必要首先引起你的观众的注意。我在这里的建议是学习各种可视化工具，以创建引人注目的突出可视化。
自从我开始为《星际争霸:育雏战争》进行数据挖掘专业回放以来，已经过去了将近十年。最近之后…
towardsdatascience.com](/visualizing-professional-starcraft-with-r-598b5e7a82ac)
创建可视化也是建立作品组合的一种有用方式。下面的博文展示了我作为数据科学家 10 多年来探索的不同工具和数据集的样本。
我在数据科学领域的职业生涯始于十年前，当时我在加州大学圣克鲁斯分校(UC Santa Cruz)上了第一门机器学习课程。自从…
towardsdatascience.com](/10-years-of-data-science-visualizations-af1dd8e443a7)
写白皮书
我最近一直提倡的数据科学技能之一是以白皮书的形式解释项目的能力，该白皮书提供了执行摘要，讨论了如何使用该工作，提供了有关方法和结果的详细信息。我们的目标是让你的研究能够被广泛的受众所理解，并且能够自我解释，以便其他数据科学家可以在此基础上进行研究。
写博客和其他形式的写作是获得提高书面交流经验的好方法。我在这里的建议是，尝试为广大受众撰写数据科学文章，以便获得在不同细节层次传达想法的经验。
数据科学家写书有很多令人信服的理由。我想更好地理解新工具，而且…
towardsdatascience.com](/data-science-for-startups-blog-book-bf53f86ca4d5)
结论
数据科学需要大量工具的实践经验。幸运的是，这些工具越来越容易获得，构建数据科学组合也变得越来越容易。
本·韦伯是 Zynga 公司的首席数据科学家，也是 T2 恶作剧公司的顾问。
六适马盖奇 R&R 分析与 R
r 代表工业工程师
探索“六西格玛”R 包
Image by William Warby available at Unsplash
测量系统分析
测量系统分析(MSA)，也称为量具 R&R(重复性和再现性)，是 DMAIC 循环中测量阶段最重要的活动。它识别并量化影响测量系统的变异源。测量系统的可变性不受控制的过程是无法改进的。一个好的测量系统只有随机可变性，这是由于被测项目的固有变化，而不是由评估者(操作员、机器等)产生的。)、零件、时间或其他因素。与此主题相关的一些重要术语包括以下概念:
- 精度:在规定条件下获得的独立测试结果之间的一致程度。
- 重复性:重复性条件下的精度(在短时间间隔内，由同一操作人员使用同一设备，在同一实验室对同一检测项目使用同一方法获得的独立检测结果)。
- 再现性:再现性条件下的精度(不同实验室不同操作员使用不同设备对相同项目使用相同方法获得测试结果)。
总之，可重复性可定义为测量系统的固有可变性(在类似条件下)，再现性可定义为不同条件(组)下的可变性，例如操作员、机器或任何其他因素。
Variability decomposition
评估测量系统
为了使测量系统准确，量具 R&R 可变性对总可变性的贡献应低于 10%。介于 10%和 30%之间的值被认为是可接受的。然而，大于 30%的值代表低效的测量系统。
为了评估这一事实，使用了可变性的平方根(即标准偏差)。然后计算每个变异源在整个研究中的贡献 Var:
同样，另一个有用的度量是不同类别的数量。它使用下面的公式计算，并舍入到最接近的较低整数(如果小于 1，则为 1)。
不同类别的数量应大于或等于 4。该值测量测量系统引起的可变性和固有可变性之间的关系。如果低于 4，则与固有可变性相比，计量器 R&R 可变性较大。否则，可以认为两个可变性之间的关系是适当的。
来自 SixSigma R 包的 ss.rr 函数自动执行量具 R & R 研究。对于下面的例子，让我们考虑这样一种情况，电池制造商用两个不同的电压表(评估员)，三个不同的时间(重复)测量三个不同电池(部件)的电压。我们来看看 R 代码！
Six Sigma Gage R&R R Code
Charts for MSA of battery example
根据 MSA 的结果，由 R&R 引起的%StudyVar 为 93.13%。不同类别的总数为 1。这种结果组合(即类别数量少，可变性百分比大)代表了测量系统的最坏情况。
图表解读
- 左上角的图表显示了每个组成部分对总方差的贡献。其目的是检测测量问题是来自重复性还是再现性。
- 中间和左下角的图表代表均值控制图和极差控制图，其相应的控制限适用于 R&R 研究。在均值控制图中，所有点都应在控制极限之外，因为该图代表了零件与零件之间的差异(相同的操作员，相同的原型)。否则，测量系统将被认为是不充分的。另一方面，在范围控制图中，所有点都应位于控制界限内。
- 右上方和中间的图表显示了图表中的每个测量点。在右上方的图表中， x- 轴代表电池，而绘制的线条连接了每个原型的平均值。在右中图中， x- 轴分别代表电压表。
- 右下角的图表显示了两个因素之间的相互作用:操作符和原型。成对电池*电压表的平均值用点表示，用线连接。好的图表有助于检测操作者和原型之间是否有任何交互，或者操作者之间是否有差异。
根据图表结果，很明显测量系统没有正常工作。这个例子中的主要问题是电压表(鉴定人)。
总结想法
测量是质量保证中的一项重要任务。错误测量的数据可能导致错误的决策，从而导致不充分的结果。因此，MSA 应作为任何六个适马项目的初始阶段进行。MSA 允许工程师分析数据，确定来源和变化。有了这些信息，他们就可以继续规划和开发改善活动，以期获得更高的收益。 SixSigma R 包是一个很棒的工具，它只需要几行代码就可以执行 gage R & R 分析。非常鼓励工业、过程和质量工程师利用这一工具，以高质量和高效率的标准满足客户的要求。
—
如果你觉得这篇文章有用，欢迎在 GitHub 上下载我的个人代码。你也可以直接在 rsalaza4@binghamton.edu 给我发邮件，在LinkedIn上找到我。有兴趣了解工程领域的数据分析、数据科学和机器学习应用的更多信息吗？通过访问我的媒体 简介 来探索我以前的文章。感谢阅读。
——罗伯特
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
