GAMEBOY 超级电脑
每秒超过 10 亿帧的总速度可以说是世界上最快的 8 位游戏主机集群。
Distributed Tetris (1989)
1 of 48 boards of IBM Neural Computer which was used for experiments
怎么建一个？
食谱
拿一把硅，倒强化学习知识、超算经验和计算机架构激情，然后加上汗水和泪水，最后搅拌大约 1000 小时，直到沸腾瞧吧……
为什么会有人想要一个呢？
简答:向强人工智能迈进。
这是更长的版本:
现在是 2016 年。深度学习无处不在。图像识别可以被认为是由卷积神经网络解决的，我的研究兴趣被吸引到具有记忆和强化学习的神经网络。
具体来说，在 Google Deepmind 展示的一篇论文中，已经表明使用一种叫做 Deep Q-Neural Network 的简单强化学习算法，在各种 Atari 2600(1977 年发布的家用游戏机)游戏上实现人类甚至超人级别的性能是可能的。所有这些都是通过观察游戏来完成的。这引起了我的注意。
一个人工智能体被开发出来，它可以直接从电脑上学习玩各种经典的 Atari 2600 电脑游戏
www.nature.com](https://www.nature.com/articles/nature14236) 
One of the Atari 2600 games, Breakout. Trained using a simple RL algorithm. After millions of iterations, the computer agent plays at a super-human level.
我开始在雅达利 2600 游戏上运行实验。虽然令人印象深刻，但《T21》并不是一个复杂的游戏。人们可以通过将你的行动(操纵杆)与结果(得分)联系起来的难度来定义复杂性。当人们需要等待很长时间才能观察到效果时，问题就出现了。
An illustration of the problem with more complex games. LEFT: Breakout (ATARI 2600) is very reactive, we get feedback very quickly; RIGHT: Mario Land (NINTENDO GAMEBOY) does not provide immediate information about consequences of an action, there can be stretches of irrelevant observations between important events
为了让学习更有效率，你可以想象尝试从简单的游戏中转移一些知识。这是目前尚未解决的问题，也是一个热门的研究课题。OpenAI 最近发布的一项挑战试图测量这一点:
OpenAI 是一家非营利性的人工智能研究公司，旨在发现和制定通往安全人工智能的道路。
contest.openai.com](https://contest.openai.com/2018-1/)
拥有迁移学习能力不仅会使训练更快，而且我甚至认为，除非有一些先验知识，否则有些问题是无法解决的。你需要的是**数据效率。**让我们来看看波斯王子:
A short clip from Prince of Persia
- 没有明显的分数
- 如果不执行任何操作，需要 60 分钟才能结束游戏(动画中还剩 58 分钟)。
Random actions
你能尝试和 Atari 2600 论文中完全一样的方法吗？如果你想一想，你随机按键到达终点的可能性有多大？
这激励我通过解决这个问题来为社区做贡献。所以，从本质上说，这是一个先有鸡还是先有蛋的问题——你需要一个更好的算法来传递知识，但这需要很长时间的研究和实验，因为我们没有更有效的算法。
Transfer learning example: Imagine first learning a very simple game, such as the one on the lefthand side, then you preserve concepts such as ‘race’, ‘car’, ‘track’, ‘win’, and learn colors or 3D models. We say that the common concepts are ‘transferrable’ between those games. One could define similarity by the amount of transferrable knowledge between two problems. I.e. Tetris and F1 racing will not be similar.
我决定退而求其次，通过让环境变得更快来避免最初的减速。我的目标是:
- 更快的环境(想象你可以用 1/100 的时间完成波斯王子)同时运行 100，000 个游戏
- 更好的研究环境(关注问题，但不做预处理，有各种游戏)
最初，我认为性能瓶颈可能在某种程度上与模拟器代码的复杂性有关(例如 Stella 的代码库相当大，而且它依赖于 C++抽象，这对于模拟器来说不是最佳选择)。
控制台
Arcade Space Invaders
总的来说，我在几个平台上工作过，从可能是有史以来最早的游戏之一(和 Pong 一起)——街机“太空入侵者”、Atari 2600、NES 和 Gameboy 开始。所有用 c 写的。
我能观察到的最大帧速率大约是 2000-3000 FPS。为了开始看到实验结果，我们需要数百万或数十亿帧，所以差距是巨大的。
My Space Invaders running inside an FPGA — low speed debug mode, the counter in the FPGA shows clock cycles elapsed
我想——如果我们可以对这些环境进行硬件加速会怎么样。例如，最初的太空入侵者运行在 1MHz 的 8080 CPU 上。我能够在 3GHz Xeon 上模拟 40MHz 8080 CPU。还不错，但作为概念验证，一旦放入 FPGA，它就能达到 400 MHz。这意味着单个实例可以达到 24000 FPS，相当于 30GHz 的至强 CPU！我有没有提到一个可以在中间层 FPGA 中安装大约 100 个 8080 CPU 的处理器？也就是 2.4M FPS。
Hardware-accelerated Space Invaders, 100MHz, 1/4 of the full speed
Over one hundred tiny cores inside one Xilinx Kintex 7045 FPGA (bright colors, the blue patch in the middle is shared logic for display)
Irregular Execution Path
此时，你可能会问，GPU 呢？简而言之，你需要 MIMD 式的平行，而不是 SIMD 式的。在学生时代，我花了一些时间在 GPU 上实现蒙特卡罗树搜索(AlphaGo 中使用了 MCTS)。
http://olab . is . s . u-Tokyo . AC . jp/~ kamil . rocki/rocki _ IP dps 11 . pdf
当时我花了无数个小时试图让 GPU 和其他种类的 SIMD 硬件(IBM Cell、Xeon Phi、AVX CPU)高效运行这类代码，但都失败了。几年前，我开始相信，如果我能设计出自己的硬件，来解决类似的问题，那就太好了。
Multiple Instruction Multiple Data (MIMD) Parallellism
雅达利 2600，NES 还是 GAMEBOY？
总之，我用《太空入侵者》、《NES》、《2600》和《游戏男孩》实现了 8080。这里有一些事实和他们各自的优势。
NES Pacman
太空入侵者只是一个热身。我们让它工作，但它只是一个游戏，所以它不是很有用。
Atari 2600 是强化学习研究的事实标准。CPU (MOS 6507)是著名的 6502 的简化版本，它的设计更优雅，比 8080 更高效。我没有选择 2600，只是因为我认为它在游戏和图像方面有一定的局限性。
我实现了 NES(任天堂娱乐系统)以及，它与 2600 共享 CPU。游戏比 2600 上的好多了。NES 和 2600 都受到过于复杂的图形处理管道和许多需要支持的墨盒格式的困扰。
与此同时，我重新发现了任天堂 Gameboy。这就是我一直在寻找的东西。
game boy 为什么这么牛逼？
1049 Classic games + 576 for Gameboy Color
- 总共超过 1000 个游戏(经典+彩色)，范围非常广，都是高质量的，其中一些非常具有挑战性(王子)，游戏可以进行某种程度的分组，并为迁移学习和课程学习的研究分配难度(例如，有俄罗斯方块、赛车游戏、Marios 的变体)。试图解决波斯王子可能需要从其他类似的游戏中转移知识，这些游戏有明确的分数(王子没有！)
Nintendo Gameboy is my transfer learning research platform of choice. In this chart I tried to group the games and plot them according to their difficulty (subjective judgement) and their similarity (concepts such as racing, jumping, shooting, variety of Tetris games; Has anyone ever played HATRIS?).
- Gameboy classic 有一个非常简单的屏幕(160×144 2 位颜色),这使得预处理更简单，我们可以专注于重要的事情。在 2600 的情况下，即使是简单的游戏也有各种颜色。除此之外，Gameboy 有更好的显示物体的方式，所以不需要闪烁或最多连续几帧。
IBM Boot ROM + Super Marioland
- 没有像 NES 或 2600 案件中那样疯狂的记忆地图绘制者。你可以用总共 2-3 个地图绘制者来完成大部分游戏工作。
- 它被证明是非常紧凑的，总的来说，我能够用不到 700 行代码完成 C 仿真器，用大约 500 行代码完成我的 Verilog 实现。
- 它带有与街机版相同的简单版《太空入侵者》
Gameboy version of Space Invaders — Basically no preprocessing needed!
这就是我 1989 年的点阵游戏机，以及在 4K 屏幕上通过 HDMI 运行的 FPGA 版本。
这是我以前的游戏机做不到的:
Hardware accelerated Tetris, this is a realtime screen recording at 1/4 the max speed.
实际有用吗？
是的，它是。到目前为止，我已经在一个简单的环境中测试过了，那里有一个外部策略网络，它与 Gameboy 的工作人员进行交互。更具体地说，我使用了分布式 A3C(优势行动者评论家)算法，我将在另一篇文章中描述这一部分。我的一个同事也将它连接到 FPGA convnet，它工作得相当好。更多将在即将到来的职位。
FPGA-NN communication
Distributed A3C setup
Mario land: Initial State. Pressing buttons at random does not get us far. Top right corner shows remaining time. If we get lucky and Goomba, we finish quickly as we touch Goomba. If not, it takes 400 seconds to ‘lose’.
Mario land: after less than 1 hour of gameplay, Mario has learned to run, jump and even discovers a secret room by going into a pipe.
Pac Man: after about one hour of training it can even finish the entire game once (eats all dots), then the games starts all over again.
结论
我喜欢将未来十年视为 HPC 和 AI 最终走到一起的时期。我希望硬件有一定程度的定制允许，这取决于人工智能算法的选择。
The Next Decade
这里的国标 C 码【https://github.com/krocki/gb
奖金 A —调试
这一部分可能值得单独写一篇文章，但是现在，我太累了。
经常有人问我:最难的部分是什么？一切……都很痛苦。一开始就没有 gameboy 的规范。我们所知道的一切都是逆向工程的结果，我的意思是:运行一些代理任务，如游戏或任何其他片段，并观察它是否执行。这与标准的软件调试非常不同，因为在这里执行代码的硬件被调试。我必须想出一些方法来完成这个过程。我有没有提到当它运行在 100MHz 时很难看到一些东西？哦，没有指纹。
One approach to implementing the CPU is to group instructions into some clusters which do more or less the same thing. This is much easier with 6502. LR35092 has much more random stuff and many edge cases. Here is the chart which I used while I worked on the Gameboy CPU. I adopted a greedy strategy of taking the largest chunk of instructions, implementing it and then crossing it out, then repeat. In this case 1/4 of the instructions are ALU, 1/4 are register loads which can be implemented relatively quickly. On the other side of the spectrum, there are some outliers such as load from HL to SP, signed which need to be handled separately
Debugging: you run a piece of code on the hardware you are debugging, record a log for your implementation and some other reference implementation (here I compared Verilog code, left vs my C emulator, right). Then you run diff on the logs to identify discrepancies (blue). One reason for using some automated way is that in many cases I found problems after millions of cycles of execution, where a single CPU flag caused a snowball effect. I tried many approaches and this seemed to be one of the most effective ones.
You need coffee! Lots of it.
These books are 40 years old — It’s incredible to go through them and see the world of computers through the eyes of the users then. I felt like a visitor from the future.
奖金 B——开放式研究申请
最初我想从内存的角度来看游戏，就像 OpenAI 在这篇文章中描述的那样。
关注的序列对序列模型取得了巨大的成功。他们使得神经网络成为可能…
openai.com](https://openai.com/requests-for-research/#q-learning-on-the-ram-variant-of-atari)
令人惊讶的是，当输入为 RAM 状态时，让 Q-learning 正常工作具有意想不到的挑战性。
该项目可能无法解决。如果事实证明 Q-learning 在 Atari 的 RAM 变体上永远不会成功，那将是令人惊讶的，但它有一些机会将被证明是具有挑战性的。
鉴于 Atari 游戏只使用 128B 的内存，处理这些 128B 的内存而不是整个屏幕帧是非常有吸引力的。我得到了不同的结果，所以我开始深入调查。
虽然我不能证明直接从记忆中学习是不可能的，但我可以证明，记忆捕捉游戏的整个状态的假设是错误的。Atari 2600 CPU (6507)使用 128B 的 RAM，但它也可以访问独立电路中的额外寄存器(电视接口适配器，类似于 GPU)。这些寄存器用于存储和处理关于物体的信息(桨、导弹、球、碰撞)。换句话说，当只考虑 RAM 时，这将不可用。同样，NES 和 Gameboy 有额外的寄存器用于屏幕操作和滚动。RAM 本身并不能反映游戏的完整状态。
只有 8080 以直接方式将数据存储到 VRAM 中，这将允许检索整个状态。在其他情况下,“GPU”寄存器插在 CPU 和屏幕缓冲区之间，在 RAM 之外。
琐事:如果你正在研究 GPU 的历史，这可能是第一个“图形加速器”——8080 有一个外部移位寄存器，使用单个命令移动入侵者，卸载 CPU。
参考
[## 我做了一个 NES 模拟器。以下是我对最初的任天堂的了解。
关于 NES 有趣的技术花絮…
medium.com](https://medium.com/@fogleman/i-made-an-nes-emulator-here-s-what-i-learned-about-the-original-nintendo-2e078c9b28fe) [## 为什么我花了 1.5 个月的时间创建一个 Gameboy 模拟器？托梅克的博客
后来我意识到这并不完全正确——创建一个工作程序是一个相当大的挑战。同时…
博客. rekawek.eu](https://blog.rekawek.eu/2017/02/09/coffee-gb/) [## 深度强化学习:来自像素的 Pong
一位计算机科学家的思考。
karpathy.github.io](http://karpathy.github.io/2016/05/31/rl/)
还有更多…
Gameboy 和 NES 是任天堂的注册商标。
文件结束
一种利用 2D 目标检测进行人脸识别的通用方法
有几种技术可以利用计算机视觉进行面部识别。在这个场景中，我们演示了一种通过面部识别来识别感兴趣的人的方法。有多种方法可以实现积极的身份识别。深度学习可用于映射和识别人脸的三维平面。或者，也可以应用聚焦在人脸上的 2D 对象检测技术。前者需要大规模扫描的 3D 人脸，这是非常昂贵的，需要目前正在研究和开发的技术，并达到远低于 80%的精度水平。用这种方法看到成功将是困难的、不确定的和昂贵的。我们认为这种方法不太实际。后者(2D 物体探测技术)需要高分辨率的照片作为数据集，这可以通过高质量的相机来实现。2D 深度学习技术，卷积神经网络(CNN)今天被应用于各种类型的识别。利用 2D 面部识别技术来识别感兴趣的人是可行的，可以实现高度的准确性(80%以上)，并且减少了总的技术债务。我们建议开发一个 2D 面部识别模型原型，作为一个实用的方法，积极的面部识别感兴趣的人。
一般方法:
数据是关键。最初，您需要评估和绘制数据采集过程和数据结构(例如，使用什么相机，照明如何，照片中通常出现多少人)。基于以上所述，您应该收集速度和准确性要求，建立什么被认为是“可接受的”。然后，评估并记录硬件限制。
一旦正确的评估完成，开始并建立使用系统的 UI/UX 流程，并定义 API 端点。建立视频/照片和标签的数据收集系统。然后将数据集解析成训练、验证和测试桶。需要开发适当的数据预处理管道来减少图像和视频数据中的模糊和噪声。
构建深度学习算法(CNN)或将迁移学习应用于现有模型；接受过特定人脸检测、嵌入和分类模型的培训。在开发了核心人脸检测和识别模型之后，您将需要开发一个跟踪模型来跟踪整个视频中的多个对象。
最后，集成和服务于模型的数据库和后端将被开发。需要定义和实现一个人在回路中的过程来持续改进模型。
型号选择/开发:
我们建议从人脸检测、关键点提取、对齐和嵌入的开源解决方案开始。这允许对感兴趣的应用的现有技术进行快速评估，并通知我们集中改进建模技术或数据采集的方向。然后，收集带标签的图像或视频，并标注人脸边界框和身份。
通过进一步的特征工程和超参数调整来调整所选模型。然后，测试系统并找出系统故障所在。一旦选择了模型，您将需要为重新培训、测试、部署和监控计划创建一个计划。
技术推荐:
对于人脸检测和识别模型，我们建议使用卷积神经网络(例如 ResNet)作为主干。对于人脸检测损失函数，我们推荐使用 SSD 损失(交叉熵和回归损失)。人脸嵌入模型的损失函数有多种选择:
- 交叉熵损失
- 三重损失
- 中心损失
对于关键点提取，我们建议使用回归损失。
当然，你需要训练数据。这可能是一个挑战。为了达到高水平的准确性，您将需要标记的训练数据(1 到 1000 万张脸)。
人脸检测
给定一幅图像，我们需要检测人脸所在的像素区域。下面是一个例子，我们的工作中，面部检测显示。
输入:图像
输出 : 1)包含正面或稍微侧面的面的边界框(上、左、下、右)；2)对应于眼睛、鼻子和嘴的面部上的关键点
采用多路最先进的卷积神经网络进行计算。
面部对齐
使用面部关键点对图像执行 2D 变换，使眼睛和月份处于大致正常的位置。
输入:检测到人脸的裁剪图像块
输出:面部对齐的图像补片
人脸嵌入(描述符)
输入:裁剪后的人脸图像补丁
输出:描述人脸的向量
人脸嵌入模型将是一个卷积神经网络。它可以基于预先训练的模型(例如:最初用于分类 1000 个对象类别-例如猫和狗)。执行迁移学习，并使用交叉熵损失或三重损失(在最近的文献中提出了新的损失)在我们的面部匹配数据集上微调模型。
图上的人脸聚类
嵌入的面部描述符存在于向量空间中。在这个向量空间中两个面的距离表明它们是多么不同。在由面部向量之间的这些成对距离指定的面部图上执行聚类算法。每个聚类包括很可能是同一个人的面部。
输入:人脸描述符
输出:人脸聚类，每个人脸聚类是一组很可能属于同一个人的人脸。
项目流程图:
预处理和标准化
需要带标签的训练数据来提供照片、人脸的带标签的边界框以及每个边界框的身份(人名或 ID)。标签需要组织成一种易于使用的格式，比如 PASCAL VOC 格式。对于原型的范围，我们将假设数据可用性和结构。
部署后改进系统
构建一个人在回路中的组件，该组件审核生产中的模型，并标注边界框和人员身份的正确标签。这包括开发一个批量训练过程，不断吸收新的训练数据。
这个人在回路组件的目的是扩大已知人脸和嵌入向量的数据库。他们会尝试在照片中加入不同的光线，不同的角度，不同的玻璃等等。提高在新照片中认出那个人的回忆率。
准确性取决于 1:N 人脸比较的干扰物比率(每个查询人脸的干扰物数量)。在 mega face(N = 100 万)上，最精确的系统可以达到大约 30%的精度。在 1:1 人脸比对上，可以达到 99%+的准确率。
潜在挑战:
用于培训和验证的图像质量将阻碍项目的成功:
- 光线太弱或太强。
- 低分辨率。
- 运动模糊。
- 鱼眼相机。
复杂的场景提供了额外的噪声，影响了准确识别面部的能力。例如，拥挤的场景和来自其他对象的遮挡必须被隔离。
此外，数据大小可能会成为技术基础架构和资源的负担。数千万张照片产生了海量数据。这对于有限的计算资源是有问题的，并且会产生长的训练周期。这导致了缓慢的迭代。
硬件提供了约束。如果模型部署在“边缘”(例如，在 Raspberry Pi 或移动电话上)而不是云实例上，则只能使用模型的一小部分。大量使用可能会有问题。如果许多相机数据流同时消耗 API，那么对面部识别 API 的请求会淹没计算能力。
会有速度限制。如果需要模型来生成实时预测，这就限制了可以使用的模型的大小和类型。
克服挑战:
图像质量
使用高质量和高分辨率的相机；安装摄像机的方式应能获得最佳的角度、距离和照明。
复杂场景
对模型的局限性设定期望。关注确定最高价值目标的更窄的用例。
数据大小
可能需要更多计算资源。我们将开发一个更好的分布式模型训练算法来减轻资源的负担。
高用量
在模型部署期间提供足够的机器资源。
培训流程:
首先，需要将数据收集到对象存储中的一个桶中(例如 AWS S3)。使用 TensorFlow/PyTorch/Caffe 框架开发培训管道。则需要供应 GPU 资源。一旦模型被训练，模型工件和度量就被持久化到对象存储中。
在模型训练期间，我们建议执行连续的超参数调整。我们调整的 nob 包括但不限于:
- 主干架构:VGG，雷斯网，DenseNet，MobileNet，NASNet。
- 学习率调度。
- 批量大小。
可能需要分布式培训。如果图像数量超过数百万，通常需要在多个 GPU 实例上进行分布式训练。我们推荐 Horovod(来自优步)。
虽然 3D 面部识别是一个可行的解决方案，但它并非没有挑战。它的计算成本很高，当前的精度基准对许多应用来说可能太低。人脸识别的 2D 对象检测应该足以处理大多数用例。以上提供了许多方法中的一种。我们鼓励其他人分享他们的！
功夫。AI 是一家人工智能咨询公司，帮助公司建立他们的战略，运营和部署人工智能解决方案。请点击 查看 www.kungfu.ai
R 中整洁统计的简明指南(第 1 部分)
While data analysis in R can seem intimidating, we will explore how to use it effectively and clearly!
介绍
Jesse Maegan 在 Twitter 上提出了以下问题，并在 R 频谱的所有点上与#RStats 用户展开了一场大讨论。
许多人加入进来，提出了缺少详细插图的包、使用其他编程语言的社会团体的影响、难以找到与你相关的例子、需要更多关于大学编程与统计的指导，以及许多其他主题。
两条推文给我留下了深刻的印象，与我在 r 中的第一个需求非常相似。关于您感兴趣的主题的多元和描述性统计！对于我们许多从事科学研究的人来说，通过 ANOVAs 和线性回归/相关性以及简单的总结(描述性统计)进行的均值测试将解决我们的许多问题。
然后，David 带来了一个很好的观点和视角，他是一个需要指导学生/同事的人，这些学生/同事表达了使用 R 的需求，但在数据科学的背景下没有 R 的需求或愿望。对于局外人来说，R 的编码本质可能是令人生畏的，让你看起来好像是在“矩阵”中执行数据分析
我肯定去过那里，还没有找到一个坚实的例子走过多变量 ANOVAs，后 hocs，绘图，和描述性统计。基本上是一个完整的 R 语言数据分析的例子，但用于学术用途。
所以我告诉大卫我会写一篇博客，我们开始吧！
如果你觉得打开和使用 R Studio 很舒服，请随意跳到第二部分！里面有所有的统计数据和图表！
We’ll get here eventually I promise!
R/R 工作室
我将大量借用切斯特·伊斯梅和艾伯特·金的精彩电子书《T4 现代潜水》。如果你对如何使用 R 感兴趣，他们的电子书是最好的之一，而且还是免费的！
如果你需要安装 R，请参见他们关于下载和安装 R and R 工作室的章节。
现在开始使用 R Studio！打个电话的比方，R 是运行你手机的电脑，而 R Studio 是主屏幕。R Studio 提供了交互的接口，告诉 R 做什么。R 分为 4 个主要部分，左边是 R 脚本(在这里输入代码)和 R 控制台，而右边是环境和文件/图/包部分。
Source: Grolemund & Wickham — R for Data Science
继续进行电话类比，我还想介绍一下包。软件包类似于手机上的应用程序。虽然你的手机/R 可以做很多开箱即用的事情，但应用程序/软件包增加了功能！要在 R 中安装这些“应用”，我们需要先明确告诉 R 安装包！
要开始用 R 编写代码，首先同时输入“ctrl + shift + N”。这将打开一个新的 R 脚本，基本上是一个文件，您可以在其中保存键入的 R 代码和您对正在做的事情的评论。或者，您可以单击左上角的绿色“+”并从下拉菜单中选择 R Script。
在那里，您可以安装我们进行这些分析所需的一个包。首先在新的 R 脚本中键入以下代码，然后在光标停留在同一行的同时单击 ctrl + enter。不管包的名字是什么，我们都需要用" "
install.packages("tidyverse")
这将首次将感兴趣的软件包安装到您的计算机上。您不需要再次安装它，除非您希望在将来下载更新的软件包。
每当你启动 R Studio 时，你也需要加载这个包。这就相当于在手机上打开一个 app！它将加载额外的功能，并允许您通过键入代码来调用它们。
library("tidyverse")
在您的脚本/分析的顶部加载包是一个好主意，因此很清楚需要什么包！我们还可以在笔记上添加注释，记录我们正在做的事情以及我们为什么要这么做！R 中的注释可以在你的注释前加一个#形成。这个#告诉 R 不要“读取”同一行上的任何内容，防止它在试图读取您的精彩评论时向您抛出错误！
# Adding the hashtag/pound sign tells R to not evaluate this line
不涉及太多细节,tidyverse包为我们提供了一个优雅的描述性统计、漂亮的出版级图表的选项，并使我们的代码更具可读性！
我想介绍的下一个概念是赋值操作符
<-
这告诉 R 你正在把右边的所有东西分配给左边的对象。让我们将变量x赋值为 5，y赋值为 6，然后做一些基本的数学运算。
x <- 5
y <- 6
x * y# output
>[1] 30
您可以在 R 中做同样的事情，一旦您运行了x * y代码，您将看到 30 在输出控制台中弹出！恭喜你写出了你的第一个代码！
我们将跳过前面的部分，直接进入分析的开始部分。如果你想了解更多关于 R 的基础知识，我再次建议你去看看Modern Dive——R 的统计介绍。此外，www.datacamp.com有很多与 R 相关的免费内容，而且经常打折！
继续第二部分！
R 中整洁统计的简明指南(第 2 部分)
While data analysis in R can seem intimidating, we will explore how to use it effectively and clearly!
第 1 部分让您踏上用 R 代码运行统计数据的旅程。
介绍
在 Jesse Maegan ( @kiersi )在 Twitter 上发起了一场大讨论之后，我决定发布一些(伪造的)实验治疗数据。这些数据与一种名为 AD-x37 的新(假)研究药物相对应，AD-x37 是一种理论药物，已被证明对阿尔茨海默病小鼠模型的认知衰退具有有益的结果。在目前的实验中，我们将对药物在减少痴呆患者的认知障碍方面是否有效进行统计测试。此处见数据。
我们将使用 MMSE(简易精神状态测验)分数来评估认知障碍的程度。在真正的临床试验中，会记录许多其他变量，但为了一个简单但多变量的例子，我们将只关注 MMSE。
Source: Folstein et al, 1975 — J Psychiatr Res 12:189–198
我们将通过tidyverse加载、绘制、分析和保存我们的分析输出，这是一个为数据分析而设计的“R 包的自以为是的集合”。我们将把依赖性限制在两个包中:tidyverse和broom，而对其余的包使用基数 R。在我看来，这两个包极大地改进了数据分析工作流程。虽然其他侧重于统计的包提供了额外的统计测试，但 base R 有相当好的能力来执行开箱即用的统计分析。我将使用knitr::kable为 markdown 文档生成一些 html 表格，但是对于工作流来说这不是必需的。
此外，我将上传本例中使用的 excel 表，以便您可以自己重新创建工作流程。您可以简单地复制粘贴此处看到的代码，它将在 R 中运行。如果您希望在 R-Markdown 文档中看到整个工作流，请参见此处的。R Markdown 是在 R 内部创建的一个文档，它允许您编写代码、内联执行代码，并在编写过程中编写注释/笔记。你可以把它想象成能够在一个基本的 Word 文档中编写 R 代码(但是它能做的远不止这些！).
虽然您可能对我提供的数据集不感兴趣，但这有望为您提供一个清晰的工作流程，让您可以交换感兴趣的数据并完成基本分析！
装载扫把、扫帚和编织机
使用库函数，我们将加载tidyverse。如果您以前从未安装过它，您也可以使用install.packages("tidyverse")调用来首次安装它。这个包包括ggplot2(图形)、dplyr / tidyr(汇总统计、数据操作)、和readxl(读取 excel 文件)以及管道%>%，这将使我们的代码更具可读性！我们还将加载broom包来整理我们的一些统计输出。最后，我们将通过knitr::kable,加载knitr来制作漂亮的 html 表格，但不需要简单地将输出保存到 Excel。
# Load libraries
library(tidyverse)
library(broom)
library(knitr)
library(readxl)
这将输出一些关于正在加载的包和任何函数调用冲突的消息。
加载数据
当我呼叫readxl::read_xlsx时，您也可以简单地使用read_xlsx, ，但是为了透明起见，我将使用完整的呼叫开始。使用::调用函数的概念很重要，因为一些包在函数上有冲突，例如多个包包括函数select 和summarize。这样，我们可以明确我们希望 R 从哪个包中调用我们的函数，所以package::function！要了解更多关于调用函数时“名称空间”的概念，请查看这里。
readxl 很不幸是一个有趣的例子，因为安装tidyverse 会安装readxl，但是通过library 调用加载tidyverse 时readxl 不会被加载。因此，我们必须像加载任何其他包一样加载readxl ，或者像readxl::read_xlsx. readxl 中允许我们读取的那样调用包和名称。xls，。或者，您可以将您的 Excel 表转换为。csv，可以通过read_csv()读取。通过使用来自dplyr的glimpse函数，我们可以看到变量是如何导入的，以及前几行。
# Read excel file
raw_df <- readxl::read_xlsx("ad_treatment.xlsx")
dplyr::glimpse(raw_df)# this is the output from glimpse(raw_df)
Observations: 600
Variables: 5
$ age <dbl> 80, 85, 82, 80, 83, 79, 82, 79, 80, 79, ...
$ sex <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, ...
$ health_status <chr> "Healthy", "Healthy", "Healthy", ...
$ drug_treatment <chr> "Placebo", "Placebo", "Placebo", ...
$ mmse <dbl> 24.78988, 24.88192, 25.10903, 23.38...
我们现在可以收集一些关于数据集的信息。也就是说，我们有 3 个分类/因素变量:性别、健康状况和药物治疗以及 1 个因变量(DV): mmse。我们也有年龄，但重要的是它被记录为一个离散的数字，而不是一个因素(如 85 岁，而不是老)。因此，我们可以考虑年龄，但我们不会把它作为方差分析的一个因素。
检查数据分布
我们将使用第一个ggplot2调用来创建一个显示年龄分布的图表。为了分解我们正在做的事情，我们需要调用 ggplot，告诉它使用什么数据，并使用aes 或“美学”调用来分配 x 坐标。然后我们添加一个+ ，告诉 ggplot 包含下一行代码。geom_density告诉 R 我们想要创建一个密度分布层，我们想要fill用蓝色！关于ggplot2的更多信息，请点击此处或此处。
ggplot(data = raw_df, aes(x = age)) +
geom_density(fill = "blue")
图表显示，年龄实际上只在 79-85 岁之间，而且实际上不存在年龄过高或过低的情况。我们可以通过一个dplyrsummary 调用或通过调用 base R 中的range来确认年龄范围 s。顺便提一下，我们现在可以讨论使用管道或%>%.管道将结果或数据从它的左边传递到右边。有关该管道的更多信息，请参见这里。
我们可以把下面的代码理解为取raw_df，然后通过取age变量的min和max来summarize它。现在，因为我们从raw_df开始，R 知道我们想要从这个数据帧中取出列age。
raw_df %>% summarize(
min = min(age),
max = max(age))# A tibble: 1 x 2
min max
<dbl> <dbl>
1 79.0 85.0
或者，我们可以使用 base R range函数，这需要使用$。美元符号表示 R 应该使用来自raw_df的age列。这两个函数给我们的结果是一样的，最小数和最大数。
range(raw_df$age)[1] 79 85
有关使用这两种语法的更多信息，请点击这里或者查看备忘单这里。
实验变量水平如何？
现在，虽然我非常清楚这个数据框架中的变量，但你可能不会不探索它！为了快速确定drug_treatment组、health_status组以及它们如何交互，我们可以做一个table调用。通过在drug_treatment和health_status上调用它，我们得到了一个很好的表格，它分解了每个变量组中有多少行。
table(raw_df$drug_treatment, raw_df$health_status)#output below Alzheimer's Healthy
High Dose 100 100
Low dose 100 100
Placebo 100 100
或者，我们可以用下面的代码在dplyr中做同样的事情。
raw_df %>%
group_by(drug_treatment, health_status) %>%
count()
现在我们知道了我们感兴趣的变量的水平，并且每个整体治疗组有 100 名患者！
因变量的数据探索
在运行我们的汇总统计数据之前，我们可以通过一个geom_boxplot调用来可视化范围、集中趋势和四分位数。
ggplot(data = raw_df, # add the data
aes(x = drug_treatment, y = mmse, # set x, y coordinates
color = drug_treatment)) + # color by treatment
geom_boxplot() +
facet_grid(~health_status) # create panes base on health status
我们将健康和阿尔茨海默病患者的数据分成单独的图形面(或窗格),并按药物治疗分成每个面内的组。这个图表告诉我们一些有趣的事情。看起来我们的(假的)神奇药物确实有效果了！让我们用描述性统计来探讨这个问题。
虽然这是一个探索性的图表，我们不一定要将其“调整”到完美，但我们可以注意到，我们的药物治疗应该是安慰剂< Low dose < High Dose and we should have Healthy patients presented first, and Alzheimer’s patients second. This is something we can fix in our next section!
汇总统计
我们希望生成 mmse 分数的平均值和标准误差，这对于测量集中趋势和创建我们的最终出版图表是有用的。我们有性别、药物治疗和健康状况的分类变量。然而，回到我们之前的glimpse调用，我们可以看到数据没有被正确“编码”。即性是一个dbl(数字)，没有描述性名称，health_status / drug_treatment是chr(人物)！这些都需要换算成因子！
Observations: 600
Variables: 5
$ age <dbl> 80, 85, 82, 80, 83, 79, 82, 79, 80, 79, ...
$ sex <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1...
$ health_status <chr> "Healthy", "Healthy", "Healthy", ...
$ drug_treatment <chr> "Placebo", "Placebo", "Placebo", ...
$ mmse <dbl> 24.78988, 24.88192, 25.10903, 23.38...
我们可以使用dplyr::mutate函数告诉 R 我们想要改变(突变)感兴趣的变量中的行。因此，我们将获取sex、drug_treatment和health_status 列中的数据，并将它们从数字或字符转换成因子变量！dplyr::mutate还可以表演数学，以及许多其他有趣的事情。更多信息，请参见此处。
我们将使用mutate 函数和基数 R factor函数将我们的变量转换成适当的因子，并给它们加上标签(表示性别)或重新排列因子的级别。
我们需要非常小心地准确键入列中显示的标签，否则它会用 NA 替换那些拼写错误的标签。例如，您是否注意到High Dose有一个大写的“D ”,而Low dose有一个小写的“D ”?
sum_df <- raw_df %>%
mutate(
sex = factor(sex,
labels = c("Male", "Females")),
drug_treatment = factor(drug_treatment,
levels = c("Placebo", "Low dose", "High Dose")),
health_status = factor(health_status,
levels = c("Healthy", "Alzheimer's"))
)
尽管 R 非常强大，但它需要明确而准确的代码输入来完成最终目标。因此，如果我们键入“High d ose ”,它会给出安娜，而“High D ose”会正确输出。我们现在将age 和mmse 视为dbl (数字)，将sex、health_status和drug_treatment 视为因子。
Observations: 600
Variables: 5
$ age <dbl> 80, 85, 82, 80, 83, 79, 82, 79, 80, 79,...
$ sex <fct> Male, Male, Male, Male, Male, Females, Mal...
$ health_status <fct> Healthy, Healthy, Healthy, Healthy, Health...
$ drug_treatment <fct> Placebo, Placebo, Placebo, Placebo, Placeb...
$ mmse <dbl> 24.78988, 24.88192, 25.10903, 24.92636,...
现在一切都已正确编码，我们可以计算我们的平均值和标准误差(se =标准偏差/样本数的平方根)！我们将使用group_by来告诉 R 我们想要…分组的因素！然后，我们将通过首先调用summarize然后用mmse_mean和mmse_se以及样本数量n()指定我们想要的汇总来创建命名汇总。最后，我们将ungroup，从数据帧中删除group_by代码。
sum_df <- sum_df %>%
group_by(sex, health_status, drug_treatment) %>%
summarize(mmse_mean = mean(mmse),
mmse_se = sd(mmse)/sqrt(n()),
n_samples = n()) %>%
ungroup() # ungrouping variable is a good habit to prevent errors
现在我们有了一个格式良好的数据框架，可以保存到 Excel 中，或用于绘图。我们需要指出我们正在写什么数据(sum_df)以及我们希望结果文件被命名为什么(“adx37_sum_stats.csv”)。
# code to save the table into a .csv Excel file
write.csv(sum_df, "adx37_sum_stats.csv")
总结图
通过调用 ggplot 函数，我们可以生成一个初步的摘要图。
ggplot(data = sum_df, # add the data
aes(x = drug_treatment, #set x, y coordinates
y = mmse_mean,
group = drug_treatment, # group by treatment
color = drug_treatment)) + # color by treatment
geom_point(size = 3) + # set size of the dots
facet_grid(sex~health_status) # create facets by sex and status
我们现在可以看到，该图按照药物治疗和健康状况进行了适当的分类。我们在最终的图表上还有一些工作要做，但是让我们先来看看 ANOVAs！
方差分析终于出来了！
我们将准备一个数据框架，通过方差分析进行分析。我们需要再次确保通过mutate将我们的因素作为因素，并且顺序正确。这是 ANOVA/事后检验工作所必需的，并使事后 hoc 和 ANOVA 输出更容易读取。
stats_df <- raw_df %>% # start with data
mutate(drug_treatment = factor(drug_treatment,
levels = c("Placebo", "Low dose", "High Dose")),
sex = factor(sex,
labels = c("Male", "Female")),
health_status = factor(health_status,
levels = c("Healthy", "Alzheimer's")))glimpse(stats_df)#output belowObservations: 600
Variables: 5
$ age <dbl> 80, 85, 82, 80, 83, 79, 82, 79, 80, 79...
$ sex <fct> Male, Male, Male, Male, Male, Male, ...
$ health_status <fct> Healthy, Healthy, Healthy, Healthy...
$ drug_treatment <fct> Placebo, Placebo, Placebo, Placebo,...
$ mmse <dbl> 24.78988, 24.88192, 25.10903...
这使我们的数据框架进入工作状态！
通过aov功能调用 ANOVA。基本语法通过下面的伪代码显示。我们首先放入因变量(在我们的例子中是 mmse ),然后是一个~,然后是我们想要测试的自变量。最后，我们指定使用什么数据。
aov(dependent_variable ~ independent variable, data = data_df)
我们可以通过下面的代码添加真实的数据集:
# this gives main effects AND interactions
ad_aov <- aov(mmse ~ sex * drug_treatment * health_status,
data = stats_df)# this would give ONLY main effects
ad_aov <- aov(mmse ~ sex + drug_treatment + health_status, data = stats_df)
因为我们有 3 个独立变量，所以我们要做出选择。我们可以简单地通过在每个变量之间添加一个+来寻找主效应，或者我们可以通过在每个变量之间添加一个*来寻找主效应和交互作用。确保不要用逗号替换+或*，因为这会导致错误！
# this throws an error because we shouldn't use commas in between!
ad_aov <- aov(mmse ~ sex, drug_treatment, health_status, data = stats_df)
通过将 ANOVA 分配给ad_aov对象，我们可以调用它的summary 来查看 ANOVA 的结果。
# look at effects and interactions
summary(ad_aov) Df Sum Sq Mean Sq F value Pr(>F)
sex 1 0 0 0.047 0.828
drug_treatment 2 3601 1801 909.213 **<2e-16**
health_status 1 10789 10789 5447.953 **<2e-16**
sex:drug_treatment 2 8 4 2.070 0.127
sex:health_status 1 5 5 2.448 0.118
drug_treatment:health_status 2 2842 1421 717.584 **<2e-16**
sex:drug_treatment:health_status 2 5 2 1.213 0.298
Residuals 588 1164 2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
总结给出了自由度、平方和、均方、F 值和 p 值。我在<2e -16, these p values are so small that R switches to scientific notation. So we see significant main effects of drug treatment, health status, and an interaction of drug treatment by health status. We can interpret that Alzheimer’s patients had different cognitive scores than healthy, and that drug treatment had an effect on cognitive scores. Importantly, sex was not a significant factor, as p = 0.828. Variables being scored as significant or non-significant can both be important!
We can also use 【 to clean up the results of the ANOVA and put them into a dataframe. This is useful for storage, or for automation of some analysis for future ANOVAs.
# this extracts ANOVA output into a nice tidy dataframe
tidy_ad_aov <- tidy(ad_aov)# which we can save to Excel
write.csv(tidy_ad_aov, "ad_aov.csv")
However, we don’t know the direction of changes, or where the changes occurred. Was it just the high dose? Low dose? Both? We need follow-up post hoc tests to determine these answers!
**后 hocs >博士后(学术界笑话<爸爸笑话)**上加了一个醒目的强调
我们有多种方式来看待后 hocs。我将在这一部分展示两个。
对于成对的，我们需要使用$从每个数据帧中选择列，并通过:查看交互。我们的第一个两两比较没有对多重比较进行校正，与无保护的 Fisher’s-LSD post-hoc 相当。这一点也不严格，考虑到我们所做的大量比较，建议要么继续进行 p.adjusting Bonferonni 校正(将p.adj = “none”改为p.adj = “bonf”)要么进行 Tukey 事后测试，如下例所示。您可以看到，由于dataset$column方法和每次交互之间对:的需要，这个方法读起来有点混乱。我们可以这样理解，我们需要对sex 和drug_treatment 和health_status的相互作用进行成对 t.test，它给出了这些因素相对于其他因素的每一次迭代。
# call and save the pair.t.test
ad_pairwise <- pairwise.t.test(stats_df$mmse, stats_df$sex**:**stats_df$drug_treatment**:**stats_df$health_status,
p.adj = "none")
此外，我们需要提取 p 值矩阵并保存到 Excel 文件中以备将来使用。
我们通过简单地用broom::tidy包装我们的ad_last posthoc 来做到这一点。
# tidy the post hoc
tidy_ad_pairwise <- broom::tidy(ad_pairwise)# save to excel
write.csv(tidy_ad_pairwise, "tidy_ad_pairwise.csv")
Tukey post-hoc 比未调整的成对 t-test 更容易调用。请注意，我们已经将 Tukey 结果封装在broom::tidy中，以保存为一个整洁的数据帧！TukeyHSD 调用合并了 ANOVA 调用的结果，并且优于前面的方法。
下面的代码可以被阅读，因为我们想要对我们的ad_aov ANOVA 的结果进行 Tukey 事后测试，该 ANOVA 跨越drug_treatment by health_status的sex 的交互。注意每个变量之间的‘sex:drug_treatment:health_status’和:的引号。这些是告诉 R 我们希望 Tukey 如何运行所必需的！一旦完成，R 就对它运行 tidy，使它成为一个好的数据帧，类似于我们之前的成对测试。然后我们可以将结果保存到 Excel 中！
# call and tidy the tukey posthoc
tidy_ad_tukey <- tidy(TukeyHSD(ad_aov, which = 'sex:drug_treatment:health_status'))# save to excel
write.csv(tidy_tukey_ad, "tukey_ad.csv")
出版图
现在，我们已经生成了 ANOVAs 和 post-hoc，并将其保存到 Excel 中进行存储，我们可以开始制作出版物级别的图表了！
图表允许极端的定制，我在图表上添加的一些内容是个人的选择，因此我建议与你所在领域的导师或有经验的成员讨论。条形图在我的领域中无处不在，虽然我认为绘制成箱线图可以告诉我们更多的数据，但我首先会从条形图开始。
我们的目标是绘制平均值、标准误差，并指出其发生的意义。我将使用tribble 函数手工制作一个定制的数据帧，而不是依赖一个包来标记重要性。除了这样做还有其他选择，但是我可以很容易地控制使用这种方法会发生什么，而且数据帧包含的内容非常明显。下面的例子显示了三层结构的基本原理。我们用~指定列，然后显式地写出我们在列的每一行中想要的内容。
[tribble](https://rdrr.io/cran/tibble/man/tribble.html)(
~colA, ~colB,
"a", 1,
"b", 2,
"c", 3
)# Output below
*# A tibble: 3 x 2*
colA colB
<chr> <dbl>
1 a 1
2 b 2
3 [c](https://rdrr.io/r/base/c.html) 3
这是我们制作自定义数据帧的实际代码。
# make the dataframe with specific points of interest to add *
sig_df <- tribble(
~drug_treatment, ~ health_status, ~sex, ~mmse_mean,
"Low dose", "Alzheimer's", "Male", 17,
"High Dose", "Alzheimer's", "Male", 25,
"Low dose", "Alzheimer's", "Female", 18,
"High Dose", "Alzheimer's", "Female", 24
)# convert the variables to factors again :)
sig_df <- sig_df %>%
mutate(drug_treatment = factor(drug_treatment,
levels = c("Placebo", "Low dose", "High Dose")),
sex = factor(sex,
levels = c("Male", "Female")),
health_status = factor(health_status,
levels = c("Healthy", "Alzheimer's")))# Output below
# A tibble: 4 x 4
drug_treatment health_status sex mmse_mean
<fctr> <fctr> <fctr> <dbl>
1 Low dose Alzheimer's Male 17.0
2 High Dose Alzheimer's Male 25.0
3 Low dose Alzheimer's Female 18.0
4 High Dose Alzheimer's Female 24.0
现在我们有了这个数据框，我们可以在一个geom_text调用中使用它，用重要性标签来标记我们的条，如*所示。
下面是最终发布图在ggplot2代码中的样子。你会注意到我把它赋给了g1，而不是直接调用它。这意味着我将不得不调用g1 来查看图表，但我现在可以保存它！为了理解我们正在做的事情，我像以前一样调用初始的ggplot 调用，但是添加了一个误差条层、一个条形图层、分为性别和健康状态的窗格、切换到另一个外观(theme_bw)、手动设置颜色、通过主题进行微小调整、添加* 以指示重要性，最后在添加图形标题的同时更改轴标签。
g1 <- ggplot(data = sum_df,
aes(x = drug_treatment, y = mmse_mean, fill = drug_treatment,
group = drug_treatment)) +
geom_errorbar(aes(ymin = mmse_mean - mmse_se,
ymax = mmse_mean + mmse_se), width = 0.5) +
geom_bar(color = "black", stat = "identity", width = 0.7) +
facet_grid(sex~health_status) +
theme_bw() +
scale_fill_manual(values = c("white", "grey", "black")) +
theme(legend.position = "NULL",
legend.title = element_blank(),
axis.title = element_text(size = 20),
legend.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 12)) +
geom_text(data = sig_df, label = "*", size = 8) +
labs(x = "\nDrug Treatment",
y = "Cognitive Function (MMSE)\n",
caption = "\nFigure 1\. Effect of novel drug treatment AD-x37 on cognitive function in healthy and demented elderly adults. \nn = 100/treatment group (total n = 600), * indicates significance at p < 0.001")g1# save the graph!
ggsave("ad_publication_graph.png", g1, height = 7, width = 8, units = "in")
保存是通过ggsave 函数完成的，在这里我们需要用" "，告诉 R 我们想要哪个 ggplot 对象(g1)，并通过高度、宽度和单位指示大小。
还有最后的图！
我认为说你可以通过简单地重现我的例子来学习ggplot是一种伤害。因此，我想为您指出数据科学教科书的 R 方向，以及现代潜水电子书。这些免费的电子书有大量的信息，可能超出你今天需要完成的，但会在你未来的努力中为你服务。他们关于数据可视化的章节对 R 绘图入门非常有帮助！
谢谢
如果你已经做到了这一步，对你有好处！我希望这有所帮助，如果你有任何问题，我建议你通过#rstats 标签联系 twitter，或者你可以在 Twitter 上找到我@thomas_mock。
此外， Jesse Maegan 有一个 R for Data Science Slack 频道，你可以在那里学习和提问。r 工作室(潮汐看护者)在 https://community.rstudio.com/举办他们自己的论坛。
购物篮分析——关联规则简介
Source: UofT
介绍
购物篮分析是大型零售商用来发现商品之间关联的关键技术之一。它通过寻找在交易中频繁出现的项目组合来工作。换句话说，它允许零售商识别人们购买的商品之间的关系。
关联规则被广泛用于分析零售购物篮或交易数据，并且旨在基于强规则的概念，使用感兴趣的度量来识别在交易数据中发现的强规则。
关联规则的一个例子
- 假设有 100 个客户
- 其中 10 人买了牛奶，8 人买了黄油，6 人两样都买了。
- 买了牛奶= >买了黄油
- 支持= P(牛奶和黄油)= 6/100 = 0.06
- 信心=支持度/P(黄油)= 0.06/0.08 = 0.75
- lift =置信度/P(牛奶)= 0.75/0.10 = 7.5
注意:这个例子非常小。在实践中，一个规则需要数百个事务的支持，才能被认为具有统计意义，数据集通常包含数千或数百万个事务。
好了，理论到此为止，让我们来看看代码。
我们今天使用的数据集来自 UCI 机器学习库。该数据集被称为“在线零售”，可以在这里找到。它包含一家总部位于英国的注册在线零售商在 2010 年 1 月 12 日至 2011 年 9 月 12 日之间发生的所有交易。
装载包裹
library(tidyverse)
library(readxl)
library(knitr)
library(ggplot2)
library(lubridate)
library(arules)
library(arulesViz)
library(plyr)
数据预处理和探索
retail <- read_excel('Online_retail.xlsx')
retail <- retail[complete.cases(retail), ]
retail <- retail %>% mutate(Description = as.factor(Description))
retail <- retail %>% mutate(Country = as.factor(Country))
retail$Date <- as.Date(retail$InvoiceDate)
retail$Time <- format(retail$InvoiceDate,"%H:%M:%S")
retail$InvoiceNo <- as.numeric(as.character(retail$InvoiceNo))glimpse(retail)
经过预处理后，数据集包括 406，829 条记录和 10 个字段:发票号、库存代码、描述、数量、发票日期、单价、客户 ID、国家、日期、时间。
人们通常在什么时候上网购物？
为了找到这个问题的答案，我们需要从时间列中提取“小时”。
retail$Time <- as.factor(retail$Time)
a <- hms(as.character(retail$Time))
retail$Time = hour(a)retail %>%
ggplot(aes(x=Time)) +
geom_histogram(stat="count",fill="indianred")
Figure 1. Shopping time distribution
一天中的小时和订单量之间有明显的偏差。大多数订单发生在 10:00-15:00 之间。
每位顾客购买多少件商品？
detach("package:plyr", unload=TRUE)retail %>%
group_by(InvoiceNo) %>%
summarize(n_items = mean(Quantity)) %>%
ggplot(aes(x=n_items))+
geom_histogram(fill="indianred", bins = 100000) +
geom_rug()+
coord_cartesian(xlim=c(0,80))
Figure 2. Number of items per invoice distribution
大多数人购买的商品少于 10 件(每张发票少于 10 件)。
十大畅销书
tmp <- retail %>%
group_by(StockCode, Description) %>%
summarize(count = n()) %>%
arrange(desc(count))
tmp <- head(tmp, n=10)
tmptmp %>%
ggplot(aes(x=reorder(Description,count), y=count))+
geom_bar(stat="identity",fill="indian red")+
coord_flip()
Figure 3. Top 10 best sellers
在线零售商的关联规则
在使用任何规则挖掘算法之前，我们需要将数据从数据帧格式转换为事务，这样我们就可以将所有购买的商品放在一行中。例如，这是我们需要的格式:
Source: Microsoft
retail_sorted <- retail[order(retail$CustomerID),]
library(plyr)
itemList <- ddply(retail,c("CustomerID","Date"),
function(df1)paste(df1$Description,
collapse = ","))
函数 ddply()接受一个数据帧，根据一个或多个因素将其分割成多个部分，对这些部分进行计算，然后将结果作为数据帧返回。我们用“，”来分隔不同的项目。
我们只需要项目交易，所以删除 customerID 和 Date 列。
itemList$CustomerID <- NULL
itemList$Date <- NULL
colnames(itemList) <- c("items")
将数据帧写入 csv 文件，并检查我们的事务格式是否正确。
write.csv(itemList,"market_basket.csv", quote = FALSE, row.names = TRUE)
完美！现在我们有了交易数据集，它显示了一起购买的商品的矩阵。我们实际上看不到他们多久一起买一次，也看不到规则。但是我们会找到答案的。
让我们仔细看看我们有多少事务，它们是什么。
tr <- read.transactions('market_basket.csv', format = 'basket', sep=',')
tr
summary(tr)
我们看到 19，296 个事务，这也是行数。有 7，881 个项目—记住项目是我们原始数据集中的产品描述。这里的事务是这 7，881 个项目的集合或子集。
总结给了我们一些有用的信息:
- 密度:稀疏矩阵中非空单元的百分比。换句话说，购买的商品总数除以矩阵中可能的商品总数。我们可以这样使用密度来计算购买了多少物品:19296 X 7881 X 0.0022
- 最常见的项目应该与图 3 中的结果相同。
- 看看交易的规模:2247 笔交易仅针对 1 件商品，1147 笔交易针对 2 件商品，一直到最大的交易:1 笔交易针对 420 件商品。这表明大多数顾客在每次交易中只购买少量商品。
- 数据的分布是右偏的。
让我们看一下项目频率图，它应该与图 3 一致。
itemFrequencyPlot(tr, topN=20, type='absolute')
Figure 4. A bar plot of the support of the 20 most frequent items bought.
创建一些规则
- 我们使用 Arules 库中的 Apriori 算法来挖掘频繁项集和关联规则。该算法对频繁项集采用逐级搜索。
- 我们通过 supp=0.001 和 conf=0.8 来返回支持度至少为 0.1%、置信度至少为 80%的所有规则。
- 我们按照置信度递减的方式对规则进行排序。
- 看一下规则的总结。
rules <- apriori(tr, parameter = list(supp=0.001, conf=0.8))
rules <- sort(rules, by='confidence', decreasing = TRUE)
summary(rules)
这些规则的总结给了我们一些非常有趣的信息:
- 规则数:89697。
- 按长度分布规则:长度为 6 项的规则最多。
- 质量度量的概要:支持、信心和提升的范围。
- 关于数据挖掘的信息:挖掘的数据总量，以及我们之前设置的最小参数。
我们有 89697 条规则。我不想把它们都印出来，所以我们来检查一下前 10 名。
inspect(rules[1:10])
解释相当直接:
- 100%买了“晃晃悠悠的鸡”的顾客也买了“装修”。
- 买了“红茶”的顾客 100%也买了“糖罐”。
并列出这 10 条最重要的规则。
topRules <- rules[1:10]
plot(topRules)
plot(topRules, method="graph")
plot(topRules, method = "grouped")
摘要
在这篇文章中，我们学习了如何在 R 中执行购物篮分析，以及如何解释结果。如果你想用 Python 实现它们， Mlxtend 是一个 Python 库，它实现了这类应用程序的 Apriori 算法。这里可以找到介绍教程。
如果你想用 R Markdown 文件来制作这篇博文，你可以在这里找到。
参考: R 和数据挖掘
用 Go 对遗传算法的简单介绍
通过进化蒙娜丽莎引入遗传算法
乍一看可能不太明显，但计算机科学算法经常受到自然和生物过程的启发。这些算法包括神经网络、粒子群优化、人工蜂群、蚁群优化、进化算法等等。事实上，你可以把生物过程看作是大自然想出的解决问题的简单算法。从这个角度来看，很容易理解为什么这些算法是优化试探法和元试探法。毕竟，大自然为了生存而优化。
试探法，如果你不熟悉这个术语，是通过做一些假设来试图更快解决问题的算法。因此，试探法通常不是最佳的，但在获得最佳结果需要很长时间的情况下更有用。元启发式把这带到下一个层次——它们是一种产生或发现启发式的启发式方法。
遗传算法
遗传算法是基于自然选择过程的元启发式算法。遗传算法是一种进化算法。
自然选择，作为一种更新，是进化中的一个关键机制。这是一个自然过程，随着时间的推移，导致(生物)种群适应它们的环境。这些群体在特征上有所不同。具有更合适特征的个体生物在环境中生存的机会更高。从这些幸存下来的生物体中繁殖的下一代将继承它们的特征，最终产生具有这些更合适特征的种群。
然而，如果整个种群都有相同的特征，而环境改变了，种群就会灭绝。幸运的是，偶尔发生的突变会导致性状的变异，这使得具有更适合变化的环境的性状的生物体得以生存并成为主导。
一个普遍使用的例子是英国胡椒蛾的颜色变化。在 19 世纪早期之前，英格兰的胡椒蛾大多是白色的，它的颜色有助于它躲避食肉鸟类，因为它与浅色的地衣和英国树木很好地融合在一起。然而，在工业革命期间，浅色的地衣因污染而死亡，许多飞蛾栖息的树木被煤烟熏黑。这使得深色的蛾子在躲避捕食者时具有优势，而浅色的蛾子很容易被发现。到了 19 世纪中期，深色蛾子的数量增加了，到了 19 世纪末，几乎所有的胡椒蛾都是深色品种。1956 年《清洁空气法案》的影响打破了这种平衡，深色的蛾子又变得罕见了。
Peppered moth (credits: https://commons.wikimedia.org/wiki/File:Lichte_en_zwarte_versie_berkenspanner_crop.jpg)
这就是自然选择。遗传算法是如何出现的？遗传算法是启发式的，使用与自然选择相同的机制——DNA、种群、变异、适应度、选择、繁殖、遗传和突变。
- DNA——定义拥有一个或多个 DNA 的生物体
- 种群——从 DNA 基因(值)不同的生物体初始种群开始
- 适合度 —确定每种生物对其环境的适合度
- 选择——选择最适合的生物，给它们更高的繁殖机会
- 繁殖 —从选择的最适合的生物体中创造下一代种群
- 遗传 —人口的下一代必须继承基因的值
- 突变——每一代人的基因值都有微小的变化
猴子、打字机和莎士比亚
无限猴子定理说的是，无限数量的猴子坐在无限数量的打字机前随机敲击键盘，只要给足够的时间，最终会重现莎士比亚的全集。假设我们只想让一只猴子复制这句话:*“生存还是毁灭”。*你觉得猴子随机锤出这个要多长时间？
Monkeying around with Shakespeare (credits: https://commons.wikimedia.org/wiki/File:Chimpanzee_seated_at_typewriter.jpg)
报价有 18 个字符(包括空格)。猴子打出*【t】(姑且说全是小盘也可以)的概率是 26 分之一。因此，键入确切序列“生存还是毁灭”*的概率是 26 的 18 次方分之一，或者是大约 29，479，510，200，013，920，000，000，000 的 1 分之一。假设猴子每秒钟打一封信，那么在 934，789，136，225，707，600 年中，它只有 1 次机会打出这句话。这是 934 万亿年的 1 倍。
显然，这种野蛮的方式不会让我们有任何进展。如果我们尝试“改进”报价会怎么样？让我们看看如何使用遗传算法来做到这一点。下面是用于解决该问题的遗传算法的步骤:
定义由一个或多个 DNA 组成的有机体
我们的莎士比亚-喷算法中的一个生物体由单个 DNA 组成，它是一个字节数组和一个代表生物体适应度的数字。
type Organism struct {
DNA []byte
Fitness float64
}
从最初的生物群体开始
我们需要为我们的初始种群创建有机体，所以这里有一个函数来做这件事。
func createOrganism(target []byte) (organism Organism) {
ba := make([]byte, len(target))
for i := 0; i < len(target); i++ {
ba[i] = byte(rand.Intn(95) + 32)
}
organism = Organism{
DNA: ba,
Fitness: 0,
}
organism.calcFitness(target)
return
}
target是我们想要实现的，在这种情况下，它是字符串“生存还是毁灭”的字节数组表示。在这个函数中，我们随机创建一个与目标长度相同的字节数组，并将其设置为新创建的生物体中基因的值。
既然我们可以创造生物体，我们需要创造一个生物体群体。
func createPopulation(target []byte) (population []Organism) {
population = make([]Organism, PopSize)
for i := 0; i < PopSize; i++ {
population[i] = createOrganism(target)
}
return
}
population是一组生物体，而PopSize是一个定义种群规模的全局变量。
找到生物体的适合度
我们需要计算我们种群中生物体的适合度。这在我们创造有机体时被称为早期，但在我们交叉有机体时也会被称为晚期。
func (d *Organism) calcFitness(target []byte) {
score := 0
for i := 0; i < len(d.DNA); i++ {
if d.DNA[i] == target[i] {
score++
}
}
d.Fitness = float64(score) / float64(len(d.DNA))
return
}
这个健身功能比较简单。我们简单地计算基因中的字节与目标匹配的次数。分数除以目标中的总字节数，以使适合度为一个百分比，即 0.0 到 1.0 之间的数字。这意味着如果适应度是 1.0，我们将进化生物体的基因以匹配“生存还是毁灭”的引用。
选择最适合的生物，给它们更高的繁殖机会
现在我们有了一个种群，我们可以找出哪些生物最适合，我们想挑选最适合的生物，让它们繁殖来创造种群的下一代。有许多不同的方法可以做到这一点，但在这种情况下，我们使用的是一种“繁殖池”机制。
func createPool(population []Organism, target []byte, maxFitness float64) (pool []Organism) {
pool = make([]Organism, 0)
// create a pool for next generation
for i := 0; i < len(population); i++ {
population[i].calcFitness(target)
num := int((population[i].Fitness / maxFitness) * 100)
for n := 0; n < num; n++ {
pool = append(pool, population[i])
}
}
return
}
我们所做的是创造一种繁殖池，在这里我根据生物的适合度将相同生物的多个副本放入池中。有机体的适应度越高，池中就有越多的有机体副本。
从所选择的最适合的生物体中创建种群的下一代
之后，我们从繁殖池中随机挑选 2 个生物体，用它们作为亲本，为种群创造下一代生物体。
func naturalSelection(pool []Organism, population []Organism, target []byte) []Organism {
next := make([]Organism, len(population)) for i := 0; i < len(population); i++ {
r1, r2 := rand.Intn(len(pool)), rand.Intn(len(pool))
a := pool[r1]
b := pool[r2] child := crossover(a, b)
child.mutate()
child.calcFitness(target) next[i] = child
}
return next
}
人口的下一代必须继承基因的价值
然后，下一代的child由两个随机挑选的生物体杂交而成，继承了两个生物体的DNA。
func crossover(d1 Organism, d2 Organism) Organism {
child := Organism{
DNA: make([]byte, len(d1.DNA)),
Fitness: 0,
}
mid := rand.Intn(len(d1.DNA))
for i := 0; i < len(d1.DNA); i++ {
if i > mid {
child.DNA[i] = d1.DNA[i]
} else {
child.DNA[i] = d2.DNA[i]
} }
return child
}
对于交叉，我简单地选择了一个中点mid，并使用第一个生物体的第一个mid字节和第二个生物体的剩余字节。
随机变异每一代
在一个新的子有机体从两个亲代有机体中复制出来后，我们观察子有机体是否发生了突变。
func (d *Organism) mutate() {
for i := 0; i < len(d.DNA); i++ {
if rand.Float64() < MutationRate {
d.DNA[i] = byte(rand.Intn(95) + 32)
}
}
}
这里的突变仅仅意味着确定一个随机生成的数是否低于MutationRate。为什么我们需要突变子生物体？如果突变从未发生，群体中的 DNA 将始终保持与原始群体相同。这意味着如果原始种群没有所需的特定基因(值),将永远无法获得最佳结果。就像在这个例子中，如果字母t在初始种群中根本找不到，那么无论我们经历多少代，我们都无法得出这个引用。换句话说，没有突变，自然选择就不起作用。
更专业地说，突变让我们摆脱局部最大值，以找到全局最大值。如果我们将遗传算法视为一种寻找最优解的机制，如果我们没有突变，一旦找到局部最大值，该机制将简单地停留在这一点上，而不会继续寻找全局最大值。突变可以使群体跳出局部最大值，因此为算法提供了继续寻找全局最大值的机会。
一旦我们检查了突变，我们就计算出子生物体的适应度，并将其插入到下一代种群中。
这就是遗传算法的全部内容！现在让我们把它们放在main函数中。
func main() {
start := time.Now()
rand.Seed(time.Now().UTC().UnixNano()) target := []byte("To be or not to be")
population := createPopulation(target) found := false
generation := 0
for !found {
generation++
bestOrganism := getBest(population)
fmt.Printf("\r generation: %d | %s | fitness: %2f", generation, string(bestOrganism.DNA), bestOrganism.Fitness) if bytes.Compare(bestOrganism.DNA, target) == 0 {
found = true
} else {
maxFitness := bestOrganism.Fitness
pool := createPool(population, target, maxFitness)
population = naturalSelection(pool, population, target)
} }
elapsed := time.Since(start)
fmt.Printf("\nTime taken: %s\n", elapsed)
}
在主函数中，我们经历了几代，每一代我们都试图找到最适合的有机体。如果最适合的生物体的基因与目标相同，我们就会找到答案。
现在运行软件程序！你花了多长时间？
因为初始群体是随机生成的，所以每次您都会得到不同的答案，但大多数情况下，我们可以在不到一秒钟的时间内完成报价！如果我们不得不强行计算的话，这与 934 万亿年有很大的不同。
进化中的蒙娜丽莎
进化莎士比亚似乎很简单。毕竟只是一串。来点不同的怎么样，比如说一个图像？还是有史以来最著名的画作，达芬奇的《蒙娜丽莎》?我们能进化吗？
Mona Lisa (public domain)
让我们给它同样的待遇。我们将从定义生物体来代表蒙娜丽莎的照片开始。
定义由一个或多个 DNA 组成的有机体
我们的 DNA 不再是字节数组，而是来自image标准库的结构。
type Organism struct {
DNA *image.RGBA
Fitness int64
}
从最初的生物群体开始
像以前一样，让我们先看看如何创造一个有机体。
func createOrganism(target *image.RGBA) (organism Organism) {
organism = Organism{
DNA: createRandomImageFrom(target),
Fitness: 0,
}
organism.calcFitness(target)
return
}
我们不创建随机字节数组，而是调用另一个函数来创建随机图像。
func createRandomImageFrom(img *image.RGBA) (created *image.RGBA) {
pix := make([]uint8, len(img.Pix))
rand.Read(pix)
created = &image.RGBA{
Pix: pix,
Stride: img.Stride,
Rect: img.Rect,
}
return
}
一个image.RGBA struct 由一个字节数组Pix ( byte和uint8是一回事)、一个Stride和一个Rect组成。对我们来说重要的是Pix，我们使用相同的Stride和Rect作为目标图像(这是蒙娜丽莎的图像)。幸运的是，math/rand标准库有一个名为Read的方法，可以用随机字节很好地填充一个字节数组。
你可能会好奇，那么我们在这里讨论的字节数组有多大呢？Pix无非是用 4 个字节代表一个像素的字节数组(R，G，B，A 各用一个字节表示)。对于一个 800 x 600 的图像，我们谈论的是每个图像中有 192 万字节！为了保持程序相对较快的速度，我们将使用一个 67 x 100 大小的较小图像，它给出了一个 26，800 字节的数组。如果你到现在还没有意识到，这和我们在上一个程序中使用的 18 字节相差甚远。
此外，你可能会意识到，因为每个像素现在是随机着色的，我们最终会得到一个彩色的静态雪花图案。
Randomly generate image
我们继续吧。
找到生物体的适合度
有机体的适合度是两个图像之间的差异。
// calculate the fitness of an organism
func (o *Organism) calcFitness(target *image.RGBA) {
difference := diff(o.DNA, target)
if difference == 0 {
o.Fitness = 1
}
o.Fitness = difference}// find the difference between 2 images
func diff(a, b *image.RGBA) (d int64) {
d = 0
for i := 0; i < len(a.Pix); i++ {
d += int64(squareDifference(a.Pix[i], b.Pix[i]))
}
return int64(math.Sqrt(float64(d)))
}// square the difference between 2 uint8s
func squareDifference(x, y uint8) uint64 {
d := uint64(x) - uint64(y)
return d * d
}
为了找到区别，我们可以回到勾股定理。如果你记得，我们可以找到两点之间的距离，如果我们平方x和y值的差，将它们相加，然后平方根结果。
Pythagorean theorem
给 2 个点a (x1，y1)和b (x2，y2)，a和b之间的距离d为:
那是在二维空间。在三维空间中，我们简单地做两次勾股定理，在四维空间中，我们做三次。一个像素的 RGBA 值本质上是 4 维空间中的一个点，因此为了找出两个像素之间的差异，我们对两个像素的r、g、b和a值之间的差异求平方，将它们相加，然后对结果求平方根。
这就是 2 个像素的区别。为了找出所有像素之间的差异，我们只需将所有结果加在一起，就可以得到最终的差异。因为Pix本质上是一个包含连续 RGBA 值的长字节数组，所以我们可以使用一个简单的快捷方式。我们简单地平方图像和目标中每个对应字节之间的差异，然后将它们全部相加，并对最终数字求平方根，以找到两个图像之间的差异。
作为参考，如果两个图像完全相同，差异将为 0，如果两个图像完全相反，差异将为 26，800。换句话说，最适合的生物体应该具有 0 的适合度，并且数字越高，生物体越不适合。
选择最适合的生物，给它们更高的繁殖机会
我们仍然使用繁殖池机制，但有所不同。首先，我们将人群从最佳健康状态到最差健康状态进行排序。然后我们把最好的生物放入繁殖池。我们使用一个参数PoolSize来表示我们希望池中有多少最适合的生物体。
为了计算出应该把什么放入繁殖池，我们减去每一个最好的生物体，最不适合的生物体在最上面。这在最好的生物体之间产生了不同的等级，并且根据该不同的等级，我们将该生物体的相应数量的拷贝放入繁殖池中。例如，如果最适合的生物体和最不适合的生物体之间的差异是 20，我们将 20 个生物体放入育种池。
如果最适合的生物之间没有差异，这意味着种群是稳定的，我们不能真正创造一个合适的繁殖池。为了克服这一点，如果差值为 0，我们将池设置为整个群体。
func createPool(population []Organism, target *image.RGBA) (pool []Organism) {
pool = make([]Organism, 0)
// get top best fitting organisms
sort.SliceStable(population, func(i, j int) bool {
return population[i].Fitness < population[j].Fitness
})
top := population[0 : PoolSize+1]
// if there is no difference between the top organisms, the population is stable
// and we can't get generate a proper breeding pool so we make the pool equal to the
// population and reproduce the next generation
if top[len(top)-1].Fitness-top[0].Fitness == 0 {
pool = population
return
}
// create a pool for next generation
for i := 0; i < len(top)-1; i++ {
num := (top[PoolSize].Fitness - top[i].Fitness)
for n := int64(0); n < num; n++ {
pool = append(pool, top[i])
}
}
return
}
从所选择的最适合的生物体中创建种群的下一代
有了游泳池之后，我们需要创造下一代。这里自然选择的代码和前面的程序没有什么不同，所以我们在这里就不展示了。
人口的下一代必须继承基因的价值
crossover函数略有不同，因为孩子的 DNA 不是字节数组，而是 image.RGBA。实际的交叉机制在Pix上工作，而是像素的字节数组。
func crossover(d1 Organism, d2 Organism) Organism {
pix := make([]uint8, len(d1.DNA.Pix))
child := Organism{
DNA: &image.RGBA{
Pix: pix,
Stride: d1.DNA.Stride,
Rect: d1.DNA.Rect,
},
Fitness: 0,
}
mid := rand.Intn(len(d1.DNA.Pix))
for i := 0; i < len(d1.DNA.Pix); i++ {
if i > mid {
child.DNA.Pix[i] = d1.DNA.Pix[i]
} else {
child.DNA.Pix[i] = d2.DNA.Pix[i]
} }
return child
}
随机变异每一代
mutate功能也相应不同。
func (o *Organism) mutate() {
for i := 0; i < len(o.DNA.Pix); i++ {
if rand.Float64() < MutationRate {
o.DNA.Pix[i] = uint8(rand.Intn(255))
}
}
}
现在我们已经有了所有的东西，我们把它们放在main函数中。
func main() {
start := time.Now()
rand.Seed(time.Now().UTC().UnixNano())
target := load("./ml.png")
printImage(target.SubImage(target.Rect))
population := createPopulation(target) found := false
generation := 0
for !found {
generation++
bestOrganism := getBest(population)
if bestOrganism.Fitness < FitnessLimit {
found = true
} else {
pool := createPool(population, target)
population = naturalSelection(pool, population, target)
if generation%100 == 0 {
sofar := time.Since(start)
fmt.Printf("\nTime taken so far: %s | generation: %d | fitness: %d | pool size: %d",
sofar, generation, bestOrganism.Fitness, len(pool))
save("./evolved.png", bestOrganism.DNA)
fmt.Println()
printImage(bestOrganism.DNA.SubImage(bestOrganism.DNA.Rect))
}
}
}
elapsed := time.Since(start)
fmt.Printf("\nTotal time taken: %s\n", elapsed)
}
现在运行它看看。你得到了什么？
Start evolving Mona Lisa
按照我设置的参数，当我运行它时，我通常会以 19，000 左右的健身值开始。平均来说，我需要 20 多分钟才能达到低于 7500 的体能。
After about 20 minutes
这是一系列随着时间推移而产生的图像:
Evolving Mona Lisa
用圆形和三角形进化蒙娜丽莎
我通过在一张图片上画圆圈和三角形来玩《蒙娜丽莎》的演变。结果没有那么快，图像也没有那么明显，但它显示了实际发生的一瞥。你可以从资源库中找到剩下的代码，自己调整参数，看看是否能得到更好的图片，但这里有一些我得到的图片。
蒙娜丽莎三角形
蒙娜丽莎的圆圈
玩得开心！
在终端上显示图像
你可能在我的截图中注意到了，我实际上在终端上显示了图像。我本来可以创建一个 web 应用程序来展示这一点，但是我想让事情简单得多，所以我想直接在终端上显示图像。虽然终端控制台不是您通常期望显示图像的地方，但实际上有几种方法可以做到这一点。
我选了一个最简单的。我正好用了优秀的 iTerm2 ，这是 MacOS 中默认终端应用的替代品，iTerm2 中有一个有趣的 hack，可以显示图像。
诀窍是这样的——如果你可以用 Base64 编码你的图像，你可以用一个特殊的命令把图像打印到终端上。下面是实现这一点的 Go 代码，但是您也可以用任何其他语言实现这一点。上面的文档中有几个脚本展示了如何使用简单的 shell 脚本来实现这一点。
func printImage(img image.Image) {
var buf bytes.Buffer
png.Encode(&buf, img)
imgBase64Str := base64.StdEncoding.EncodeToString(buf.Bytes())
fmt.Printf("\x1b]1337;File=inline=1:%s\a\n", imgBase64Str)
}
不幸的是，这意味着如果你在 iTerm2 之外的任何地方运行这段代码，你将看不到图像的变化。但是，您可以随时调整输出，以便每隔几代就捕获一次输出。
密码
这篇文章中的所有代码和图片都可以在这里找到:https://github.com/sausheong/ga
参考
示例代码受到了以下工作的启发:
- 丹尼尔·席夫曼的优秀著作代码的本质【http://natureofcode.com】T2——这是一本很棒且容易理解的读物！我还用 Java 编写了一些丹尼尔的代码，并将其转换成莎士比亚引用算法
- 罗杰·约翰逊在《遗传编程:蒙娜丽莎的进化》一文中的出色工作https://rogerjohansson . blog/2008/12/07/Genetic-Programming-evolution-of-Mona-Lisa/——尽管我最终使用了一种完全不同的方式来完成遗传算法，但他的原创工作给了我灵感，让我使用蒙娜丽莎，并尝试用三角形和圆形来完成它
最大似然估计简介
我第一次听到有人使用术语最大似然估计时，我去了谷歌，并找出了它的意思。然后我去维基百科了解它真正的意思。我得到了这个:
在统计学中，最大似然估计 ( MLE )是一种估计统计模型给定观测值的参数的方法，通过寻找使给定参数下进行观测的似然最大化的参数值。MLE 可以被视为最大后验概率估计 (MAP)的特殊情况，其假设参数的均匀 先验分布，或者被视为忽略先验分布的 MAP 的变体，因此其不规则。
好吗？？
为了避免你为理解 MLE 并将其融入你的数据科学工作流程、理念和项目而绞尽脑汁，我编辑了这个指南。下面，我们将:
- 将概率上下文设置为 MLE
- 钻研所需的数学
- 看看 MLE 在 Python 中是如何工作的
- 利用 MLE 探索数据科学的最佳实践
不过先来点 xkcd :
常客 vs .贝叶斯
这很有趣(如果你遵循这个奇怪的幽默领域)，而且这两个阵营之间的差异基本上是正确的。不要介意我们的太阳进入新星并不是一个真正可重复的实验——抱歉，常客们！—我们可以概括为，对于真实的观察和研究，两个阵营通常会得出类似的结论，但当研究设计或数据开始变得棘手时，就会有很大的不同。
简而言之，MLE 帮助我们回答了这个问题:
哪些参数/系数最适合我的模型？
有趣的是，你可以用这两种观点中的任何一种来解释 MLE 为什么有效！因为，虽然最大似然估计给出了一个现场估计——这在 frequentist 输出中很常见——但它可以被视为最大后验概率 (MAP)估计的一个特例，在这里我们使用了一个天真的先验知识，并且从不费心去更新它。
设置我们的问题
今天为了接近 MLE，让我们从贝叶斯的角度出发，使用贝叶斯定理将我们的问题框定如下:
P(β∣y) = P(y∣β) x P(β) / P(y)
或者，用英语说:
posterior = likelihood x prior / evidence
我们可以有效地忽略prior和evidence,因为——给定均匀先验分布的 Wiki 定义——所有系数值都是同等可能的。而且所有数据值的概率(假设连续)都是相等的，而且基本上为零。
所以，在实际的英语中:给定一些特定系数的概率，我看到一些结果，涉及到以完全相反的方式提出问题。这很有帮助，因为这个问题更容易解决。
概率和可能性
从今以后，我们将在代码中引入可能性的概念，或L。为了理解其中的区别，我将从 Randy Gallistel 的精彩文章中摘录:
概率和可能性之间的区别是非常重要的:概率与可能的结果相联系；可能性附属于假设。
可能的结果是互斥的和穷尽的。假设我们让一个受试者预测 10 次投掷硬币的结果。只有 11 个可能的结果(0 到 10 个正确的预测)。实际的结果总是一个又一个可能的结果。因此，可能结果的概率总和必须为 1。
与结果不同，假设既不是互斥的，也不是穷尽的。假设我们测试的第一个对象正确预测了 10 个结果中的 7 个。我可能会假设受试者只是猜测，你可能会假设受试者可能有某种千里眼，你的意思是受试者可能被期望以稍高于概率的长期正确预测结果。这些是不同的假设，但并不互相排斥，因为你说“可能是”的时候就对冲了。你因此允许你的假设包括我的。用专业术语来说，我的假设嵌套在你的假设中。其他人可能会假设受试者有很强的透视能力，观察到的结果低估了她下一次预测正确的可能性。另一个人可以假设一些完全不同的东西。人们可以接受的假设是无限的。
我们赋予可能性的一系列假设受到我们虚构它们的能力的限制。实际上，我们很少能确信我们已经设想了所有可能的假设。我们关心的是估计实验结果在多大程度上影响我们和其他人目前所接受的假设的相对可能性。因为我们通常不考虑所有的替代假设，而且因为一些假设嵌套在其他假设中，我们附加到假设上的可能性本身没有任何意义；只有相对可能性——也就是两个可能性的比率——才有意义。
太神奇了！谢谢你兰迪。
MLE 是频率主义者，但可以从贝叶斯的角度出发:
- 常客可以声称 MLE，因为它是一个逐点估计(不是一个分布)，并且它假设没有先验分布(技术上，不知情或一致)。
- 此外，MLE 没有给出真实参数值的 95%概率区域。
- 然而，最大似然法是一种特殊形式的映射，并且使用了似然的概念，这是贝叶斯哲学的核心。
当心天真或一致的先验假设！！您可能会将数据错误地归因于一个极不可能的模型。你可能会成为辛普森悖论的受害者，如下图。你很容易被小样本欺骗。
Simpson’s Paradox
以上都是常客和数据科学家必须处理或意识到的问题，所以 MLE 没有什么本质上更糟糕的。
回到我们的问题
所以如果p(y|*β*)相当于 L(*β*|y)，那么p(y_1,y_2,...,y_n|*β*)相当于L(*β*|y_1,y_2,...,y_n)。另外，请记住，我们可以将独立概率相乘，就像这样:
p(A,B) = p(A)p(B)
我们越来越近了！这是我们当前的设置:
L(*β*|y1,y2,…,yn) = p(y1|*β*)p(y2|*β*),…,p(yn|*β*) = ∏p(yi|*β*)
右边的部分看起来像是我们可以最大化的东西:
Initial Cost Function
但是我们可以做得更好！用自然对数把我们的积函数变成和函数怎么样？日志是单调变换，所以我们将简化我们的计算，但保持我们的最佳结果。
Halfway there!
我们最终的成本函数如下所示:
Ready to roll!
为了简单起见，让我们假设我们有一个回归问题，所以我们的结果是连续的。最大似然法对于离散结果的分类问题非常有效，但是我们必须使用不同的分布函数，这取决于我们有多少个类，等等。
现在，记住普通最小二乘法(OLS)等模型的一个中心假设是残差正态分布在均值零附近，我们拟合的 OLS 模型实际上成为了最大期望值y的体现。而我们的概率分布是… 正态！
y is normally distributed around our ŷ
因为计算机在计算概率方面比我们好得多，所以我们将从这里转向 Python！
Python 中的 MLE
在您的数据科学建模管道中实现 MLE 可能非常简单，有多种方法。下面是一个你可以偷着开始的方法。
设置
如果导入正确的包，MLE 很容易:
# import libraries
import numpy as np, pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
import scipy.stats as statsimport pymc3 as pm3
import numdifftools as ndt
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel
%matplotlib inline
在此基础上，我们将生成遵循围绕基本事实函数的正态分布误差的数据:
# generate data
N = 100
x = np.linspace(0,20,N)
ϵ = np.random.normal(loc = 0.0, scale = 5.0, size = N)
y = 3*x + ϵdf = pd.DataFrame({‘y’:y, ‘x’:x})
df[‘constant’] = 1
最后，让我们使用 Seaborn 的 regplot 进行可视化:
# plot
sns.regplot(df.x, df.y);
我得到了下面的，你应该看到类似的东西。但是，请记住这里有随机性，我们没有使用种子:
Scatter plot with OLS line and confidence intervals
使用 Statsmodels 为 OLS 建模
由于我们创建了类似回归的连续数据，我们将使用sm.OLS来计算最佳系数和对数似然(ll)作为基准。
# split features and target
X = df[[‘constant’, ‘x’]]# fit model and summarize
sm.OLS(y,X).fit().summary()
我得到这个，并将记录拟合模型的系数:
注意constant接近于零，对于我们使用的地面真实发生器，特征x的beta接近于 3。
最大化 LL 以求解最佳系数
从这里开始，我们将使用软件包和自定义函数的组合，看看我们是否可以使用 MLE 方法计算相同的 OLS 结果。
因为scipy.optimize只有一个minimize方法，我们将最小化对数似然的负值。这甚至是他们推荐的！数学欺骗通常比重新发明轮子更快更容易！
我们可以构建一个简单的函数，一次完成回归输出的所有工作:
# define likelihood function
def MLERegression(params):
intercept, beta, sd = params[0], params[1], params[2] # inputs are guesses at our parameters
yhat = intercept + beta*x # predictions# next, we flip the Bayesian question
# compute PDF of observed values normally distributed around mean (yhat)
# with a standard deviation of sd
negLL = -np.sum( stats.norm.logpdf(y, loc=yhat, scale=sd) )# return negative LL
return(negLL)
现在我们有了一个成本函数，让我们初始化并最小化它:
# let’s start with some random coefficient guesses and optimize
guess = np.array([5,5,2])results = minimize(MLERegression, guess, method = ‘Nelder-Mead’,
options={‘disp’: True})--------------------------------------------------------------------
Optimization terminated successfully.
Current function value: 311.060386
Iterations: 111
Function evaluations: 195
让我们来看看结果:
results # this gives us verbosity around our minimization
# notice our final key and associated values…--------------------------------------------------------------------
final_simplex: (array([[0.45115297, 3.03667376, 4.86925122],
[0.45123459, 3.03666955, 4.86924261],
[0.45116379, 3.03667852, 4.86921688],
[0.45119056, 3.03666796, 4.8692127 ]]), array([300.18758478, 300.18758478, 300.18758478, 300.18758479]))
fun: 300.18758477994425
message: 'Optimization terminated successfully.'
nfev: 148
nit: 80
status: 0
success: True
x: array([0.45115297, 3.03667376, 4.86925122])
我们可以进一步清理:
# drop results into df and round to match statsmodels
resultsdf = pd.DataFrame({'coef':results['x']})
resultsdf.index=['constant','x','sigma']
np.round(resultsdf.head(2), 4)# do our numbers match the OLS model?
--------------------------------------------------------------------
你会注意到 OLS 和我很相配！你的结果会有所不同，因为我们没有使用随机种子。
MLE 的最佳实践
在我们更进一步之前，这可能是一个加强我们对 MLE 信任的好时机。作为我们的回归基线，我们知道，根据定义，普通最小二乘法是具有正态分布残差并满足线性回归的其他假设的连续结果的最佳线性无偏估计量。使用最大似然法来寻找我们的系数是否稳健？
是的！
- MLE 与 OLS 一致。
- 对于无限的数据，它将估计最佳的 β，并很好地近似它用于小而健壮的数据集。
- MLE 是高效的；没有一致的估计量具有比 MLE 更低的渐近均方误差。
所以看起来它完全复制了 OLS 的做法。那么…为什么用 MLE 而不是 OLS？
因为！
- MLE 对于回归和分类是可推广的！
- MLE 是高效的；如果你使用正确的分布，没有一致的估计比 MLE 有更低的渐近误差。
我们可以将 MLE 视为通过优化概率成本函数来拟合模型的模块化方式！
应用 MLE 的四个主要步骤:
- 定义可能性，确保您对回归或分类问题使用正确的分布。
- 取自然对数，将乘积函数简化为和函数。
- 最大化——或最小化目标函数的负值。
- 验证制服前科是一个安全的假设！否则，你可以将数据归因于一个生成函数或世界模型，它不符合简约法则。
在最大似然估计空间中有更多的东西，包括分类分布，使用贝叶斯统计软件包如PyMC3等。但是今天我们就讲到这里。
您如何在数据科学工作流程中使用 MLE？在下面评论，或者在 LinkedIn 或者 Twitter 上联系我！
特别感谢 查德·谢勒 的优秀同行点评。
神经网络简介
神经网络系列
神经网络系列—第 0 章
序言注释
这个故事是我正在创作的关于神经网络的新系列的第一部分。作为一名前研究人员，我喜欢让自己跟上研究领域的最新发展，我也喜欢写作。我最近一直在努力做到这一点，直到我想到了一个主意:如果我强迫自己通过写作来强迫自己阅读和调查更多的东西，会怎么样？所以我想出了这个从零开始的想法，深入研究神经网络领域的历史，从基础直到今天应用的复杂的神经网络架构。这个系列对我个人来说有两个挑战:填补知识上的空白(我很确定我不是什么都知道)，以及确保我可以开始一个想法并坚持到底而不放弃。
我决定将这一章命名为第 0 章，因为在这篇文章中，我将只解释该领域的最初历史步骤，为什么神经网络存在，它们用于什么，以及神经生物学和计算机科学之间的桥梁是什么。
神经网络
神经网络是一种人工智能连接主义方法，能够学习如何执行模式识别或过程控制任务，而无需显式编程。神经网络应用的一些例子可以是(1) 控制房间的温度，(2) 识别图片/视频中的人脸，(3) 生成文本序列。它们受到我们大脑工作方式的启发，通过连接和引起多个基本实体之间的交互来实现最终目标。这些实体被称为神经元(根据这些算法和我们大脑工作方式之间的整体相似性)。神经元可以与其他神经元相互作用，这取决于它们相互连接的方式以及它们的连接有多强。
这个过程是如何从生物学转化到计算机科学的？
从生物学上讲，神经元由 3 个主要部分组成:树突、核/细胞体和轴突。其完整结构如图 1 所示。
Figure 1: Structure of a neuron. Image taken from this lifehacker article
树突以电脉冲的形式接收来自其他神经元的信息。细胞体是接收和处理来自树突的信息的地方。这种处理的结果然后通过轴突传导，直到它到达末端(突触)。然后，信息通过一个称为突触过程的化学过程传递给另一个树突神经元。在整个人类生活中，许多神经元被创建和死亡，神经元之间的许多连接被创建，而其他连接则不复存在(这是一种过于简化的说法，但它让你了解了我们大脑的动态)。
根据赫布边理论，当它们之间的突触过程重复发生时，这些连接的强度会更大。换句话说，两个神经元一起放电越多，它们之间的联系就越强。另一方面，如果两个神经元从未一起被激发，它们的连接就会变弱。这条规则被称为赫布规则。
1943 年，麦卡洛克和皮茨提出了一个数学模型，能够捕捉神经元的这些生物特征。图 2 描述了神经元的这种数学模型的模式。这个被称为麦卡洛克-皮茨神经元。
Figure 2: Schema of the McCulloch-Pitts neuron
麦卡洛克-皮茨神经元有三个组成部分:
- 一组加权输入(x和 w )。它们分别代表输入和连接强度。赫布规则是关于刺激或减弱神经元之间的连接，取决于它对两个神经元一起放电的重要性。 w 的角色就是这个意思:描述关系的优势。它们被称为权重，因为它们量化了连接的权重。 x 映射到从突触中一个神经元的轴突末端传递到其他神经元的树突的输入。
- 一个加法器(即 z )。这代表负责将所有输入信号整理和聚集成统一电脉冲的核/细胞体，该电脉冲稍后将沿着轴突传输。用数学术语来说，麦卡洛克和皮茨将这种聚合转化为输入( x )和权重( w )之间的点积。但这意味着什么呢？让我们使用图 2 中的例子。向量 x 由 x₀ 、 x₁ 、 x ₂组成，而向量 w 具有元素 w₀ 、 w₁ 、 w ₂ 。那么 x 和 w 的点积就是:x₀×w₀﹢x₁×w₁﹢x₂×w₂.让我们用数学符号将这个逻辑形式化。
- 一个激活功能(*f(z)*块)。激活功能的作用是决定整理后的信号的输出，并基于整理后的信号检查神经元是否应该触发。神经元的这一部分是一个数学函数，它抓取 z 并基于该函数计算输出 y 。麦卡洛克和皮茨提出的方法是阈值函数。基本上定义一个阈值电平𝜃，如果 z 大于𝜃，输出将是 1 ，否则输出将是 0 。
让我们尝试一个端到端的例子。想象你有一个神经元，它有输入: x₀﹦2 ， x₁﹦-0.75 ， x ₂﹦ 2 ，还有权重 w₀﹦-1 ， w₁﹦-1 ， w ₂﹦ 1 。所以，要做的第一件事就是计算 z。根据之前定义的方程，z﹦2×(-1)﹢(-0.75)×(-1)﹢2×1﹦0.75。最后一步是计算激活函数 f(z) 。例如，让我们假设一个𝜃﹦为 0 的阈值函数。因为 0.75 > 0 我们的输出会是 1 。
我已经解释了最基本的神经网络是如何工作的。但真的是全部吗？你现在已经理解了网络如何基于外部输入或来自环境的传感器来计算其输出/预测/结果，但是在这个解释中描述的学习过程在哪里呢？答案是:无处！神经网络学习的方式是通过找到向量 w 的正确值，这些值将增加网络根据提供的输入正确猜测答案的次数。向量 w 在很大程度上是神经网络中的“内部”,因此尝试不同的值是我们可以用来训练网络的。那么这些权重 w 如何修改才能提高网络质量呢？我将在我的下一章:感知机网络中谈论这个。敬请期待！
感谢阅读！你喜欢这篇文章吗？非常感谢你的反馈，🗣。你可以随时在 Twitter 或 LinkedIn 上联系我，或者如果你对下一章的最新消息感兴趣，就在 Medium 上关注我😀。
一些相关阅读:
- 【1】神经网络简介、优势及应用作者Jahnavi Mahanta；
- 【2】神经网络简介作者约翰·奥拉芬瓦；
- [3] 对神经网络的可视化介绍作者shik har Sharma；
- [4] 人工神经网络——简介作者 Tharani Gnanasegaram 。
- 【5】麦卡洛克-皮茨神经元——人类第一个生物神经元的数学模型作者阿克谢·钱德拉·拉甘杜拉
神经网络系列简介—第 1 部分
神经网络系列简介(GINNS)。
介绍
神经网络和深度学习是计算机科学和技术行业的重要话题，它们目前为图像识别、语音识别和自然语言处理中的许多问题提供了最佳解决方案。最近，许多论文发表了人工智能，它可以学习绘画，建立 3D 模型，创建用户界面(pix2code)，有些可以根据一个句子创建图像，还有许多令人难以置信的事情正在使用神经网络进行日常工作。
我正在撰写这一系列关于神经网络和深度学习的帖子，我将指导您学习人工神经网络(ANN)的基本概念，向您展示从简单网络模拟与门到使用卷积神经网络(CNN)、递归神经网络(RNN)等解决图像识别任务的示例。代码将始终用 python 编写，有时会借助 Tensorflow(我不期望您是使用 Tensorflow 的专家，因为我将尝试详细解释代码)。
议程
- 神经网络导论(本帖)
- 与门神经网络(感知器)和异或门前馈神经网络(2 层)。
- 使用 CNN 进行首次数字识别
- 用 RNN 识别第一位数字
神经网络
神经网络，更确切地说是“人工”神经网络(ANN)的定义是由第一台神经计算机的发明者 Robert Hecht-Nielsen 博士提出的。他将神经网络定义为:
"…由许多简单、高度互联的处理元件组成的计算系统，这些元件通过对外部输入的动态响应来处理信息
或者，你也可以将人工神经网络视为计算模型，它是受人脑中生物神经网络处理信息的方式的启发。
生物动机和联系
大脑的基本计算单元是一个神经元。在人类神经系统中可以找到大约 860 亿个神经元，它们与大约 10 个⁴-10 个⁵ 突触相连。下图显示了一个生物神经元(左)和一个常见数学模型(右)的卡通图。
biological neuron (left) and a common mathematical model (right)
神经网络中的基本计算单元是神经元，通常称为节点或单元。它接收来自其他节点或外部源的输入，并计算输出。每个输入都有一个相关的
权重(w ),该权重根据其对其他输入的相对重要性进行分配。该节点对其输入的加权和应用一个函数。
这个想法是，突触强度(权重 w )是可以学习的，并控制影响的强度和方向:一个神经元对另一个神经元的兴奋性(正权重)或抑制性(负权重)。在基本模型中，树突将信号传递到细胞体，并在那里进行汇总。如果最终总和高于某个阈值，神经元就会触发，沿着其轴突发出一个尖峰信号。在计算模型中，我们假设尖峰脉冲的精确时间并不重要，只有放电的频率才能传递信息。我们用一个激活函数 (e.x sigmoid function) 来模拟神经元的放电频率，它代表了沿着轴突的尖峰的频率。
神经网络体系结构
从上面的解释中，我们可以得出结论，神经网络由神经元组成，从生物学上讲，神经元通过突触连接，信息在突触中流动(计算模型的权重)，当我们训练神经网络时，我们希望神经元在从数据中学习到特定模式时触发，我们使用激活函数对触发率进行建模。
但这并不是一切……
- 输入节点(输入层):这一层内不做任何计算，它们只是把信息传递给下一层(大部分时间是隐藏层)。节点块也被称为层。
- 隐藏节点(Hidden layer):中的隐藏层是进行中间处理或计算的地方，它们执行计算，然后将权重(信号或信息)从输入层传递到下一层(另一个隐藏层或输出层)。有一个没有隐藏层的神经网络是可能的，我稍后会解释这一点。
- ****输出节点(输出层):这里我们最后使用一个激活函数，它映射到所需的输出格式(例如用于分类的 softmax)。
- ****连接和权重:网络由连接组成，每个连接将神经元 I 的输出传输到神经元 j 的输入。从这个意义上来说， i 是 j 的前任，而 j 是 i 的继任者，每个连接都被赋予一个权重 Wij。
- 激活函数:给定一个输入或一组输入，节点的激活函数**定义该节点的输出。一个标准的计算机芯片电路可以被看作是一个激活功能的数字网络，激活功能可以是“开”(1)或“关”(0)，这取决于输入。这类似于神经网络中线性感知器的行为。然而，正是非线性激活函数允许这样的网络仅使用少量节点来计算非平凡的问题。在人工神经网络中，这个函数也称为传递函数。
- **学习规则:学习规则是修改神经网络参数的规则或算法，目的是使网络的给定输入产生有利的输出。这个学习过程通常相当于修改权重和阈值。
神经网络的类型
神经网络有许多类别，这些类别也有子类，这里我将列出最常用的类别，并使学习神经网络的过程变得简单。
1.前馈神经网络
前馈神经网络是一种人工神经网络，其中单元之间的连接不也不形成循环。在这个网络中，信息只在一个方向上移动，即从输入节点向前，通过隐藏节点(如果有的话)到达输出节点。网络中没有循环或环路。
我们可以区分两种类型的前馈神经网络:
1.1。单层感知器
这是最简单的前馈神经网络，不包含任何隐藏层，这意味着它只由单层输出节点组成。之所以说这是单一的，是因为当我们计算层数时，我们不包括输入层，原因是因为在输入层不进行任何计算，输入通过一系列权重直接馈入输出。
Simple Perceptron
1.2.多层感知器(MLP)
这类网络由多层计算单元组成，通常以前馈方式互连。一层中的每个神经元都与下一层的神经元有直接连接。在许多应用中，这些网络的单元应用 sigmoid 函数作为激活函数。MLP 非常有用，一个很好的原因是，它们能够学习非线性表示(大多数情况下，呈现给我们的数据不是线性可分的)，我们将在下一篇文章中展示的示例中回来分析这一点。
MLP
1.3.卷积神经网络(CNN)
卷积神经网络与普通神经网络非常相似，它们由具有可学习权重和偏差的神经元组成。在卷积神经网络(CNN 或 ConvNet 或移位不变量或空间不变量)中，单元连接模式受到视觉皮层组织的启发，单元在称为感受野的受限空间区域中对刺激做出反应。感受野部分重叠，覆盖整个视野。单位响应可以通过卷积运算进行数学近似。它们是使用最少预处理的多层感知器的变体。它们广泛应用于图像和视频识别、推荐系统和自然语言处理。CNN 需要大量数据进行训练。
CNN for image classification
2.递归神经网络
在递归神经网络(RNN)中，单元之间的连接形成了一个有向循环(它们向前传播数据，但也向后传播数据，从后面的处理阶段到前面的阶段)。这允许它展示动态的时间行为。与前馈神经网络不同，rnn 可以使用其内部存储器来处理任意输入序列。这使得它们适用于诸如未分段的、连接的手写识别、语音识别和其他通用序列处理器的任务。
常用的激活功能
每个激活函数(或非线性)取一个数字，并对其执行某个固定的数学运算。以下是一些您在实践中经常会发现的激活功能:
- 乙状结肠
- 谭
- ReLU
- 泄漏的 ReLU
有关更多详细信息，请遵循此激活功能帖子以及此cs 231n(CNN)课程笔记。
人工神经网络和大脑
人工神经网络不像我们的大脑那样工作，人工神经网络只是简单粗糙的比较，生物网络之间的连接比那些由人工神经网络架构实现的要复杂得多，记住，我们的大脑要复杂得多，我们需要从中学习更多。关于我们的大脑，有许多事情我们不知道，这也使得我们很难知道我们应该如何在人类层面上模拟一个人工大脑来进行推理。每当我们训练一个神经网络时，我们希望我们的模型学习最佳权重(w ),该权重在给定输入信号或信息(x)的情况下最好地预测期望的结果(y)。
接下来是什么？
在下一篇文章中，我将带来一些代码示例，并解释我们如何训练神经网络。敬请期待！
参考资料:
- 神经网络和深度学习
- http://cs231n.github.io/neural-networks-1/
如果你喜欢这些作品，请留下你的掌声👏推荐这篇文章，让其他人也能看到。**
OCR 简介
想了解更多？参观www.Shibumi-ai.com
在这里阅读这个帖子的重访版本
介绍
光学字符识别(OCR)是最早解决的计算机视觉任务之一，因为在某些方面它不需要深度学习。因此，在 2012 年深度学习热潮之前就有不同的 OCR 实现，有些甚至可以追溯到 1914 年(！).
这让很多人以为 OCR 挑战是**“解决了**”，就不再有挑战性了。另一个来自类似来源的信念是，OCR 不需要深度学习，或者换句话说，使用深度学习进行 OCR 是一种矫枉过正。
任何实践计算机视觉或一般意义上的机器学习的人都知道，没有所谓的已解决的任务，这个案例也不例外。相反，OCR 仅在非常具体的用例上产生非常好的结果，但是一般来说，它仍然被认为是具有挑战性的。
此外，对于某些不需要深度学习的 OCR 任务，确实有很好的解决方案。然而，要真正迈向更好、更通用的解决方案，深度学习将是强制性的。
我为什么要写 OCR？
像我的许多作品/文章一样，这也是作为客户的项目开始的。我被要求解决一个特定的 OCR 任务。
在从事这项工作期间和之后，我得出了一些值得分享的结论和见解。此外，在集中精力完成一项任务后，很难停下来扔掉它，所以我继续我的研究，并希望获得一个更好、更通用的解决方案。
你会在这里找到什么
在这篇文章中，我将探索一些用于解决不同 OCR 任务的策略、方法和逻辑，并分享一些有用的方法。在最后一部分，我们将用代码解决一个现实世界的问题。这不应该被认为是一个详尽的评论(不幸的是),因为这种方法的深度、历史和广度对于这种博客文章来说太广了。
然而，和往常一样，我不会放过你参考文章、数据集、知识库和其他相关的博客文章。
OCR 的类型
如前所述，OCR 有不止一种含义。从最普遍的意义上来说，它指的是从每一个可能的图像中提取文本，无论是一本书的标准印刷页，还是一张随机的带有涂鸦的图像( in the wild )。在这之间，你可能会发现许多其他任务，比如阅读车牌，无人机器人验证码，街道标志等。
尽管这些选项中的每一个都有自己的困难，但显然“在野外”的任务是最难的。
Left: Printed text. Right: text in the wild
从这些例子中，我们可以得出 OCR 任务的一些属性:
- 文字密度:在打印/书写的页面上，文字密度大。然而，给定一个只有一个街道标志的街道图像，文本是稀疏的。
- 文本的结构:页面上的文本是结构化的，大部分是严格的行，而野生的文本可能会以不同的旋转方式散落在各处。
- 字体:印刷字体更容易，因为它们比嘈杂的手写字符更有结构。
- **字符类型:**文本可能以不同的语言出现，这些语言可能彼此非常不同。此外，文本的结构可能不同于数字，例如门牌号等。
- 伪像:很明显，室外照片比舒适的扫描仪要嘈杂得多。
- 位置:一些任务包括裁剪/居中的文本，而在其他情况下，文本可能位于图像中的随机位置。
数据集/任务
SVHN
一个很好的起点是街景门牌号数据集。顾名思义，这是从谷歌街景中提取的一组门牌号数据。任务难度中等。数字以各种形状和书写风格出现，然而，每个门牌号位于图像的中间，因此不需要检测。这些图像的分辨率不是很高，它们的排列可能有点奇特。
车牌
另一个常见的挑战是车牌识别，这在实践中不是很难，也不是很有用。这个任务和大多数 OCR 任务一样，需要检测车牌，然后识别它的字符。由于车牌的形状相对恒定，一些方法在实际识别数字之前使用简单的整形方法。以下是一些来自网络的例子:
OpenALPR example. with car type a s abonus
- OpenALPR 是一个非常强大的工具，不需要深度学习，可以识别不同国家的车牌
- 这个 repo 提供了 CRNN 模型的实现(将进一步讨论)来识别韩国车牌。
- 监督. ly，一家数据应用公司，写了关于使用他们的工具生成的人工数据来训练一个车牌识别器(人工数据也将被进一步讨论)
验证码
由于互联网上到处都是机器人，区分它们和真正人类的常见做法是视觉任务，特别是文本阅读，又名验证码。这些文本中的许多是随机的和扭曲的，这应该使计算机更难阅读。我不确定开发验证码的人是否预测到了计算机视觉的进步，但是今天大多数文本验证码并不难解决，尤其是如果我们不试图一次解决所有的验证码。
Facebook knows how to make challenging CAPTCHAs
亚当·盖特基提供了一个不错的教程，用深度学习解决一些验证码，其中包括再次合成人工数据。
PDF OCR
OCR 最常见的情况是印刷/pdf OCR。打印文档的结构化特性使得解析它们变得更加容易。大多数 OCR 工具(例如 Tesseract )主要是为了解决这个任务，并获得良好的结果。因此，我不会在这篇文章中对这个任务做过多的阐述。
野外的 OCR
这是最具挑战性的 OCR 任务，因为它将所有一般的计算机视觉挑战(如噪声、光照和伪像)引入 OCR。此任务的一些相关数据集是 coco-text 和 SVT 数据集，它们再次使用街景图像来提取文本。
COCO text example
合成文本
SynthText 不是一个数据集，甚至可能不是一个任务，但提高训练效率的一个好主意是人工数据生成。由于文本的扁平性质，在图像上随意添加字符或单词看起来比其他任何物体都要自然。
我们之前已经看到了一些为简单任务生成的数据，比如验证码和车牌。在野外生成文本稍微复杂一点。该任务包括考虑图像的深度信息。幸运的是，SynthText 是一个很好的作品，它接收带有上述注释的图像，并智能地散布文字(来自新闻组数据集)。
SynthText process illustration: top right is the segmentation of an image, bottom right is the depth data. Bottom left is a surface analyses of the image, which according to text is sprinkled on the image.
为了使“分散”的文本看起来真实而有用，SynthText 库为每幅图像设置了两个遮罩，一个是深度遮罩，另一个是分段遮罩。如果您喜欢使用自己的图像，也应该添加这些数据
- 建议查看 repo ，自己生成一些图像。你应该注意，回购使用了 opencv 和 maptlotlib 的一些过时版本，所以一些修改可能是必要的。
Mnist
虽然这并不是一个真正的 OCR 任务，但是要写 OCR 而不包括 Mnist 例子是不可能的。最广为人知的计算机视觉挑战实际上并不是一项经过深思熟虑的 OCR 任务，因为它一次只包含一个字符(数字)，而且只有 10 个数字。然而，它可能暗示了为什么 OCR 被认为是容易的。此外，在一些方法中，每个字母将被单独检测，然后 Mnist like(分类)模型变成 relevantץ
战略
正如我们已经看到和暗示的，文本识别主要是一个两步任务。首先，您希望检测图像中的文本外观，可能是密集的(如在打印的文档中)或稀疏的(如在野外的文本)。
在检测行/字级别之后，我们可以再次从一大组解决方案中进行选择，这些解决方案通常来自三种主要方法:
- 经典的计算机视觉技术。
- 专业化的深度学习。
- 标准深度学习方法(检测)。
让我们逐一检查一下:
1.经典的计算机视觉技术
如前所述，计算机视觉解决各种文字识别问题由来已久。你可以在网上找到很多例子:
经典 CV 方法通常声称:
- 应用滤镜让人物从背景中凸显出来。
- 应用轮廓检测逐个识别字符。
- 应用图像分类识别字符
显然，如果第二部分做得好，第三部分很容易，无论是模式匹配还是机器学习(例如 Mnist)。
然而，轮廓检测对于概括来说是相当具有挑战性的。它需要大量的手动微调，因此在大多数问题中变得不可行。例如，让我们将一个简单的计算机视觉脚本从应用到来自 SVHN 数据集的一些图像上。第一次尝试我们可能会取得很好的结果:
但是当角色之间的距离越来越近时，事情就开始变得不可收拾了:
我艰难地发现，当你开始摆弄参数时，你可能会减少这样的错误，但不幸的是会引起其他错误。换句话说，如果你的任务不简单，这些方法就不合适。
2.专门的深度学习方法
大多数成功的深度学习方法都在通用性方面表现出色。然而，考虑到上述属性，专用网络可能非常有用。
我将在这里检查一些著名方法的非详尽示例，并对介绍它们的文章做一个快速总结。与往常一样，每篇文章都以“任务 X(文本识别)最近获得关注”为开头，并继续详细描述他们的方法。仔细阅读这些文章会发现，这些方法是从以前的深度学习/文本识别工作中收集的。
结果也描述得很透彻，但是由于设计上的许多差异(包括数据集的微小差异)，实际比较是不可能的。真正了解这些方法在您的任务中的性能的唯一方法，是获取它们的代码(最好到最坏:找到官方的回购，找到非官方但评价很高的回购，自己实现)并在您的数据上尝试。
因此，我们总是更喜欢附有好的回购的文章，如果可能的话，甚至是演示。
东
EAST (高效准确的场景文本检测器)是一种简单而强大的文本检测方法。使用专门的网络。
与我们将要讨论的其他方法不同，只限于文本检测(不是实际的识别)，但是它的健壮性值得一提。
另一个优势是它也被添加到 open-CV 库(从版本 4 开始)所以你可以很容易地使用它(见教程这里)。这个网络实际上是众所周知的 U 形网的一个版本，它可以很好地检测大小不同的特征。这个网络的底层前馈“stem”(如文中杜撰，见下图)可能非常— 文中使用 PVANet ，然而 opencv 实现使用 Resnet 。显然，它也可以被预先训练(例如使用 imagenet)。与 U-Net 一样，特征是从网络的不同层次提取的。
最后，网络允许两种类型的输出旋转边界框:带有旋转角度(2X2+1 个参数)的标准边界框或“四边形”,它只是带有所有顶点坐标的旋转边界框。
如果现实生活的结果就像上面的图片一样，那么识别这些文字并不需要太多的努力。然而，现实生活的结果并不完美。
CRNN
卷积回归神经网络是 2015 年的一篇文章，其中提出了一种混合(或三混合？)端到端架构，旨在以三步方法捕获单词。
其思想如下:第一层是标准的全卷积网络。网络的最后一层被定义为特征层，并划分为“特征列”。在下图中可以看到，每个这样的特征列是如何表示文本中的某一部分的。
然后，特征列被输入到深度双向 LSTM 中，后者输出一个序列，用于查找字符之间的关系。
最后，第三部分是转录层。它的目标是获取杂乱的字符序列，其中一些字符是冗余的，而另一些是空白的，并使用概率方法来统一和解释它。
这种方法叫做 CTC 丢失，在这里可以读到。该层可以与/或不与预定义的词典一起使用，这可以促进单词的预测。
本文在固定文本词典的情况下达到了高准确率(> 95%)，在没有固定文本词典的情况下达到了不同的成功率。
STN-net/参见
参见 —半监督的端到端场景文本识别，是 Christian Bartzi 的作品。他和他的同事应用真正的端到端策略来检测和识别文本。他们使用非常弱的监督(他们称之为半监督，与通常意义不同)。因为他们训练网络只使用文本注释(没有边界框)。这允许他们使用更多的数据，但使他们的训练过程非常具有挑战性，他们讨论了不同的技巧来使其工作，例如，不要在超过两行文本的图像上训练(至少在训练的第一阶段)。
这份文件有一个更早的版本，叫做 STN OCR 。在最后的论文中，研究人员改进了他们的方法和陈述，此外，由于结果的高质量，他们更加强调他们方法的通用性。
SEE strategy
STN-OCR 的名字暗示了使用空间转换器的策略。
他们在中训练两个级联网络，其中第一个网络，即变换器，学习对图像的变换，以输出更容易解释的子图像。
然后，另一个 LSTM 在前的前馈网络(嗯…好像我们以前见过)识别文本。
研究人员在这里强调了使用 resnet 的重要性(他们使用了两次),因为它提供了到早期层的“强”传播。然而，这种做法如今已被广泛接受。
不管怎样，这都是一个有趣的尝试。
3.标准深度学习方法
正如标题所暗示的，在检测到“单词”后，我们可以应用标准的深度学习检测方法，如 SSD、YOLO 和 Mask RCNN。因为网上有太多的信息，所以我不打算详细说明这些方法。
我必须说这是我目前最喜欢的方法，因为我喜欢深度学习的“端到端”哲学，在这种哲学中，你可以应用一个强大的模型，通过一些调整可以解决几乎所有问题。在这篇文章的下一部分，我们将看到它实际上是如何工作的。
然而，当涉及到密集、相似的类时，SSD 和其他检测模型会受到挑战，如这里的所述。我觉得这有点讽刺，因为事实上，深度学习模型发现识别数字和字母比识别更具挑战性和复杂的物体(如狗、猫或人)要困难得多。它们往往达不到预期的精度，因此，专门的方法蓬勃发展。
实际例子
所以说了这么多之后，是时候动手了，试着自己做些模特。我们将尝试解决 SVHN 任务。SVHN 数据包含三个不同的数据集:训练、测试和额外。差异不是 100%清楚，但是最大的额外的数据集(约 50 万样本)包括更容易识别的图像。因此，为了这次拍摄，我们将使用它。
要准备该任务，请执行以下操作:
- 你需要一台 Tensorflow≥1.4、Keras≥2 的基本 GPU 机器
- 从这里克隆 SSD_Keras 项目。
- 从这里下载 coco 数据集上预先训练好的 SSD300 模型。
- 在此从克隆该项目的回购…
- 下载extra.tar.gz文件，其中包含 SVHN 数据集的额外图像。
- 更新此项目报告中 json_config.json 中的所有相关路径。
为了有效地遵循流程，您应该阅读以下说明，并从项目的 repo 中运行 ssd_OCR.ipynb 笔记本。
而且…你已经准备好开始了！
步骤 1:解析数据
不管你喜不喜欢，但是在探测任务中没有“黄金”格式的数据表示。一些众所周知的格式有:coco、via、pascal、xml。还有更多。例如，SVHN 数据集用晦涩的进行了注释。mat 格式。幸运的是，这个要点提供了一个巧妙的 read_process_h5 脚本来转换。mat 文件转换为标准 json，您应该更进一步，将其进一步转换为 pascal 格式，如下所示:
def json_to_pascal(json, filename): #filename is the .mat file
# convert json to pascal and save as csv
pascal_list = []
for i in json:
for j in range(len(i['labels'])):
pascal_list.append({'fname': i['filename']
,'xmin': int(i['left'][j]), 'xmax': int(i['left'][j]+i['width'][j])
,'ymin': int(i['top'][j]), 'ymax': int(i['top'][j]+i['height'][j])
,'class_id': int(i['labels'][j])})
df_pascal = pd.DataFrame(pascal_list,dtype='str')
df_pascal.to_csv(filename,index=False)*p = read_process_h5(file_path)*json_to_pascal(p, data_folder+'pascal.csv')
现在我们应该有一个更加标准的 pascal.csv 文件，它将允许我们继续前进。如果转换太慢，你应该注意到我们不需要所有的数据样本。~10K 就够了。
第二步:看数据
在开始建模过程之前，您最好对数据进行一些探索。我只为健全性测试提供了一个快速功能，但是我建议您做一些进一步的分析:
def viz_random_image(df):
file = np.random.choice(df.fname)
im = skimage.io.imread(data_folder+file)
annots = df[df.fname==file].iterrows() plt.figure(figsize=(6,6))
plt.imshow(im) current_axis = plt.gca() for box in annots:
label = box[1]['class_id']
current_axis.add_patch(plt.Rectangle(
(box[1]['xmin'], box[1]['ymin']), box[1]['xmax']-box[1]['xmin'],
box[1]['ymax']-box[1]['ymin'], color='blue', fill=False, linewidth=2))
current_axis.text(box[1]['xmin'], box[1]['ymin'], label, size='x-large', color='white', bbox={'facecolor':'blue', 'alpha':1.0})
plt.show() viz_random_image(df)
A representative sample form SVHN dataset
对于接下来的步骤，我在 repo 中提供了一个 utils_ssd.py ，以方便训练、负重等。部分代码取自 SSD_Keras repo，也是广泛使用的。
第三步:选择策略
如前所述，对于这个问题，我们有许多可能的方法。在本教程中，我将采用标准的深度学习检测方法，并将使用 SSD 检测模型。我们将从这里的使用 SSD keras 实现。这是 PierreLuigi 很好的实现。虽然它的 GitHub 明星比 rykov8 的实现少，但它看起来更新了，也更容易集成。当你选择使用哪个项目时，这是一件非常重要的事情。其他好的选择将是 YOLO 模型和掩模 RCNN。
步骤 4:加载并训练 SSD 模型
一些定义
要使用 repo，您需要验证您有 SSD_keras repo，并填写 json_config.json 文件中的路径，以允许笔记本找到这些路径。
从导入开始:
**import** **os**
**import** **sys**
**import** **skimage.io**
**import** **scipy**
**import** **json****with** open('json_config.json') **as** f: json_conf = json.load(f)ROOT_DIR = os.path.abspath(json_conf['ssd_folder']) *# add here mask RCNN path*
sys.path.append(ROOT_DIR)
**import** **cv2**
**from** **utils_ssd** **import** *
**import** **pandas** **as** **pd**
**from** **PIL** **import** Image
**from** **matplotlib** **import** pyplot **as** plt
%matplotlib inline
%load_ext autoreload
% autoreload 2
以及更多的定义:
task = 'svhn'labels_path = f'{data_folder}pascal.csv'input_format = ['class_id','image_name','xmax','xmin','ymax','ymin' ]
df = pd.read_csv(labels_path)
车型配置:
**class** **SVHN_Config**(Config):
batch_size = 8
dataset_folder = data_folder
task = task
labels_path = labels_path
input_format = input_format
conf=SVHN_Config()
resize = Resize(height=conf.img_height, width=conf.img_width)
trans = [resize]
定义型号，装载重量
与大多数深度学习案例一样，我们不会从头开始训练，但我们会加载预训练的权重。在这种情况下，我们将加载 SSD 模型的权重，在 COCO 数据集上训练，该数据集有 80 个类。很明显，我们的任务只有 10 个类，因此在加载权重后，我们将重建顶层以获得正确数量的输出。我们在 init_weights 函数中完成。补充说明:在这种情况下，正确的输出数量是 44:每个类 4 个(边界框坐标),另外 4 个用于背景/无类。
learner = SSD_finetune(conf)
learner.get_data(create_subset=**True**)
weights_destination_path=learner.init_weights()
learner.get_model(mode='training', weights_path = weights_destination_path)
model = learner.model
learner.get_input_encoder()
ssd_input_encoder = learner.ssd_input_encoder
*# Training schedule definitions*
adam = Adam(lr=0.0002, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0)
ssd_loss = SSDLoss(neg_pos_ratio=3, n_neg_min=0, alpha=1.0)
model.compile(optimizer=adam, loss=ssd_loss.compute_loss)
定义数据加载器
train_annotation_file=f'**{conf.dataset_folder}**train_pascal.csv'
val_annotation_file=f'**{conf.dataset_folder}**val_pascal.csv'
subset_annotation_file=f'**{conf.dataset_folder}**small_pascal.csv'batch_size=4
ret_5_elements={'original_images','processed_images','processed_labels','filenames','inverse_transform'}train_generator = learner.get_generator(batch_size, trans=trans, anot_file=train_annotation_file,
encoder=ssd_input_encoder)val_generator = learner.get_generator(batch_size,trans=trans, anot_file=val_annotation_file,
returns={'processed_images','encoded_labels'}, encoder=ssd_input_encoder,val=True)
5.训练模型
现在模型准备好了，我们将设置一些与上次训练相关的定义，并开始训练
learner.init_training()history = learner.train(train_generator, val_generator, steps=100,epochs=80)
作为奖励，我在训练脚本中包含了 training_plot 回调，以在每个时期后可视化一个随机图像。例如，下面是第六个纪元后的预测快照:
SSD_Keras repo 在几乎每个时期后处理保存模型，因此您可以稍后加载模型，只需更改 weights_destination_path 行使其等于路径
weights_destination_path = <path>
如果你按照我的指示做，你应该能训练出这个模型。ssd_keras 提供了更多的特性，例如数据扩充、不同的加载器和评估器。经过短暂的训练，我已经达到了> 80 地图。
你达到了多高？
Training for 4X100X60 samples, from tensorboard
摘要
在这篇文章中，我们讨论了 OCR 领域的不同挑战和方法。正如深度学习/计算机视觉中的许多问题一样，它比乍看起来要复杂得多。我们已经看到了它的许多子任务，以及解决它的一些不同的方法，目前没有一个是银弹。另一方面，我们已经看到，没有太多的争论，达成初步结果并不困难。
希望你喜欢！
扫一扫
722
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
