TowardsDataScience 博客中文翻译 2016~2018（一百六十九）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

自然语言处理导论

原文：https://towardsdatascience.com/introduction-to-natural-language-processing-part-1-777f972cc7b3?source=collection_archive---------6-----------------------

用数学破译语言

Photo by Franki Chamaki on Unsplash

自然语言处理(NLP)是人工智能领域，涉及人类语言的处理和理解。自 20 世纪 50 年代问世以来，机器对语言的理解已经在翻译、主题建模、文档索引、信息检索和实体抽取中发挥了关键作用。

如今，它被用来驱动搜索引擎，过滤垃圾邮件，并以快速和可扩展的方式获得分析。随着计算效率和机器学习的新发展，这些系统的性能呈指数增长。今天的几个 NLP 系统甚至可以夸耀接近人类水平的对等性。

有了专门用于 NLP 的剩余工具和技术，现在是所以很容易开始。我的目标是涵盖自然语言处理的基本理论，并向您展示如何构建简单的工具来展示它们的魔力。查看这些后续文章，获取更多 NLP 阅读资料。

• 垃圾邮件还是火腿——一个简单的文本分类问题
• 探索单词嵌入的神经网络

顺便提一下:这篇文章以及后续的文章主要关注文本处理环境中的 NLP。

我们到底是如何对语言建模的？

语言是流动的。语言充满了抽象和歧义。成年后学习第二语言很难，那么机器到底是怎么做到的呢？

事实是，大多数自然语言理解系统都是非常结构化的。他们将文本数据转换成数字表示，并使用复杂的数学模型来学习每个单词、短语和文档之间的依赖关系。最终用户所看到的对语言的智能理解实际上是一个大型语料库中的数百万个参数。

把文字变成数字的过程俗称矢量化或者嵌入。这些技术是将单词映射到实数向量的函数。向量形成了向量空间，这是一个代数模型，其中应用了向量加法的所有规则和相似性度量。

Using the word embedding technique word2vec, researchers at Google are able to quantify word relationships in an algebraic model

向量空间中两个单词之间的距离度量了它们之间的某种相似性。还可以使用向量加法来量化单词关系。上图展示了国王和王后的关系就像男人和女人的关系一样。理论上，这将意味着King Vector — Queen Vector = Man Vector - Woman Vector

每个向量的粒度不限于单词，我们也可以将字符、短语或文档映射到向量

“文档”是指与问题相关的完整的文本实体。例如，垃圾邮件分类系统会将每封邮件视为一个文档。

文档可以映射到向量的向量，其中每个向量对应一个单词。在这种情况下，向量的向量被算法摄取，其中单词关系和单个向量的连接被用于导出含义。或者，可以将整个文档映射到一个向量。这创建了对语言的更简化的理解，但是更少依赖于数据量和计算资源。

文档级嵌入通常通过检查每个单词的频率来创建。这些嵌入技术依赖于分布假设，假设:

分布假设:在相同语境中使用和出现的词倾向于表达相似的意思

让我们来看看最简单的嵌入技术之一——单词包，它用于将每个文档映射到它自己的向量上。

一袋单词

由于语言过于简单，这种技术在实践中并不常用。它的使用案例主要是作为一种教育工具，旨在使 NLP 的学生轻松地进入该主题。

让我们考虑以下文件。如果有帮助的话，你可以想象它们是朋友之间分享的荒谬简单的短信。

文件一:击掌

文件二:我老了。

文件 3:她五岁了。

我们从这组文档中获得的词汇是(高，五，我，我，老，她，是)。我们现在将忽略标点符号，尽管根据我们的使用情况，将它们纳入我们的词汇表也很有意义。

我们可以创建一个矩阵来表示词汇表中的每个术语和文档之间的关系。矩阵中的每个元素表示该术语在特定文档中出现的次数。

使用这个矩阵，我们可以获得每个单词以及文档的向量。我们可以将“五”矢量化为[1，0，1]，将“文档 2”矢量化为[0，0，1，1，1，0，0]。

单词袋并不能很好地代表语言，尤其是当你的词汇量很小的时候。它忽略了词序和词的关系，产生了稀疏向量，其中大部分是零。从这个小例子中我们还可以看到，单词“我”、“am”、“old”被映射到同一个向量。这暗示这几个字差不多，其实没什么意义。

这里有一些关于 NLP 阅读的文章。 Spam 或 Ham 引入了一种新的文档嵌入类型，它考虑了每个术语的相对重要性。探索单词嵌入的神经网络介绍了使用神经网络创建的嵌入，这两种方法都是单词袋的更实用的替代方法。

• 垃圾邮件还是火腿—一个简单的文本分类问题
• 探索单词嵌入的神经网络

在我矢量化/嵌入我的文本后，现在该怎么办？

数字表示允许我们使用数学模型来分析我们的文本。下一步是将这些嵌入作为分类问题中的特征。

By one hot encoding the classes, we can plug this data into any type of classifier

文本数据丰富且维数高，这使得更复杂的分类器(如神经网络)成为 NLP 的理想选择。

感谢您的阅读！

如果你喜欢这篇文章，可以看看我关于数据科学、数学和编程的其他文章。通过 Medium 关注我的最新更新。😃

作为一个业余爱好项目，我还在www.dscrashcourse.com建立了一套全面的免费数据科学课程和练习题。

如果你想支持我的写作，下次你报名参加 Coursera 课程时，可以考虑使用我的会员链接。完全公开—我从每一次注册中获得佣金，但不会对您产生额外费用。

再次感谢您的阅读！📕

神经网络导论

原文：https://towardsdatascience.com/introduction-to-neural-networks-ead8ec1dc4dd?source=collection_archive---------2-----------------------

神经网络基础

什么是神经网络

神经网络，通常称为人工神经网络(ANN ),是机器学习(ML)问题中对人脑功能的模拟。人工神经网络不能解决所有出现的问题，但可以和其他技术一起为各种 ML 任务提供更好的结果。人工神经网络最常见的用途是聚类和分类，它们也可用于回归任务，但在这方面还有更好的方法。

人工神经网络的构建模块和功能

神经元

这是神经网络的构建单元，它模仿人类神经元的功能。典型的神经网络使用 sigmoid 函数，如下所示。使用这个函数主要是因为它能够根据 f(x) 本身写出导数，这在最小化误差时很方便。

sigmoid function

Neuron

*z* =  ∑ *w*×*x
y = sigmoid(z)**w* = weights
*x* = inputs

神经元按层连接，因此一层可以与其他层通信，形成神经网络。除了输入和输出层之外的内层称为隐藏层。一层的输出被馈送到另一层的输入。

调整重量

人工神经网络的学习任务是调整权重以最小化误差。这是通过误差反向传播实现的。对于使用 Sigmoid 函数作为激活函数的简单神经元，误差可以如下所示。让我们考虑一种一般情况，其中权重称为向量 W，输入称为向量 x。

Error calculation and weight adjustment

从上面的等式中，我们可以概括出权重调整，令人惊讶的是，您会注意到这仅需要相邻神经元级别的细节。因此，这是一种稳健的学习机制，称为反向传播算法。从输出节点开始反向传播，更新前一个神经元的权重。

带有简单 javascript 库的演示应用程序

让我们编写一个简单的应用程序，它将使用两个图像进行训练，并对给定的图像应用过滤器。以下是训练过程的源图像和目标图像。

Source image(left) and target image(right)

我使用了一个人工神经网络，它使用反向传播来调整误差。训练的意图是找到一个函数 *f(红、绿、蓝、阿尔法)*来匹配目标颜色变换。使用源图像的几种颜色调整来制作目标图像。让我们看看代码。

Package.json

Index.ts

$ npm install
$ npm start

源图像应为input_image_train.jpg，目标图像名称应为output_image_train.jpg。要应用滤镜的图像文件应该是test.jpg，新的图像文件将保存为out.jpg。以下是我用训练好的模型过滤的一些示例图像。

样本输出

Input images (left) and output images(right)

酷吧？训练需要几秒钟，但过滤是即时的。如果需要，您可以保存模型以供将来使用。这在现实世界的应用中是非常智能的。

用于数据分析的 Numpy 简介

原文：https://towardsdatascience.com/introduction-to-numpy-for-data-analysis-3aaded0d0996?source=collection_archive---------3-----------------------

查看更新的 DevOps 课程。

课程注册链接:

[## 创建一个帐户

编辑描述

www.101daysofdevops.com](https://www.101daysofdevops.com/register/)

课程链接:

[## 101 天的开发周期

http://100daysofdevops.com/day-100-100-days-of-devops/好消息是从话题开始包括(请让我知道…

www.101daysofdevops.com](https://www.101daysofdevops.com/courses/101-days-of-devops/)

YouTube 链接:

[## 普拉尚·拉赫拉

主要想法是通过视频分享我的 Linux 知识。我的 Linkedin 个人资料…

www.youtube.com](https://www.youtube.com/user/laprashant/videos)

NumPy 是 Python 的线性代数库，这也是 PyData 生态系统中所有库都依赖 NumPy 作为主要构建模块的重要原因。

安装 Numpy

*# pip2 install numpy**Collecting numpy**Using cached numpy-1.12.1-cp27-cp27mu-manylinux1_x86_64.whl**Installing collected packages: numpy**Successfully installed numpy-1.12.1*

强烈建议使用 Anaconda 发行版安装 Python，以确保所有底层依赖项( 如线性代数库 )都与 conda 安装同步。

如果你有康达安装【https://www.continuum.io/downloads】T21

*conda install numpy*

Numpy 数组是我们使用 Numpy 的主要原因，它们有两种类型

• 向量(一维数组)
• 矩阵(二维数组)

*# 1-D Array**>>> test = [1,2,3]**>>> import numpy as np**# We got the array
>>> np.array(test)**array([1, 2, 3])**>>> arr = np.array(test)*

我们来看看 二维数组

*>>> test1 = [[1,2,3],[4,5,6],[7,8,9]]**>>> test1**[[1, 2, 3], [4, 5, 6], [7, 8, 9]]**>>> np.array(test1)**array([[1, 2, 3],**[4, 5, 6],**[7, 8, 9]])*

但是生成 NumPy 数组最常见的方法是使用range函数(类似于 Python 中的 range)

*#Similar to range(start,stop,step),stop not included and indexing start with zero 
>>> np.arange(0,10)**array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])*

但是如果我们在寻找特定类型的数组

*>>> np.zeros(3)**array([ 0., 0., 0.])**#We are passing Tuple where* ***first value represent row*** *and* ***second represent column****>>> np.zeros((3,2))**array([[ 0., 0.],**[ 0., 0.],**[ 0., 0.]])*

类似的还有那些

*>>> np.ones(2)**array([ 1., 1.])**>>> np.ones((2,2))**array([[ 1., 1.],**[ 1., 1.]])*

现在我们来看看linspace

*#It will give 9 evenly spaced point between 0 and 3(It return 1D vector)
>>> np.linspace(0,3,9)**array([ 0\. , 0.375, 0.75 , 1.125, 1.5 , 1.875, 2.25 , 2.625, 3\. ])**>>> np.linspace(0,10,3)**array([  0.,   5.,  10.])*

让我们创建 单位矩阵 (二维方阵，其中行数等于列数，对角线为 1)

创建数组 随机数

*#1-D, it create random sample uniformly distributed between 0 to 1
>>> np.random.rand(3)**array([ 0.87169008, 0.51446765, 0.65027072])**#2-D
>>> np.random.rand(3,3)**array([[ 0.4217015 , 0.86314141, 0.14976093],**[ 0.4348433 , 0.68860693, 0.88575823],**[ 0.56613179, 0.56030069, 0.51783999]])*

现在如果我想要 随机整数

*#This will give random integer between 1 and 50
>>> np.random.randint(1,50)**27**#In case if we need 10 random integer
>>> np.random.randint(1,50,10)**array([39, 34, 30, 21, 18, 30,  3,  6, 37, 11])*

我们可以 重塑 我们现有的阵列

*>>> np.arange(25)**array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,**17, 18, 19, 20, 21, 22, 23, 24])**>>> arr = np.arange(25)**>>> arr**array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,**17, 18, 19, 20, 21, 22, 23, 24])**>>> arr.reshape(5,5)**array([[ 0, 1, 2, 3, 4],**[ 5, 6, 7, 8, 9],**[10, 11, 12, 13, 14],**[15, 16, 17, 18, 19],**[20, 21, 22, 23, 24]])*

让我们来看看其他一些方法

*>>> np.random.randint(0,50,10)**array([10, 40, 18, 30, 6, 40, 49, 23, 3, 18])**>>> ranint = np.random.randint(0,50,10)**>>> ranint**array([18, 49, 6, 28, 30, 10, 46, 11, 40, 16])**#It will return* ***max value*** *of the array
>>> ranint.max()**49****#Minimum value*** *>>> ranint.min()**6**>>> ranint**array([18, 49,  6, 28, 30, 10, 46, 11, 40, 16])**#To find out the position
>>> ranint.argmin()**2**>>> ranint.argmax()**1*

找出一个数组的形状

*>>> arr.shape**(25,)**>>> arr.reshape(5,5)**array([[ 0, 1, 2, 3, 4],**[ 5, 6, 7, 8, 9],**[10, 11, 12, 13, 14],**[15, 16, 17, 18, 19],**[20, 21, 22, 23, 24]])**>>> arr = arr.reshape(5,5)**>>> arr.shape**(5, 5)*

找出 的数据类型

*>>> arr.dtype**dtype(‘int64’)*

NumPy 情况下的分度

*>>> import numpy as np**>>> arr =np.arange(0,11)**>>> arr**array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])**>>> arr[0]**0**>>> arr[0:4]**array([0, 1, 2, 3])*

Numpy 数组与 Python 列表的不同之处在于它们的广播能力

*>>> arr[:] = 20**>>> arr**array([20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20])**# Now let's try to slice this array**>>> arr1 = arr[0:5]**>>> arr1**array([0, 1, 2, 3, 4])**>>> arr1[:] = 50**>>> arr1**array([50, 50, 50, 50, 50])**#But as you can see the side effect it change the original array too(i.e data is not copied it's just the view of original array)* ***>>> arr******array([50, 50, 50, 50, 50,  5,  6,  7,  8,  9, 10])****#If we want to avoid this feature, we can copy the array and then perform broadcast on the top of it**>>> arr2 = arr.copy()**>>> arr2**array([50, 50, 50, 50, 50,  5,  6,  7,  8,  9, 10])**>>> arr2[6:10] = 100**>>> arr2**array([ 50,  50,  50,  50,  50,   5, 100, 100, 100, 100,  10])**>>> arr**array([50, 50, 50, 50, 50,  5,  6,  7,  8,  9, 10])*

索引 二维数组(矩阵)

*>>> arr = ([1,2,3],[4,5,6],[7,8,9])
>>> arr**([1, 2, 3], [4, 5, 6], [7, 8, 9])**>>> arr1 = np.array(arr)**>>> arr1**array([[1, 2, 3],**[4, 5, 6],**[7, 8, 9]])**>>> arr[1]**[4, 5, 6]**# To grab 5(Indexing Start with zero)**>>> arr1[1][1]**5**#Much shortcut method
>>> arr1[1,1]**5*

从二维数组中抓取元素

*>>> arr**array([[1, 2, 3],**[4, 5, 6],**[7, 8, 9]])**#This will grab everything from Row 1 except last element(2) and staring from element 1 upto the end from Row 2* *>>> arr[:2,1:]**array([[2, 3],**[5, 6]])*

条件选择 :这样会返回布尔值

*>>> arr**array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])**>>> arr > 5**array([False, False, False, False, False, True, True, True, True, True], dtype=bool)**# We can save this value to an array and perform boolean selection**>>> my_arr = arr > 5**>>> my_arr**array([False, False, False, False, False,  True,  True,  True,  True,  True], dtype=bool)**>>> arr[my_arr]**array([ 6,  7,  8,  9, 10])****#OR much easier way****>>> arr[arr > 5]**array([ 6,  7,  8,  9, 10])**>>> arr[arr < 5]**array([1, 2, 3, 4])*

操作

*# It's the same operation as we are doing with Normal Python**>>> arr = np.arange(0,10)**>>> arr**array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])**>>> arr**array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])****#Addition*** *>>> arr + arr**array([ 0, 2, 4, 6, 8, 10, 12, 14, 16, 18])****#Substraction*** *>>> arr — arr**array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])****#Multiplication*** *>>> arr * arr**array([ 0, 1, 4, 9, 16, 25, 36, 49, 64, 81])****#Broadcast(It add's/substract/multiply 100 to each element)*** *>>> arr + 100**array([100, 101, 102, 103, 104, 105, 106, 107, 108, 109])**>>> arr - 100**array([-100,  -99,  -98,  -97,  -96,  -95,  -94,  -93,  -92,  -91])**>>> arr * 100**array([  0, 100, 200, 300, 400, 500, 600, 700, 800, 900])*

在 Python 的情况下，如果我们试图用零除一，我们将得到被零除的异常

*>>> 0/0**Traceback (most recent call last):**File "<stdin>", line 1, in <module>**ZeroDivisionError: division by zero****OR****>>> 1/0**Traceback (most recent call last):**File “<stdin>”, line 1, in <module>**ZeroDivisionError: division by zero*

在 Numpy 的情况下，如果我们试图除以零，我们不会得到任何异常，但它会返回nan*(不是一个数字)*

*#Not giving you 
>>> arr/arr**__main__:1: RuntimeWarning: invalid value encountered in true_divide**array([* ***nan****, 1., 1., 1., 1., 1., 1., 1., 1., 1.])*

如果被零除，它将返回 无穷大

*>>> 1/arr**array([* ***inf****, 1\. , 0.5 , 0.33333333, 0.25 ,**0.2 , 0.16666667, 0.14285714, 0.125 , 0.11111111])*

万能数组功能

*#Square root
>>> np.sqrt(arr)**array([ 0\. , 1\. , 1.41421356, 1.73205081, 2\. ,**2.23606798, 2.44948974, 2.64575131, 2.82842712, 3\. ])**#Exponential
>>> np.exp(arr)**array([ 1.00000000e+00, 2.71828183e+00, 7.38905610e+00,**2.00855369e+01, 5.45981500e+01, 1.48413159e+02,**4.03428793e+02, 1.09663316e+03, 2.98095799e+03,**8.10308393e+03])**#Maximum
>>> np.max(arr)**9**#Minimum
>>> np.min(arr)**0**#Logarithmic
>>> np.log(arr)**__main__:1: RuntimeWarning: divide by zero encountered in log**array([       -inf,  0\.        ,  0.69314718,  1.09861229,  1.38629436,**1.60943791,  1.79175947,  1.94591015,  2.07944154,  2.19722458])*

遗传算法优化导论

原文：https://towardsdatascience.com/introduction-to-optimization-with-genetic-algorithm-2f5001d9964b?source=collection_archive---------1-----------------------

为机器学习任务选择最佳参数是具有挑战性的。一些结果可能不是因为数据有噪声或者所使用的学习算法弱，而是因为参数值的错误选择。本文简要介绍了进化算法，并描述了一种最简单的基于随机的进化算法&遗传算法。

From Pixabay by qimono

简介

假设一位数据科学家将一个图像数据集划分为多个类别，并且要创建一个图像分类器。在数据科学家研究了数据集之后，K-最近邻(KNN)似乎是一个不错的选择。为了使用 KNN 算法，有一个要使用的重要参数是 k。假设选择了初始值 3。科学家从选定的 K=3 开始 KNN 算法的学习过程。生成的训练模型达到了 85%的分类准确率。这个百分比可以接受吗？换句话说，我们能得到比目前更好的分类精度吗？在进行不同的实验之前，我们不能说 85%是最佳精度。但是要做另一个实验，我们必须改变实验中的一些东西，比如改变 KNN 算法中使用的 K 值。我们不能肯定地说 3 是在这个实验中使用的最佳值，除非尝试对 K 应用不同的值并注意分类准确度如何变化。问题是“如何找到使分类性能最大化的 K 的最佳值？”这就是所谓的优化。

在最优化中，我们从实验中使用的变量的某种初始值开始。因为这些值可能不是最好的，所以我们应该改变它们，直到得到最好的值。在某些情况下，这些值是由复杂的函数生成的，我们无法轻松地手动求解。但是进行优化是非常重要的，因为分类器可能产生差的分类精度，这不是因为例如数据有噪声或者所使用的学习算法弱，而是由于学习参数初始值的差的选择。因此，运筹学(OR)研究人员提出了不同的优化技术来完成这种优化工作。根据[1]，优化技术分为四个主要类别:

1. 约束优化
2. 多模态优化
3. 多目标优化
4. 组合最优化

观察各种自然物种，我们可以注意到它们是如何进化和适应环境的。我们可以从这种已经存在的自然系统及其自然进化中受益，来创造我们的人工系统做同样的工作。这叫仿生学。例如，飞机是基于鸟类如何飞行，雷达来自蝙蝠，潜艇是基于鱼类发明的，等等。因此，一些优化算法的原理来自自然。比如遗传算法(GA)，其核心思想来自查尔斯·达尔文的自然进化论“适者生存”。

在进入遗传算法工作的细节之前，我们可以对进化算法有一个总体的了解。

进化算法

我们可以说，优化是使用进化算法(EAs)执行的。传统算法和进化算法的区别在于进化算法不是静态的，而是动态的，因为它们可以随着时间的推移而进化。

进化算法有三个主要特征:

1. 基于群体的:进化算法是优化一个过程，在这个过程中，当前的解决方案是不好的，以产生新的更好的解决方案。从中产生新解的当前解的集合称为群体。
2. 健身导向:如果有几个方案，怎么说一个方案比另一个方案好？从适应度函数计算出的每个单独的解决方案都有一个相关的适应度值。这样的适应值反映了解决方案有多好。
3. 变异驱动:如果根据每个个体计算出的适应度函数，在当前种群中没有可接受的解，我们就要做出一些东西来生成新的更好的解。因此，单个解决方案将经历许多变化以生成新的解决方案。

我们将转到正式发布并应用这些条款。

遗传算法

遗传算法是一种基于随机的经典进化算法。这里的随机是指为了使用 GA 找到一个解决方案，对当前解决方案进行随机更改以生成新的解决方案。请注意，遗传算法可以被称为简单遗传算法(SGA)，因为它比其他进化算法简单。

遗传算法是基于达尔文的进化论。这是一个缓慢渐进的过程，通过做出微小而缓慢的改变来实现。此外，遗传算法缓慢地对其解进行细微的改变，直到得到最佳解。

下面是遗传算法工作原理的描述:

GA 处理由一些解组成的群体，其中群体大小(popsize)是解的数量。每个解决方案都被称为个体。每个溶液都有一条染色体。染色体被表示为定义个体的一组参数(特征)。每条染色体都有一组基因。每个基因都以某种方式表示，比如表示为一串 0 和 1，如图 1 所示。

Figure 1. Chromosome and gene.

而且，每个个体都有一个适合度值。为了选择最佳个体，使用适应度函数。适应度函数的结果是代表解的质量的适应度值。适应值越高，解决方案的质量越高。基于质量的最佳个体的选择被用于产生所谓的交配池，其中质量较高的个体在交配池中被选择的概率较高。

交配池中的个体被称为父母。从交配池中选出的每两个父母将产生两个后代(孩子)。仅仅通过优质个体的交配，就有望获得比其父母更优质的后代。这将杀死不良个体，防止产生更多的不良个体。通过不断选择和交配高质量的个体，将有更高的机会只保留个体的良好特性，而忽略不好的特性。最后，这将以期望的最优或可接受的解决方案结束。

但是目前使用所选择的父母产生的后代仅仅具有其父母的特征，而没有更多变化。没有新的东西加入其中，因此它的双亲中的同样的缺点实际上会存在于新的后代中。为了克服这样的问题，一些变化将被应用于每个后代以产生新的个体。所有新生成的个体的集合将是替换先前使用的旧群体的新群体。每创造一个群体称为一代。新人口取代旧人口的过程称为更替。图 2 总结了 GA 的步骤。

Figure 2. Genetic algorithm steps.

要全面了解 GA，需要回答两个问题:

1. 两个后代是如何从双亲中产生的？
2. 每一个后代是如何稍微改变成为一个个体的？

这些问题我们以后再来回答。

染色体表示和评估

染色体有不同的表示法，选择正确的表示法要视具体问题而定。好的表示使搜索空间更小，从而更容易搜索。

可用于染色体的表示包括:

二进制:每个染色体被表示为一串 0 和 1。

排列:适用于排序问题，如旅行商问题。

值:实际值按原样编码。

例如，如果我们用二进制编码数字 7，它可能如图 3 所示。

Figure 3. Binary encoding example.

上述染色体的每一部分被称为基因。每个基因都有两种特性。第一个是它的值(等位基因),第二个是在染色体内的位置(基因座),它是高于它的值的数字。

每条染色体有两种表现形式。

1. 基因型:代表染色体的一组基因。
2. 表现型:染色体的实际物理表现。

在上面的例子中，0111 的二进制是基因型，7 是表型表示。

在用正确的方式表示每个染色体后，用来搜索空间，接下来是计算每个个体的适应值。假设我们示例中使用的适应度函数为:

f(x)=2x+2

其中 x 是染色体值

那么前一个染色体的适应值是:

f(7)=2(7)+2=16

计算染色体适应值的过程称为评估。

初始化

在获得如何表示每个个体之后，下一步是通过在群体中选择适当数量的个体来初始化群体。

选择

下一步是从交配池的种群中选择一些个体。基于先前计算的适应值，基于阈值的最佳个体被选择。在这一步之后，我们将结束选择交配池中种群的子集。

变异算子

基于交配池中选择的个体，选择亲本进行交配。每两个亲本的选择可以通过依次选择亲本(1-2、3-4 等)来进行。另一种方法是随机选择父母。

对于每两个选定的父代，都有许多变异运算符要应用，例如:

1. 交叉(重组)
2. 变化

图 4 给出了这些操作符的一个例子。

Figure 4. Crossover and mutation.

交叉

遗传算法中的交叉与自然变异一样产生新的一代。通过改变老一代父母，新一代的后代携带父母双方的基因。每个父母携带的基因数量是随机的。记住 GA 是基于随机的 EA。有时后代从父母一方获得一半基因，从另一方获得另一半基因，有时这种百分比会发生变化。对于每两个亲本，通过选择染色体中的随机点并交换亲本该点前后的基因来进行交叉。产生的染色体就是后代。因此算子被称为单点交叉。

请注意，交叉很重要，没有它，后代将与其父代相同。

突变

下一个变异算子是变异。对于每个后代，选择一些基因并改变其值。突变因染色体表现形式而异，但如何应用突变取决于您。如果编码是二进制的(即每个基因的值空间只有两个值 0 和 1)，则翻转一个或多个基因的位值。

但如果基因值来自 1、2、3、4、5 等两个以上值的空间，那么二元突变就不适用，我们应该另寻他法。一种方法是从如图 5 所示的一组值中选择一个随机值。

Figure 5. Mutation by randomly updating some genes.

请注意，如果没有突变，后代将具有其父母的所有特性。为了给这样的后代增加新的特征，基因发生了突变。但是因为突变是随机发生的，所以不建议增加应用于突变的基因数量。

突变后的个体称为突变体。

参考

[1]艾本、阿戈斯顿和詹姆斯·史密斯。进化计算导论。第 53 卷。海德堡:施普林格，2003 年。

原文在 LinkedIn 可在此链接:https://www . LinkedIn . com/pulse/introduction-optimization-genetic-algorithm-Ahmed-gad

在 KDnuggets 上也分享到了这个链接:https://www . kdnugges . com/2018/03/introduction-optimization-with-genetic-algorithm . html

联系作者:

领英:【https://linkedin.com/in/ahmedfgad

电子邮件:ahmed.f.gad@gmail.com

介绍我们的放射学/人工智能系列

原文：https://towardsdatascience.com/introduction-to-our-radiology-ai-series-51d0e54cd3b9?source=collection_archive---------17-----------------------

Photo by rawpixel on Unsplash

放射科医生的机器学习介绍— 10 部分系列

第一部

作者: 达尼洛泰·瓦雄

编辑: 迈克尔博士

在 2017 年，Kaggle 数据科学碗旨在使用机器学习和人工智能来抗击美国男性和女性癌症死亡的主要原因。参赛者被要求使用成千上万张高分辨率肺部 CT 图像的数据集来创建新的肺癌检测算法。这些算法是为了提高诊断和减少假阳性率。

在 394 个参赛队中，哪个队获得了头奖？一个由中国清华大学医学和计算机科学系成员组成的团队。

诸如此类的比赛是将国际人才与全球问题结合起来的好方法。这种团队合作的风格仅仅是通过医学专业人员和计算机科学之间的互动，在我们的领域内的无限发展潜力的皮毛。

在放射学培训中，我们了解到一个 3 厘米长、呈毛刺状、软组织变薄的肺部肿块很有可能是癌症。同样，一个 5 毫米的光滑钙化结节患癌的概率也很低。

然而，我们也知道许多肺结节介于我们准确预测恶性肿瘤的能力之间。弗莱舍纳协会非常努力地在 2017 年更新了后续标准，提供了一个解决方案，其中包括大小和密度的变化。然而，我们仍然不能看着一个 8 毫米的边缘稍微不规则的结节就说它有多大可能是癌症。

为了使 Kaggle 竞赛更进一步，Fleischner 标准(或其替代标准)非常有可能是非常可定制的，肺结节跟踪将得到改善。我们将在第 2 篇文章中对此进行更深入的探讨。

这些团队在脸书的面部识别背后使用了类似的技术。你有没有想过他们怎么能确定你照片里的人是谁？这项技术被称为深度学习，是机器学习的子集，也是人工智能的子集，我们也会在第 3、5、6 篇文章中更深入地探讨这些话题。

这些努力证明了开源社区，以及人们如何决心通过合作和共享数据来找到重要问题的新颖解决方案。

让我们采取相反的观点，机器学习和人工智能可能预示着放射学专业的过时。我想起了杰弗里·辛顿教授。

杰弗里·辛顿是个非常聪明的家伙，但他缺乏关于放射学专业的细微差别的医学训练，如图像引导活检、肿瘤板和与我们外科同事的讨论等。，也许让他对我们的职业只有一个肤浅的看法。放射科医生的工作将随着新工具的出现而改变，但只要我们继续帮助我们的临床同事，他们的工作就会一直存在。

Kaggle 竞赛的结果也得益于时间和可用资源。由于巨大的视频游戏市场，我们第一次有了经济高效的高性能计算，称为图形处理单元(GPU)。详见第三条。另一个我们可能想当然的领域是语音识别。虽然我们可能看到也可能看不到每天的改善，但在过去的十年里，情况确实有所改善。详见第 4 条。

这个难题的另一个非常重要的部分是数据。大量数据。大量数据被正确标记。美国放射学院(ACR)和斯坦福大学目前正在进行这方面的研究。详见第五条。(或者 7)。

计算机科学家和医学专家的合作团队在开发改变领域的算法方面有着惊人的潜力。但是，当我们谈论将这些技术融入我们的日常工作流程，或者在隐私或数据管理的背景下，请注意。详见第八条。

在我们继续之前，我们想正式介绍一下我们自己。

泰·瓦雄

这是旨在指导我的同事的 10 篇系列文章中的第一篇。自 2012 年以来，我一直在跟踪人工智能和医疗应用的发展，并追随德雷耶博士和米哈尔斯基博士等伟大的导师。信息学，尤其是医学图像的利用，是我背景的很大一部分。在海军陆战队担任飞行外科医生后，我在美国海军接受了放射学培训。我的最后一次海军之旅是在冲绳，包括一次放射科主任之旅，之后我完成了我的海军使命，并作为天使投资人、企业家、信息学顾问和咨询顾问返回圣地亚哥。在撰写本文时，我还没有关于这个系列的相关财务披露。”

达尼洛·佩纳

“我有化学工程背景，做过两年工程师。在工作期间，我意识到我需要对社会产生更大的影响，我也想学习编码。因此，我申请了学校，被录取了，然后辞去了工作。我目前是休斯顿德克萨斯大学健康科学中心的生物医学信息学硕士生，也是艾伯特-史怀哲研究员。我总是在学习，我很高兴能帮助别人学习我所知道的。我希望通过这一系列文章，从医疗领域到机器学习领域的人们，乃至普通人，都可以利用这些信息来了解放射学的当前格局及其与人工智能技术进步的关系。”

我们相信，通过我们不同但互补的技能组合，我们可以在这个令人兴奋的领域教育他人。

既然你对我们有所了解，你可能会问为什么你要花时间在我们的系列节目上。

在本系列中，我们将从慢开始，回顾关键术语。我们还将讨论足够多的近期历史进展，以提供背景和应对新趋势。如果你更高级一点，请在评论中提供来自个人经历的清晰想法。当然，如果你注意到任何错误的文字，我们也会虚心审查这些评论。

本系列并不打算面面俱到，但这些文章旨在为放射学和人工智能提供公平的竞争环境。这个领域变化很快，有很多移动的部分。我们正在尽自己的一份力量通过这一过程进行教育和学习。

加入我们一组有趣的文章，这些文章由一位经验丰富的放射科医生策划，他专注于技术，是一名对理解 ML 和 AI 将如何影响下一代医疗保健感兴趣的学生。

第二部分——放射学和人工智能将如何走到一起

语法演变的 PonyGE2 介绍

原文：https://towardsdatascience.com/introduction-to-ponyge2-for-grammatical-evolution-d51c29f2315a?source=collection_archive---------1-----------------------

语法进化是算法优化的一种强有力的方法。给定一个目标函数和一个搜索空间(语法)，有可能使用进化计算方法来进化一个算法，以最优地(或至少有效地)最大化目标函数。

PonyGE2 是 Python3 的一个实现，它让 GE 变得简单。本演练讨论了 GE 的优点和缺点，并提供了一个如何用 GE 改进您自己的算法的示例。一篇全面描述 PonyGE2 运作的论文可以在这里找到。

不幸的是，在当前计算机上用 PonyGE2 进行语法进化还不够快，不能在合理的时间内从大量语法中进化出函数，但这种技术肯定是未来需要关注的。

语法演变

这种方法最初是由的迈克尔·奥尼尔(付费墙)博士在 2001 年提出的，作为一种将进化计算方法应用于具有正式语法结构的问题的方式，比如算法。

一个正式的语法可以被认为是定义一种语言的一套规则，或者仅仅是一种规定什么样的单词组合有意义的官方方式。在英语中，我们在学说话时从父母那里学到这些，在学校里也是如此。学习第二语言时，通常可以分为两个步骤:

1. 学习词汇
2. 学习语法

编程并没有太大的不同，编程语言之间的主要区别不是你试图实现什么，而是你如何告诉计算机去做。

那么我们如何向计算机教授语法呢？那么这些词是如何用于进化优化的呢？

这个过程大概是这样的:

1. 定义语法的巴克斯诺尔形式(BNF)表示；
2. 随机生成线性染色体来编码信息；
3. 通过将语法映射到染色体并对照适应度函数进行测试来评估染色体的适应度；
4. 使用交叉、变异、繁殖等遗传算子生成下一代染色体；
5. 重复第 3 步和第 4 步一定数量的世代。

巴克斯诺尔形式

BNF 语法是一种编写计算机可以理解的语法的方式(上下文无关)。在 Python 中，它由包含 4 个集合的元组表示:

• t:终端设备
• n:非终结集
• 生产规则集
• s:开始符号(是 N 的成员)。

这可能不能澄清事实，但是记住这一点是有好处的。让我们看一个例子:

T = {+, -, /, *, x, y}

所以终端集是你想在程序中使用的所有操作符和变量的集合。

N = {<e>, <o>, <v>}

非终结集是程序可以承担的所有生产活动的集合。在这种情况下，<e>是一个表达式，<o>是一个运算符，<v>是一个变量。

P接下来是生产规则的设定:

<e> ::= <e><o><e> | <v>
<o> ::= + | - | / | *
<v> ::= x | y

注意，P中的所有条目都是T或N的条目。管道操作符|代表or。因此，当在染色体中遇到<e>时，它要么被替换为<e><o><e>，要么被替换为<v>，对于<o>和<v>也是如此。

最后，开始符号:

S = <e>

这可能仍然很令人困惑，但是在评估染色体时，以这种形式表示语法的原因变得很明显。

染色体

描述染色体最简单的方式是作为一个非负整数数组。值并不重要，但长度很重要(稍后会详细介绍)。

举个例子，我们称我们的染色体为:

C: [4, 15, 75, 8, 41, 12]

将染色体映射到语法

现在一切都在一起了。记得我们的开始符号是:

<e>

现在，根据我们的生产规则P，当遇到一个<e>时，要么用<e><o><e>替换，要么用<v>替换。但是我们怎么知道是哪一个呢？染色体告诉我们。染色体的第一个值是4，产生式规则中有<e>的2选项。4 mod 2是0，所以我们从产生式规则<e><o><e>中取第0 -th(零索引，所以真的是1 -st 条目)选项。

所以现在我们的表达式已经扩展到了<e><o><e>，我们正在评估第一个条目<e>(因为它不是终端集的成员，所以我们一直评估每个条目，直到它被终端集的成员替换)。

我们知道在产生式规则中有<e>的2选项，所以我们取染色体中的第二个条目15，并且15 mod 2给出了1，所以我们从产生式规则<v>中取第一个(零索引，所以实际上是第二个条目)选项。

我们继续以下步骤:

• 完整表达式:<v><o><e>
  当前表达式:<v>
  语法中<v>选项个数:2
  染色体值:75``75 mod 2 = 1
  用y替换<v>
• 完整表达式:y <o><e>
  当前表达式:<o>
  语法中<o>选项个数:4
  染色体值:8
8 mod 4 = 0
  用+替换<o>
• 完整表达式:y + <e>
  当前表达式:<e>
  语法中<e>选项个数:2
  染色体值:41
41 mod 2 = 1
  用<v>替换<e>
• 完整表达式:y + <v>
  当前表达式:<v>
  语法中<v>选项个数:2
  染色体值:12
12 mod 2 = 0
  用x替换<v>

最终表情:y + x

所以我们已经用我们的语法表示把我们的染色体(数组)变成了一个函数！请随意重读几遍，因为理解表达式和运算符的来源需要一些时间。

评估适应度函数

现在我们有了一个语法上正确的函数，我们可以根据一些适应度函数来评估这个函数的性能。那些表现良好的功能有更大的机会将它们的 DNA 传递给下一代染色体。从这一点来看，GE 与其他遗传算法没有太大的不同，因为染色体被交叉和变异以产生下一代染色体，然后这些染色体被映射到语法并根据适应度函数进行测试。

这将持续用户指定的代数，最后将最有效的(适应度函数的最佳分数)算法返回给您。

通用电气公司分析

现在你(希望)对语法进化有了更多的了解，你可能会问为什么它这么好？

首先是快，真的快。通过代代相传信息，它大大改进了随机搜索。

其次，因为语法是你写的，你可以很容易地将领域知识编码到函数中。你可能会经常听到这种说法，并问“这到底是什么意思？”(我知道我有)。以下是一些例子:

• 允许访问你所有的信息，不管它是否有用。这可以通过在语法中使用一个<v>表达式来实现，并且包含你想要的所有变量。
• 根据需要使程序尽可能复杂或简单。在<e>中包含许多选项(包括重复的选项，以改变特定操作符被选中的机会)允许容易地操纵进化程序的复杂性。
• 通过设置<n> ::= 0|1|2|3|4|5|6|7|8|9然后<const> ::= <n><n><n>.<n><n><n>来定义常量的精度等级，设置精度为 3 位小数的 0 到 999 之间的常量。
• 包含表达式包含<fn>(<e>)和<fn> ::= sin | max的函数

这一切都很酷，但通用电气也有一些缺点。首先，如上所述，染色体的长度确实会导致程序无效。对于上面的例子，如果染色体中的倒数第二个条目是 42 而不是 41，那么<e>将被替换为<e><o><e>，但是没有留下染色体数据来评估这些表达式。由于该染色体无效，因此不能根据适应度函数对其进行评估。因此，虽然一条染色体可能包含接近好的解决方案，但它永远不会为人所知。

与此相反的是，程序进化时没有使用整个染色体。剩余的数据称为尾部，可以忽略。因此，通常最好将染色体构造为长尾，因为这可以减少无效解出现的机会。

然而，可能最大的缺点是，语法必须非常非常仔细地构建，为问题构建适当的语法和适应度函数需要一点艺术。除此之外，“大”语法可能需要很长时间才能进化。包括单级 for 循环、if 语句、列表函数和字典操作的语法可能需要 12 个小时以上才能完成。希望这在未来会有所改善，但这是阻止语法进化取代人类代码优化器的主要原因吗…

所以现在你有希望对语法演变有更多的了解，我们可以尝试在 PonyGE2 中实现一个简单的例子。

PonyGE2

PonyGE2 可以从这里克隆或者分叉。不需要任何设置，演化将从带有适当参数的命令行运行。

研究这个包的最好方法是通过一个例子来完成。我们将尝试开发一个函数来查找列表中的最大值。这很容易在 Python 中用max()或函数实现:

for val in list:
    if val >= current_max:
        current_max = val

第一步是创建一个适应度函数，或者我们想要最大化或最小化的函数。这里的术语有点混乱，因为我们正在进化一个函数，其结果将优化适应度函数。目前我们只是在讨论适应度函数。

这些例子的所有代码都包含在 https://github.com/Padam-0/NC-GA 的中

在PonyGE2/src/fitness中有很多例子，要运行你自己的，在那个文件夹中创建一个名为max_in_list.py的文件，结构如下:

*from* fitness.base_ff_classes.base_ff *import* base_ff
import random
import time*class* max_in_list(base_ff):
    *def* __init__(self):
        # Initialise base fitness function class.
        super().__init__() *def* evaluate(self, *ind*, ***kwargs*):
        p = *ind*.phenotype print("\n" + p) fitness = 0 for trial in range(50):
            self.test_list = generate_list()
            m = max(self.test_list) d = {'test_list': self.test_list} *try*:
                t0 = time.time()
                exec(p, d)
                t1 = time.time() guess = d['return_val'] fitness += len(p) v = abs(m - guess)
                if v <= 10**6:
                    fitness += v
                else:
                    fitness = self.default_fitness
                    break if t1 - t0 < 10:
                    fitness = self.default_fitness
                    break
                else:
                    fitness += (t1 - t0) * 1000 *except*:fitness = self.default_fitness
                break return fitness *def* generate_list():
    return [random.randint(0, round(random.random() * 90 + 10, 0)) for i in range(9)]

让我们一行一行地来。

第 1–3 行:

• 进口。base-ff是我们将要继承的类，random 用于生成随机列表，time 用于跟踪执行时间

第 5 行:

• 定义类，max_in_list，继承自base_ff。这需要与您正在处理的文件同名(在本例中是max_in_list.py)。

第 6–9 行:

• 该类的实例初始化信息。

第 10 行:

• 创建带有三个参数的evaluate函数，self、ind 和 kwargs。

到目前为止，这是每个健身功能的标准配置。下一步是独立于每一个正在进化的算法。

第 11–13 行:

• 从进化中返回当前算法作为p。我们可以将它打印到命令行(如我们在第 13 行中所做的那样)，并查看它。

第 15 行:

• 将适应性初始化为0

第 17 行:

• 将 for 循环设置为运行 50 次。如果我们不这样做，每个算法将运行在同一个“随机”列表上，并可能演化出一行，如随后返回的m=49。49可能恰好在那个列表里，而且是最大值，所以计算机认为找到了完美的算法。短、快、准！但是对于这个例子来说。在任何其他列表中，如果49不在其中(或者不是最大值)，那么算法就是错误的。因此，我们给出了适应度函数的 50 个列表，每个列表都有不同的值(实际上是最大可能值)，以避免过度适应。

第 18 行:

• 从函数generate_list()创建随机列表self.test_list
• 跳到最后几行，这个函数在0之间生成 10 个随机整数，以及 10 到 100 之间的某个数字。这是避免过度拟合的最好方法。

第 19 行:

• 找出当前使用列表的最大值。重要的是在这里完成，而不是在算法执行之后，因为在执行过程中列表可能会改变(如果语法中已经包含列表操作)，所以可能不会返回正确的结果。
• 此时你可能(正确地)会问，如果我们已经有了一个函数来寻找列表的最大值，那么做这些有什么意义呢？这是一个很好的问题，通用电气的重点是找到一种潜在的更有效的方法来做到这一点。对于简单的函数来说，这并不值得，但却是一个很好的学习练习。

第 21 行:

• 创建一个字典d，其中有一个条目test_list，用于保存当前列表。这在我们写语法的时候会更有意义。

第 23/45–47 行:

• Try/Except 语句。一般来说，拥有一个开放式 Except 语句不是一个好主意，因为你永远不会捕捉到错误，但在这种情况下，所有错误都可能发生在进化的算法中(因为它只是半随机地将单词放在一起，所以它们不可能都是完美的代码)，所以该语句是合理的。
• 如果返回一个错误，该算法的适合度被设置为默认值(NA ),循环被中断，下一个算法被测试。

第 24 行:

• 获取开始时间。

第 25 行:

• 执行算法p并在d内传递参数，以便可以访问它们。这是实际评估算法的地方。返回的值作为return_val写入字典d(我们在语法中定义了这一点)。

第 26 行:

• 获取结束时间。

第 28 行:

• 从当前算法中检索最大值猜测。

第 30 行:

• 将当前适应度分数增加算法的长度。这是为了避免臃肿和不必要的 if 语句。进化试图让事情尽可能的短。

第 32 行:

• 求列表的实际最大值和猜测值之差。最好的情况是这是0，但如果不是，这是一个算法如何接近得到正确答案的代理。

第 33–37 行:

• 有一些算法可能会产生完全超大的猜测，这将破坏进化，所以最好保留一个 if 语句来检查值是否高于某个永远不可行的大(但不太大)数。
• 同样，如果该检查失败，则将适应度设置为 NA，并继续下一个算法。否则，将猜测值和正确值之间的差加到当前的适应性分数上。

第 39–43 行:

• 一些算法(具有多层 for 循环，或for i in range(1000)语句)需要很长时间来评估。除非这是预期的，否则这表示有错误，因此建议使用 if 语句来捕获那些花费了更多时间(在本例中为 10 秒)的错误。
• 如果算法花费的时间少于此，则可以将该时间添加到适应度分数中，这样，工作更快的算法将获得较低的适应度。由于时间尺度很小，建议将这部分分数放大 100 或 1000 倍。

第 49 行:

• 还健身，就完事了！

有很多东西需要理解，有一些可能没有意义，但我们会继续写语法，希望有了完整的描述，会更清晰一些。

语法

这是困难的部分，因为你必须对你写的东西非常精确和深思熟虑。我们将逐步介绍发展 max 函数所需的基本语法，但是包括附加逻辑、列表函数和字典函数的语法扩展可在https://github.com/Padam-0/NC-GA获得。

同样，PonyGE2 在 PonyGE2/grammars 目录中提供了语法示例。在那里创建你的(它可以被称为任何东西，但是给它扩展名.pybnf)。

第一行包含程序的起点:

<fc> ::= <deff>{::}<callf>

用英语来说，这表示:

• <fc>是起始符号，包含一个选项:<deff>后面跟一个换行符({::}代表换行符，{:是缩进(tab 或 4 个空格):}关闭缩进)，后面跟<callf>。

下一行，我们定义<deff>:

<deff> ::= def fun(li):{:m = 0{::}<code>{::}return m:}

这看起来有点像 Python。我们可以看到我们有一个名为fun()的函数的函数定义，有一个参数li。

接下来我们打开一个缩进(像 Python 要求的那样)，初始化我们的变量m=0，开始一个新行，写一些<code>，开始一个新行，返回m，然后关闭缩进。

还是那句话，<deff>只有一个选项，所以这个肯定会写。接下来，<callf>:

<callf> ::= return_val = fun(test_list)

同样，这看起来更像 Python。而且只有一个选项，所以不管染色体如何，我们的算法看起来总是这样:

def fun(li):
    m = 0
    <code>
    return m
return_val = fun(test_list)

<code>部分是进化完成的地方，但是在大多数情况下，我们现在将有一个工作函数。酷嘿！天气变凉了。

还记得在我们的健身函数中，我们给了字典d一个名为test_list的键吗？嗯，因为我们调用了exec(p, d)，所以程序可以看到那个字典d的内容，所以当它在算法中看到test_list的时候，就用字典d中test_list对应的值来替换它，或者在我们的例子中，就是我们想要找到的最大值的随机列表！很酷吧，现在我们有办法将数据从我们的适应函数传递到我们的进化函数中。

同样，在我们的适应度函数中，我们设置guess = d['return_val']？事实证明，我们不仅可以从d读取，还可以向它写入，从fun(test_list)返回的值作为return_val传递给d，这样我们就可以在适应度函数中访问进化算法的返回值，这样我们就可以访问它的适应度了！

因此，我们现在可以双向发送数据，我建议重新阅读适应度函数部分，看看这是否更有意义。对我来说绝对是。

至于我们剩下的语法:

<code> ::= <stmt> | <stmt>{::}<code>

代码可以是一条语句，也可以是一条语句后在新的一行上跟随更多的代码。这样我们就可以给自己多线功能了！这种语法的递归性质是非常有用的，将会出现很多，使我们能够增加进化的复杂性。

<stmt> ::= <var> = <expr> | <for> | <if>

然后，语句要么将表达式赋给变量，要么赋给 for 循环，要么赋给 if 语句。让我们以相反的顺序来看这些:

<if> ::= if <cond>:<if-opt>
<if-opt> ::= {:<code>:} | {:<code>:}else:{:<code>:} | {:<code>:}elif <cond>:{:<if-opt>:}<cond> ::= <expr> <c-op> <expr>
<c-op> ::= ">=" | "<=" | ">" | "<"

所以 if 语句接受一个条件，然后是一个选项。条件采用表达式、运算符和表达式的形式。

为了使这一点更清楚，现在有必要讨论一下<expr>标签:

<expr> ::= <number> | <var> | <expr> <op> <expr>

所以一个表达式要么是一个数字，要么是一个变量，要么是一个操作(后面会详细介绍)。

回到我们的 if 语句，我们的条件采用形式if <cond>:,它扩展为:

if <expr> <c-op> <expr>:

然后可以扩展为:

if <var> >= <number>:

我们有一个我们认可的 if 语句。<if-opt>稍微复杂一点，因为它们必须包括elif和else功能，但是它们都归结为一个 if 语句、一个条件，然后一些<code>被写入其中。

接下来，对于循环:

<for> ::= for i in <list-var>:{:<fl-code>:}
<fl-stmt> ::= <var> = <expr> | <fl-if>
<fl-code> ::= <fl-stmt> | <fl-stmt>{::}<fl-code>

因此，for 循环采用形式for i in list:,然后为某个语句创建一个缩进(在新的一行上)。语句本身可能只是我们上面的<stmt>，但是这将允许嵌套 for 循环，这会成倍地增加计算时间。除非您认为嵌套的 for 循环是绝对必要的，否则请避免使用它们(一般来说，这是一个非常好的编程规则)。

所以我们的<fl-stmt>只包括变量赋值，和一个特殊的 if 语句，<fl-if>:

<fl-if> ::= if <cond>:<fl-if-opt>
<fl-if-opt> ::= {:<fl-code>:} | {:<fl-code>:}else:{:<fl-code>:} | {:<fl-code>:}elif <cond>:{:<fl-if-opt>:}

这实际上与上面的 if 语句语法相同，但是它不调用<code>，而是调用<fl-code>，基本上是为了避免 for 循环内部的 if 语句内部出现 for 循环。有点麻烦，但比允许嵌套 for 循环要好！

最后，变量赋值。我们已经看过<expr>，所以让我们看看组成它的零件:

<var> ::= m | i
<list-var> ::= li<number> ::= <num><num><num> | <num><num> | <num>
<num> ::= 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9<op> ::= + | - | * | / | // | %

这应该很简单。我们的变量是开始时初始化的m和运行 for 循环时初始化的i。数字可以从 0 到 999，我们有所有的 python 运算符(除了指数运算符，它会导致整数溢出和一大堆执行非常慢的算法。像嵌套的 for 循环一样，只在需要时才包含它)。

信不信由你，这就是我们的语法！让我们看看如何从一个<code>标签到我们的最终函数。像这样的染色体:

[0, 1, 0, 1, 0, 1, 1, 2, 1, 0, 0, 0, 0, 1, 1]

绘制出来后，我们可以得到:

def fun(li):
    for i in li:
        if i > m:
            m = i
return_val = fun(test_list)

太棒了。现在我们有了适应度函数和语法，我们如何用 PonyGE2 运行它呢？

好吧，如果 fitness 函数在 PonyGE2/src/fitness 中(并且与它保存的类具有相同的文件名，在本例中是 max_in_list)，并且语法在 PonyGE2/grammars 中，我们可以从命令行运行:

python ponyge.py --fitness_function max_in_list --grammar_file max_in_list.pybnf

这将运行您的语法和健身功能的演变，与其他选项的默认参数。有很多，值得一读参考文献，但两个真正重要的是人口规模和世代数。

群体大小是在开始时产生的染色体的数量，并且在每一代进化。这个数字越大，搜索范围越广，但进化时间越慢。

代数是染色体进化的次数。同样，它扩大了搜索范围，但以进化速度为代价。

您可以通过以下方式设置这些参数:

python ponyge.py --fitness_function max_in_list --grammar_file MIL_e2.pybnf --population_size 500 --generations 100

现在，对于简单的语法，这些参数的值在 100-500 的范围内就足够了，但是随着语法变得更加复杂，它们需要更多的群体和世代来确保进化出足够的算法。一个好的测试是给你的语法增加复杂性，用和上面一样的适应度函数来运行它。如果你得到一个在列表中找到最大值的函数，你可以确信你的参数是合适的，然后把它转换到你要解决的问题的另一个适应度函数上。

不幸的是，运行时间变得很快。对我来说，2000 代和 750 代的人口需要大约 6 个半小时，这不足以用上面的语法进化出一个合适的算法(但也包括列表和字典操作)。

这就是通用电气的问题所在。承诺是，如果你能写一个适应度函数，和一个包含大量 Python 操作的健壮语法，那么 ge 将能进化出一种更有效的做事方式。不幸的是，目前这种处理能力和时间在公共机器上是不可行的。

不管怎样，这是一个值得关注的领域，就像我们要教机器编码一样，很有可能语法进化与此有关。

其他示例

UCD 自然计算研究和应用集团的团队是这项研究的前沿，他们发表了许多展示通用电气力量的伟大论文(尽管是在超级计算机上)。其中包括:

• 寻找正则表达式的改进性能；
• 用语法进化优化复杂的桥塔结构；和
• 结合结构分析和多目标标准的进化建筑设计。

如果你喜欢这篇文章，请点击❤按钮让你的追随者知道，或者让我知道你的想法。

Python 中的功耗分析简介

原文：https://towardsdatascience.com/introduction-to-power-analysis-in-python-e7b748dfa26?source=collection_archive---------2-----------------------

Source: unsplash.com

了解概念的重要性，如显著性水平、效应大小、统计功效和样本大小

如今，许多公司——网飞(T1)、亚马逊(Amazon)、T2()、优步(T3)，但也有一些规模较小的——不断进行实验(A/B 测试)，以测试新功能并实现那些用户认为最好的功能，最终带来收入增长。数据科学家的角色是帮助评估这些实验，换句话说，验证这些测试的结果是否可靠，是否可以/应该在决策过程中使用。

在本文中，我介绍了功耗分析。简而言之，权力是用来报告从实验结果中得出的结论的可信度。它还可以用于估计实验所需的样本量，，即，在给定的置信度下，我们应该能够检测到一个效应的样本量。实际上，人们可以理解许多事情，例如，在一个组内更频繁的转换，以及在网上商店中经历某一注册流程的顾客的更高的平均花费，等等。

首先，我介绍了一些理论，然后用 Python 实现了一个功耗分析的例子。你可以在文末找到我回购的链接。

简介

为了理解功效分析，我认为理解三个相关概念是很重要的:显著性水平、I/II 类错误和效应大小。

在假设检验中， 显著性水平 (通常表示为希腊字母 alpha)是拒绝零假设(H0)的概率，当它实际上为真时。与显著性水平密切相关的指标是p-值，它是在 H0 为真的情况下，获得至少是极端结果(离零假设更远的结果)的概率。这在实践中意味着什么？在从总体中随机抽取样本的情况下，观察到的效应可能仅仅是由于抽样误差而产生的。

当相关的p-值小于所选的α时，实验结果(或例如线性回归系数)在统计上是显著的。显著性水平应在设置研究之前指定，并取决于研究领域/业务需求。

第二个值得一提的概念是在统计测试假设时我们可能犯的种错误。当我们拒绝一个真正的 H0 时，我们谈论的是一个I 型错误(假阳性)。这是与显著性水平相关的误差(见上文)。另一种情况发生在我们未能拒绝假 H0 时，这被认为是一个II 型错误(假阴性)。你可以从混乱矩阵中回忆起这些概念！

Types of error in a nutshell 😃

最后要考虑的是 效应大小 ，这是一种现象在人群中出现的量化幅度。根据具体情况，可以使用不同的指标来计算效果大小，例如:

• 皮尔逊相关
• 回归系数
• 两组之间均值的差异，例如 Cohen 的 d

我在另一篇文章的中描述了效果大小的各种度量。

统计能力

现在我们已经修改了与功效分析相关的关键概念，我们终于可以谈论统计功效了。假设检验的统计功效简单来说就是当备选项事实上为真时，给定检验正确拒绝零假设(这意味着与接受 H1 相同)的概率。

更高的实验统计能力意味着犯第二类错误的概率更低。这也意味着当存在要检测的效应时，检测到效应的概率更高(真阳性)。这可以用下面的公式来说明:

Power = Pr(拒绝 H0 | H1 为真)= 1 - Pr(拒绝 H0 失败| H0 为假)

在实践中，功率太小的实验结果会导致错误的结论，进而影响决策过程。这就是为什么只能考虑具有可接受功率水平的结果。设计 80%功效水平的实验是很常见的，这相当于犯第二类错误的概率为 20%。

功率分析

功耗分析基于以下构建模块:

• 显著性水平
• 效果大小
• 力量
• 样本量

我以前没有讨论过样本大小，因为这是不言自明的。唯一值得一提的是，一些测试联合考虑两组的样本量，而对于另一些测试，样本量必须单独指定(在它们不相等的情况下)。

这四个指标是相互关联的。例如:降低显著性水平会导致功效降低，而较大的样本会使效果更容易检测。

功率分析的思想可以归结为以下几点:通过四个度量中的三个，我们估计缺少的一个。这在两个方面很方便:

• 当我们设计一个实验时，我们可以假设什么样的显著性水平、功效和效果大小是我们可以接受的，并且——因此——估计我们需要收集多大的样本来进行这样的实验以产生有效的结果。
• 当我们验证一个实验时，我们可以看到，在给定所用样本大小、效应大小和显著性水平的情况下，从商业角度来看，犯第二类错误的概率是否可以接受。

除了为给定指标计算一个值之外，我们还可以通过多次执行功效分析(针对组件的不同值)并在图上显示结果来执行一种敏感性分析。这样我们可以看到——例如——必要的样本量是如何随着显著性水平的增加或减少而变化的。这自然可以扩展到 3 个度量的 3D 平面。

Python 中的例子

在本例中，我对独立双样本 t 检验(样本大小和方差相等)的情况进行了功效分析。库statsmodels包含对一些最常用的统计测试进行功耗分析的函数。

让我们从一个简单的例子开始，假设我们想知道我们的实验需要收集多大的样本，如果我们接受 80%的功效水平，5%的显著性水平和预期的效应大小为 0.8。

首先，我们导入所需的库。

然后，我们需要运行以下命令，并达到所需的样本大小 25。

做了这些之后，是时候更进一步了。我们希望看到当我们修改其余的构建模块时，功率如何变化。为此，我们绘制了功率与其他参数的关系。我从检查样本大小如何影响功效开始分析(同时将显著性水平和效应大小保持在特定水平)。我已经选择了[0.2，0.5，0.8]作为考虑的效果大小值，因为它们对应于小/中/大的阈值，如在科恩的 d 的情况中所定义的。

从图中，我们可以推断，样本/效果大小的增加会导致功效的增加。换句话说，样本越大，功率越高，保持其他参数不变。下面，我还在 x 轴上展示了两个剩余构建模块的图，结果不言自明。

Power vs. effect size

Power vs. significance level

最后，我想把分析扩展到三个方面。为此，我将显著性水平固定在 5%(这是实践中经常使用的)，并创建一个可能的样本和效应大小组合的网格。然后我需要获得每个组合的功率值。为此，我使用 NumPy 的meshgrid和vectorize。

为了创建 3d 情节，我选择了plotly，因为它很容易快速获得好的，互动的情节，然后可以嵌入到这个帖子中。在代码中，我使用plotly的离线模式，不需要注册。如果你想创建可以共享/嵌入的图形，你必须在这里创建一个免费账户来获得一个 API 密匙。

结论

总的来说，功耗分析现在主要用于 A/B 测试，在计划实验/研究或评估结果时都可以使用。由于许多公司使用频率主义方法进行假设检验，知道如何进行权力分析以及如何展示其含义无疑是件好事。

一如既往，我们欢迎任何建设性的反馈。可以在 Twitter 或者评论里联系我。文章的代码可以在这里找到。

喜欢这篇文章吗？成为一个媒介成员，通过无限制的阅读继续学习。如果你使用这个链接成为会员，你将支持我，不需要你额外付费。提前感谢，再见！

SQL 中的过程和游标介绍

原文：https://towardsdatascience.com/introduction-to-procedures-and-cursors-in-sql-f9d9b9ea1fe7?source=collection_archive---------6-----------------------

Photo by Caspar Rubin on Unsplash

学习如何为一个 RDBMS 编写程序和游标。

如果你想从数据科学的角度了解更多关于 SQL 的知识，你可以参加 DataCamp 的免费课程“数据科学的 SQL 介绍”。

SQL 是任何现代软件工程师的必备技能。因为大多数软件依赖于某种类型的数据，并且与 RDBMS(关系数据库管理系统)集成得很好。无论是 web 应用程序、API 还是内部应用程序，RDBMS 都在那里。SQL 是查询 RDBMS 的语言。

作为一名数据科学家，了解 SQL 及其相关技术是非常初级的。为了能够查询 RDBMS 并获得关于您正在处理的数据的特定问题的答案，SQL 是最低要求。

在他与 DataCamp 的最新视频中，大卫·罗宾逊(首席数据科学家@ DataCamp) 向我们展示了他如何在一个数据科学问题中使用 SQL。请检查一下，他的工作流程非常有趣。

在本教程中，您将学习编写过程和游标；SQL 的另一个重要方面。您是否曾经希望 RDBMS 在执行特定操作时自动执行某些操作？例如，假设您已经在名为Employees的表中创建了一个新的雇员记录，并且您希望它反映在其他相关的表中，如Departments。好吧，你将选择正确的教程。

在本教程中，您将学习:

• RDBMS 中的过程是什么？
• 如何编写一个过程呢？
• 不同类型的程序
• RDBMS 中的游标是什么？
• 如何编写不同类型的游标？
• 不同类型的光标

听起来很刺激？让我们开始吧。

RDBMS 中的过程是什么？

在继续学习过程和游标之前，您需要了解一些关于PL/SQL的知识，它是一种块结构语言，使像您这样的开发人员能够将 SQL 的强大功能与过程语句结合起来。但是你不会以传统的方式学习，你会随着你的前进和需要而学习。

所以如果你有一个 SQL 查询，你想多次执行它。程序是解决方法之一。通常在这种情况下调用过程，因为它们保持存储状态，并在特定操作或一系列操作时被触发。程序也被称为Procs。

现在您将看到如何编写一个过程。

写作程序:

编写过程的一般语法如下:

**CREATE** **PROCEDURE** procedure_name
**AS**
sql_statement
**GO**;

请注意，这些语法适用于几乎所有 RDBMS，无论是 Oracle、PostgreSQL 还是 MySQL。

创建过程后，您必须执行它。下面是它的语法。

**EXEC** procedure_name;

现在让我们写一个简单的程序。考虑下面来自 RDBMS 的快照，它包含一个名为Customers的表。

Source: W3Schools

您将编写一个名为SelectAllCustomers的程序，它将从Customers中选择所有的客户。

**CREATE** **PROCEDURE** SelectAllCustomers
**AS**
**SELECT** * **FROM** Customers
**GO**;

执行SelectAllCustomers的人:

**EXEC** SelectAllCustomers;

过程也可以是独立的语句块，这使得它们独立于任何表，不像前面的表。下面的示例创建了一个显示字符串“Hello World！”的简单过程作为执行时的输出。

**CREATE** **PROCEDURE** welcome
**AS**
**BEGIN**
dbms_output.put_line('Hello World!');
**END**;

有两种执行独立过程的方法。

• 使用EXEC关键字
• 从 PL/SQL 块调用过程名

可以使用关键字EXEC调用上述名为“welcome”的过程，如下所示:

**EXEC** welcome;

现在您将看到下一个方法，即从另一个 PL/SQL 块调用过程。

**BEGIN**
welcome;
**END**;

程序也可以被替换。你只需要在创建程序时添加REPLACE关键字。这将替换已经存在的过程(如果),否则将创建一个新的过程。

**CREATE** **OR** **REPLACE** **PROCEDURE** welcome
**AS**
**BEGIN**
dbms_output.put_line('Hello World!');
**END**;

删除存储过程没什么大不了的:

**DROP** **PROCEDURE** **procedure**-name;

根据参数的不同，程序也会有所不同。可以有单参数过程，也可以有多参数过程。现在你将研究这些变体。

为此，您将使用同一个表Customers。为了方便起见，下面的部分再次给出了快照。

Source: W3Schools

您将编写一个存储过程，从表中选择特定城市的客户:

**CREATE** **PROCEDURE** SelectAllCustomers @City nvarchar(30)
**AS**
**SELECT** * **FROM** Customers **WHERE** City = @City
**GO**;

让我们在这里剖析一下共同的原则:

• 您编写了第一个@City，并将其类型和大小定义为程序执行时将给出的参数之一。
• 第二个@City被分配给条件变量City，它只是Customers表中的一列。

该过程执行如下:

**EXEC** SelectAllCustomers City = "London";

现在让我们看看另一个变体。

多参数的编写程序与前面的完全相同。你只需要添加它们。

**CREATE** **PROCEDURE** SelectAllCustomers @City nvarchar(30), @PostalCode nvarchar(10)
**AS**
**SELECT** * **FROM** Customers **WHERE** City = @City **AND** PostalCode = @PostalCode
**GO**;

按照以下方式执行程序:

**EXEC** SelectAllCustomers City = "London", PostalCode = "WA1 1DP";

上面的代码可读性不是很好吗？当代码可读时，做起来更有趣。手续到此为止。现在您将学习光标。

RDBMS 中的游标是什么？

Oracle 等数据库创建了一个内存区域，称为上下文区域，用于处理 SQL 语句，其中包含处理该语句所需的所有信息，例如处理的行数。

光标是指向该上下文区域的指针。PL/SQL 通过一个光标控制上下文区域。游标是执行 SQL 语句时在系统内存中创建的临时工作区。游标包含有关 select 语句及其所访问的数据行的信息。因此，游标用于加快大型数据库中查询的处理时间。您可能需要使用数据库游标的原因是，您需要对单独的行执行操作。

光标有两种类型:

• 隐式光标
• 显式光标

现在您将看到如何编写不同类型的游标。

写入光标:

您将从理解什么是隐式游标开始这一部分。

当没有为 SQL 语句定义显式游标时，每当执行该语句时，Oracle 都会自动创建隐式游标。程序员不能控制隐式光标和其中的信息。每当发出 DML(数据操作语言)语句(INSERT、UPDATE 和 DELETE)时，都会有一个隐式游标与该语句相关联。对于插入操作，光标保存需要插入的数据。对于更新和删除操作，光标标识将受影响的行。

您可以将最近的隐式游标称为 SQL 游标，它始终具有如下属性:

• 找到%个，
• %ISOPEN，
• %未找到，
• %ROWCOUNT

下图简要描述了这些属性:

Source: TutorialsPoint

让我们考虑一个数据库的快照，它包含一个名为Employee的表:

现在，您将编写一个游标，将年龄小于 30 岁的人的薪水增加 1000 英镑。

**DECLARE**
total_rows number(2);**BEGIN**
**UPDATE** Employee
**SET** salary = salary + 1000
**where** age < 30;
**IF** **sql**%notfound **THEN**
    dbms_output.put_line('No employees found for under 30 age');
**ELSIF** **sql**%**found** **THEN**
    total_rows := **sql**%rowcount;
    dbms_output.put_line( total_rows || ' employees updated ');
**END** **IF**;
**END**;

现在让我们回顾一下你所写的一切:

• 您定义了一个名为total_rows的变量，用于存储将受光标操作影响的雇员数量。
• 您用BEGIN开始了 cursors 块，并编写了一个简单的 SQL 查询来更新那些年龄小于 30 岁的人的工资。
• 如果数据库中没有雇员年龄小于 30 岁的条目，您可以处理输出。您为此使用了%notfound属性。注意，这里的隐式光标sql存储了所有相关信息。
• 最后，您使用%rowcount属性打印出受游标影响的记录数。

太好了！你做得很好！

当在 SQL 提示符下执行上述代码时，它会产生以下结果:

更新了 2 名员工(假设有 2 条记录，其中年龄< 30)

现在，您将学习显式游标。

显式光标对上下文区域进行更多定义的控制。它是在返回多行的 SELECT 语句上创建的。

创建显式游标的语法是

**CURSOR** **cursor_name** **IS** select_statement;

如果您正在使用显式游标，您需要遵循如下一系列步骤:

• 在内存中声明用于初始化的光标
• 打开光标分配内存区域
• 获取用于检索数据的光标
• 关闭光标以释放内存

下图表示典型显式游标的生命周期:

现在，您将进一步了解这些步骤。

声明光标:

您声明了一个游标和一个SELECT语句。例如:

**CURSOR** **C** **IS** **SELECT** id, name, address **FROM** Employee **where** age > 30;

打开光标:

当您打开游标时，CPU 会为游标分配内存，并准备好获取 SQL 语句返回的行。例如，我们将打开上面定义的光标，如下所示:

**OPEN** **C**;

获取光标:

获取游标涉及到从游标所需的 SQL 中的关联表中一次访问一行。

**FETCH** **C** **INTO** C_id, C_name, C_address;

关闭光标:

关闭游标意味着释放分配的内存。您将关闭上面打开的光标，如下所示:

**CLOSE** **C**;

你现在要以一种有意义的方式把所有这些部分组合在一起。

组装这些零件:

**DECLARE**
C_id Employee.ID%**type**;
C_name Employee.NAME%**type**;
C_address Employee.ADDRESS%**type**;
**CURSOR** **C** **is**
**SELECT** id, name, address **FROM** Employee **where** age > 30;
**BEGIN**
**OPEN** **C**;
LOOP
**FETCH** **C** **INTO** C_id, C_name, C_address; 
dbms_output.put_line(ID || ' ' || NAME || ' ' || ADDRESS); 
EXIT **WHEN** **C**%notfound;
**END** LOOP;
**CLOSE** **C**;
**END**;

您还学习了声明游标变量 C_id、C_name 和 C_address。C_id Employee.ID%type; -这确保了创建 C_id 的数据类型与在Employee表中 id 的数据类型相同。

通过使用LOOP,您可以通过光标循环读取记录并显示出来。如果光标没有找到记录，您也可以处理这种情况。

当在 SQL 提示符下执行代码时，它会产生

结论:

恭喜你。你已经坚持到最后了。您讨论了数据库世界中最流行的两个主题——过程和游标。这些在处理大量事务的应用程序中很常见。是的，你猜对了！银行自古以来就在使用这些。您学习了如何编写一个过程，它有哪些不同的类型以及为什么会这样。您还学习了游标及其几种变体，以及如何编写它们。

太神奇了！

以下是撰写本教程的一些参考资料:

• Oracle PL/SQL 编程
• 光标上的 TutorialsPoint 博客
• SQL 存储过程— W3Schools

Q 学习简介

原文：https://towardsdatascience.com/introduction-to-q-learning-88d1c4f2b49c?source=collection_archive---------1-----------------------

想象自己在迷宫中寻宝。游戏如下:

你从一个给定的位置开始，起始状态。从任何状态你都可以向左、向右、向上或向下走，或者停留在同一个地方，只要你不穿过迷宫的前提。每个动作将带你到网格的一个单元(不同的状态)。现在，在其中一个州(目标州)有一个宝箱。此外，迷宫在某些位置/状态有一个蛇坑。您的目标是沿着一条没有蛇的路径从起始状态行进到目标状态。

Grid outline of the maze

当您将一个代理放入网格(我们称之为我们的环境)时，它将首先探索。它不知道什么是蛇，也不知道宝藏在哪里。所以，为了给它蛇和宝箱的概念，我们会在它完成每个动作后给它一些奖励。它每踏上一个蛇穴，我们将给予它-10 的奖励。对于宝藏，我们将给予+10 的奖励。现在我们希望我们的代理尽快完成任务(走最短的路线)。为此，我们将给予其他州-1 的奖励。然后我们会告诉它最大化分数。现在，随着代理探索，它了解到蛇对它是有害的，宝藏对它是有益的，它必须尽快得到宝藏。图中的’-'路径表示奖励最大的最短路径。

Q-Learning 试图了解处于特定状态的价值，并在那里采取特定的行动。

我们要做的是开发一个表格。其中行是状态，列是它可以采取的动作。因此，我们有一个 16x5 (80 个可能的状态-动作)对，其中每个状态是迷宫网格的一个单元。

我们首先将该表初始化为统一的(全零)，然后当我们观察到我们从各种行为中获得的奖励时，我们相应地更新该表。我们将使用 贝尔曼方程来更新表格。

eqn.1

“s”代表当前状态。‘a’表示代理从当前状态采取的动作。s”表示该操作产生的状态。“r”是你采取行动得到的奖励，“γ”是折扣系数。所以，采取行动 a 的状态的 Q 值是即时奖励和未来奖励的折扣之和(结果状态的值)。折扣系数‘γ’决定了你对未来奖励的重视程度。比方说，你去了一个离目标州更远的州，但是从那个州开始，遇到有蛇的州的机会就更少了，所以，这里未来的奖励更多，尽管即时奖励更少。

我们将每次迭代(代理所做的尝试)称为一个情节。对于每一集，代理将试图达到目标状态，并且对于每个转换，它将继续更新 Q 表的值。

让我们看看如何计算 Q 表:

为此，为了方便起见，我们将采用一个较小的迷宫网格。

最初的 Q 表看起来像这样(状态沿着行，动作沿着列) :

Q Matrix

上、下、左、右

奖励表应该是这样的:

R Matrix

这里，E 表示零，即，不可能有这样的转变。

算法:

1. 用全零初始化 Q 矩阵。设置“γ”的值。填写奖励矩阵。
2. 每集。选择一个随机的起始状态(这里我们将起始状态限制为 state-1)。
3. 从当前状态的所有可能操作中选择一个。
4. 作为动作(a)的结果，移动到下一个状态(S’)。
5. 对于状态中所有可能的动作，选择具有最高 Q 值的一个。
6. 使用公式 1 更新 Q 表。
7. 将下一个状态设置为当前状态。
8. 如果达到目标状态，则结束。

例如:假设我们从状态 1 开始。我们可以选择 D 或 r。比如说，我们选择了 D。然后我们将到达 3(蛇坑)。因此，我们可以选择 U 或 R，取γ = 0.8，我们有:

Q(1，D) = R(1，D) + γ*[max(Q(3，U) & Q(3，R))]

Q(1，D) = -10 + 0.8*0 = -10

这里，max(Q(3，U) & Q(3，R)) = 0，因为 Q 矩阵尚未更新。-10 代表踩蛇。因此，新的 Q 表看起来像:

现在，3 是起始状态。从 3 开始，假设我们去 r。所以，我们去 4。从 4 开始，我们可以去 U 或 L。

Q(3，R) = R(3，R) + 0.8*[max(Q(4，U) & Q(4，L))]

Q(3，R) = 10 + 0.8*0 = 10

所以，现在我们已经达到了目标状态 4。所以，我们终止和更多的传递，让我们的代理理解每个状态和动作的值。继续传递，直到值保持不变。这意味着您的代理已经尝试了所有可能的状态-动作对。

python 中的实现:

最后一个 q_matrix 的输出:

在下一篇文章中，我将介绍使用神经网络进行 Q 学习，以及使用表格方法的缺点。此外，我们将致力于开放人工智能健身房的游戏测试。在那之前，再见。

量化金融导论:资产回报的程式化事实

原文：https://towardsdatascience.com/introduction-to-quantitative-finance-part-i-stylised-facts-of-asset-returns-5190581e40ea?source=collection_archive---------6-----------------------

作为一名具有定量金融背景的数据科学家，我一直对探索两个领域结合的可能性感兴趣。我认为这是一个令人着迷的探索领域，这就是为什么我想开始一系列的文章来描述量化金融的基础知识。

在本系列的最后，我打算提出一个简单的分配策略，并希望表明，通过使用数据科学/量化金融知识，有可能超越基本的基准策略。当然，没有人在谈论建立一个完美的模型，准确预测未来的股票价格，并使我们——潜在的投资者——赚取数百万美元。

有一个理论解释了为什么这是不可能的，即有效市场假说(EMH)。它指出，资产价格充分反映了所有可用的信息。这意味着持续跑赢市场是不可能的，因为市场价格只会对新信息做出反应。你可以在这里了解更多关于的信息。但是我们仍然可以尝试一下，看看短期内是否有可能赚到一些钱，至少在理论上是这样的:)

在直接进入机器学习和建立资产配置策略之前，我认为花一些时间在基础知识上并理解我们试图建模的过程是至关重要的。在本文中，我研究了资产回报的程式化事实，并展示了如何使用 Python 来验证它们的存在。一些统计学的基础知识会有所帮助，但是当我试图直观地解释正在发生的事情时，不要气馁。

回报以及我们为什么与他们合作

我首先使用 Python 的quandl库下载历史股票价格。这非常简单，你只需要创建一个免费账户来获得 API 密匙。

在本文中，我使用调整后的收盘价，因为它们说明了可能的公司行为，如股票分割等。

我选择微软(股票代码:MSFT)作为例子，以数据帧的形式下载时间序列。然后，我将价格转换成对数回报，以便进一步分析:

其中 P_t 表示资产在时间 t 的价格，log 函数代表自然对数(有时也称为 ln )。

有一篇很棒的文章描述了简单返回和日志返回之间的区别，因此您可以在那里找到使用这两种方法的优缺点。

但是当我们已经有了价格，为什么还要处理退货呢？原因是价格通常是非平稳的，即均值和方差(数学矩)等统计数据随时间变化。这也可能意味着在我们的价格序列中观察一些趋势或季节性。正如你所想象的，通过处理回报率，我们使时间序列稳定，这是统计建模中的一个期望属性。现在，就让它这样吧。

Decomposition of a time series into trend/seasonal components.

下面我展示了微软价格和回报随时间的演变。

在图中可以直接观察到的一个重要事实是“波动聚集”的存在——大回报周期与小回报周期交替出现，表明波动不是恒定的。我稍后会回来。

一般来说，程式化的事实是出现在许多经验资产回报(跨时间和市场)中的统计属性。了解它们是很重要的，因为在构建代表资产价格动态的模型时，模型必须能够捕捉/复制这些属性。

事实 1:收益分布不正常

Standard Normal (Gaussian) distribution

据观察，回报表现为:

• 负偏度(三阶矩)——大的负回报比大的正回报更常见。目测:左尾较长；分布的质量集中在分布图的右侧。
• 超额峰度(四阶矩)——大(和小)回报比预期的更经常发生。目测:厚尾、尖峰分布。

下面，我展示了一个直方图，直观地显示了微软的对数回报的分布，以及一条代表正态概率密度曲线的线(平均值和标准偏差等于样本平均值)。我们看到，回报确实呈现出更高的峰值，而且更多的质量位于尾部(比正常情况下预期的要多)。

为了进一步检验，我还看了 Q-Q 图。红线代表标准的正态分布。在收益服从高斯分布的情况下，这两条线会对齐。然而，我们看到有差异，主要是在尾部。这进一步证实了上述发现。

关于如何创建和解释 Q-Q 图的更多细节，你可以阅读本文。

最后，我查看了所考虑的回报的描述性统计数据(它们以日价值表示，实际上，通常以年价值表示)。 Jarque-Bera 正态性检验证实了我们的怀疑，p 值小到足以拒绝声称数据遵循高斯分布的零假设。

Descriptive statistics of MSFT’s log returns

事实 2:回报没有(或几乎没有)显著的自相关性

自相关测量给定时间序列与连续时间间隔内同一序列的滞后版本之间的相似程度。这类似于两个时间序列之间的相关性:第一个是原始形式，另一个滞后了 n 个周期。

例如:当某种资产的回报表现出历史上的正自相关性，并且在过去几天里价格正在上涨，人们可能会合理地预期进一步的正运动(当然，预测股票价格并不那么简单，否则成为百万富翁将是一件相当容易的事情……EMH 再次出击)。

事实 3:平方收益和绝对收益的自相关性小且缓慢下降

在对回报建模时，考虑波动性在决策(买/卖)过程中可能是至关重要的。波动性通常被理解为回报的标准差(方差的平方根)。

现在，让我们考虑误差，而不是回报，即实际值——由某个模型预测/解释的值。方差基本上是误差平方的平均值，而绝对偏差是绝对误差的平均值。通过绘制随时间变化的平方/绝对误差，我们可以看到方差(或绝对偏差，也是波动性的一种度量)是否随时间保持不变。在资产回报的情况下，情况并非如此，我们可以观察到高/低波动的时期。这被称为“波动聚类”，可以在对数收益的时间序列图上观察到。

另一方面，长期短期(每日)平均回报预计为零(EMH)。这就是为什么通过查看平方和绝对回报，我们可以有效地测量与预期均值的偏差，而不用查看误差的方向——平方和绝对函数都会抵消误差的方向。

下面我展示了 MSFT 对数收益的自相关图，以及平方值和绝对值(检验事实 2 和 3)。蓝色区域表示 95%的置信区间，超出该区间的点具有统计学意义。我们看到，对于对数回报，只有几个重要的点(这与事实 2 一致)。至于事实 3，我们看到相关性是显著的，平方收益比绝对收益更容易观察到相关性的下降。总之，这让我们相信，我们可以尝试利用自相关结构来进行波动建模。

ACF Plots of returns, squared and absolute returns

结论

在本文中，我简要介绍了资产回报的程式化事实。当然，这仅仅是触及了话题的表面。还有更多的东西需要学习，不同的资源也增加了更多的程式化事实。如果你感兴趣，我在参考资料中添加了一些额外的资源。

需要注意的一件重要事情是，不能保证所有的事实在每个收益序列中都清晰可见。这很容易受到分析中考虑的时间范围的影响。在这种情况下，我从 2000 年开始收集数据。但是没有严格的规定。使用提供的代码，您可以轻松地调查不同的时间范围和资产。

介绍完之后，在下一篇文章中，我将介绍一些预测资产回报的常用方法(也许还有它们的波动性)。在最后一篇文章中，我将试图说明，通过(或多或少准确的)预测，我们可以尝试建立资产配置策略，并调查我们是否在历史上赚了一些钱！

一如既往，我们欢迎任何建设性的反馈。你可以在推特上或评论中联系我。用于调查程式化事实的代码可以在我的 GitHub 上找到。

我最近出版了一本关于使用 Python 解决金融领域实际任务的书。如果你感兴趣，我在贴了一篇文章介绍这本书的内容。你可以在亚马逊或者 Packt 的网站上买到这本书。

参考资料:

[1]https://orfe.princeton.edu/~jqfan/fan/FinEcon/chap1.pdf

[2]资产回报的经验属性:程式化的事实和统计问题—http://finance . martinse well . com/typed-facts/dependency/cont 2001 . pdf

推荐系统介绍。第 2 部分(神经网络方法)

原文：https://towardsdatascience.com/introduction-to-recommender-system-part-2-adoption-of-neural-network-831972c4cbf7?source=collection_archive---------1-----------------------

1.介绍

在我本系列的上一篇博文: 推荐系统介绍。第一部分(协同过滤，奇异值分解) ，我谈到了协同过滤(CF)和奇异值分解(SVD)如何用于构建推荐系统。随着神经网络的兴起，你可能会好奇我们如何利用这种技术来实现推荐系统。这篇博文将介绍由 PyTorch 支持的推荐系统框架 Spotlight，以及我借用单词嵌入的思想创建的 Item2vec。最后，我将比较我讨论过的模型的性能。

2。聚光灯

Spotlight 是一个实现良好的 python 框架，用于构建推荐系统。它包含两大类模型，因子分解模型和序列模型。前者利用了 SVD 背后的思想，将效用矩阵(记录用户和项目之间交互的矩阵)分解成用户和项目矩阵的两个潜在表示，并将它们馈入网络。后者是用长短期记忆(LSTM)和一维卷积神经网络(CNN)等时间序列模型建立的。由于 Spotlight 的后端是 PyTorch，所以在使用之前，请确保您已经安装了正确版本的 PyTorch。

相互作用

效用矩阵被称为聚光灯下的交互。为了创建隐式交互，我们为所有用户-项目交互对指定用户和项目的 id。额外的评级信息将这种交互转变为明确的交互。

因子分解模型

因式分解模型包含隐式或显式的交互。为了便于说明，我们将使用隐式方法。

Implicit Feedback (https://github.com/maciejkula/netrex/)

它的思想与 SVD 非常相似，在 SVD 中，用户和项目被映射到一个潜在的空间中，这样它们就可以直接进行比较。本质上，我们使用两个嵌入层来分别表示用户和项目。目标是我们传入的交互(效用矩阵)。为了计算用户-项目对的分数，我们获取该用户和项目的潜在表示的点积，并将其传递给 sigmoid 激活函数。通过计算所有用户-项目对相对于真实交互的损失(稍后将详细介绍),我们可以反向传播并优化嵌入层。网络结构如下图所示。

Neural Network Structure (https://github.com/maciejkula/netrex/)

我们只需要几行代码来用 Spotlight 训练这样的模型，它看起来非常类似于 scikit-learn toolkit:

顺序模型

顺序模型将推荐问题视为顺序预测问题。给定过去的交互，我们想知道在下一个时间步骤中用户最有可能喜欢哪个项目。例如，假设用户 A 与 id 顺序为[2，4，17，3，5]的项目进行交互。那么我们将有下面的扩展窗口预测。

[2] -> 4
[2, 4] -> 17
[2, 4, 17] -> 3
[2, 4, 17, 3] -> 5

左边的数组存储过去的交互，而右边的整数表示用户 A 接下来将要交互的项目。

为了训练这样的模型，我们简单地将原始交互对象转换成顺序交互对象。剩下的都一样。

注意 to_sequence() 函数在长度不够长的序列前面填充零，以确保每个序列都具有相同的长度。

因此，id 为 0 的项目应该被更改为其他任意未使用的 id 号，以便此功能能够工作。

损失函数的选择

当指定模型时，我们可以灵活地改变损失函数。不同损失函数的模型在性能上可能有显著差异。我将简要描述 Spotlight 中定义的两种主要类型的损失函数。

• “逐点”:与其他形式的损失函数相比，这是最简单的形式。由于样本的稀疏性(效用矩阵中有很多 0)，将所有项目都考虑进去在计算上是不可行的。因此，不是计算给定用户的所有项目的损失，我们只考虑随机选择的一部分否定样本(用户没有与之交互的项目)和所有肯定样本。
• “bpr”:贝叶斯个性化排名(BPR)给每个用户的每个项目一个排名。它试图确保正样本的等级高于负样本的等级，公式如下。

ranking loss

现在，您已经了解了如何使用 Spotlight 构建推荐系统。它使用起来非常简单，并且有足够的灵活性来满足你的需求。虽然对于大多数问题，序列模型优于因式分解模型，但训练序列模型需要更长的时间。此外，在应用序列模型时，如果数据没有明确的序列相关性，也不会有很大帮助。

3.项目 2Vec

Item2Vec 的想法是上个月我参加一个竞赛时产生的，国际数据分析奥林匹克 (IDAO)。它提出的挑战要求参与者为 Yandex 建立一个推荐系统。因为那时我正在学习 Word2Vec，所以我想类似的概念也可以用在推荐系统中。我不确定是否有论文或文章指出过这个想法，但是我没有看到 Word2Vec 的概念在这个领域有类似的应用。Word2Vec 背后的大致思想是，我们利用分布式表示来编码每个单词。也就是说，每个单词由该单词周围的其他单词所确定的向量来表示。如果你想了解更多关于 Word2Vec 的内容，可以参考我之前的博文: Word2Vec 和 FastText 用 Gensim 进行单词嵌入。类似地，我试图使用分布式表示，根据用户在交互之前和之后交互的项目，对每个项目进行编码。

对于每个用户，我首先按照时间顺序创建了一个条目列表。然后，在这些项目列表上训练 Gensim 的 Word2Vec 模型。经过训练的项目向量存储在磁盘中，以便我们可以加载它供以后使用。

之后，我们将训练好的项目向量加载到嵌入矩阵中。

然后，我们定义了预测用户未来交互的模型。基本上，它是由 CuDNN 加速的 GRU 模型。如果你没有英伟达 GPU，也不用担心。你可以简单的把 CuDNNGRU 换成 GRU。

请注意，at 我将预训练的嵌入矩阵加载到模型的嵌入层中，并通过将 trainable 设置为 false 来冻结它。

4.表演

我已经在 IDAO 的公共电路板测试集上测试了我在本系列的两个部分中提到的一些模型。下表显示了它们各自的公开委员会分数。

Public Board Score

看起来基于神经网络的模型不一定能打败传统的推荐系统构建方法。虽然易于理解和实现，但 SVD 的得分与需要更长时间训练的 Spotlight 序列模型相当。令人惊讶的是，Item2Vec 是我测试过的所有型号中最好的型号。诚然，我们不能只用一个测试集来判断所有这些模型。这让你大致了解每个模型有多好。

5.结论

我已经讨论了在 Spotlight toolkit 中实现的两种类型的模型和我自己创建的依赖于 Word2Vec 概念的 Item2Vec。我还基于 IDAO 的公共电路板测试集比较了我提到的模型。事实证明，SVD 是推荐系统的有效解决方案，Item2Vec 已经证明了它能够更准确地推荐项目。如果你对这篇文章有任何问题，请不要犹豫，在下面留下你的评论或者给我发电子邮件:khuangaf@connect.ust.hk。如果你喜欢这篇文章，请确保你在 twitter 上关注我，以获得更多机器学习/深度学习博文！

递归神经网络导论

原文：https://towardsdatascience.com/introduction-to-recurrent-neural-network-27202c3945f3?source=collection_archive---------2-----------------------

有许多专门解决许多任务的深度学习模型。这里我们讨论深度学习模型处理序列的能力。

什么是顺序数据？

如果相关事物之间有特定的顺序，我们称之为序列。

《我是个好孩子》****《我是个好孩子吗》。

你认为两个句子意思相同吗？不！也就是说单词的位置很重要！它们是一系列单词。

想象一个正在播放的视频。如果你已经看过了，你可以很容易地预测下一个场景。但是考虑到你困了，不记得框架的位置了(脑子里全是冗杂的框架)。那你能预测下一个场景吗？？？当然不是！！！

递归神经网络

他们和我们一样！！！他们能记住序列。如果你让 RNN 预测下一个场景，它会告诉你。

当前馈神经网络仅基于当前输入做出决策时，RNN 基于当前和先前输入做出决策。

Recurrent Neural Network

绿框代表一个神经网络。箭头表示记忆或简单地反馈到下一个输入。

第一张图显示了 RNN。第二幅图展示了同一幅 RNN 在时间中的展开。考虑一个序列**【我是个好孩子】**。我们可以说序列是按时间排列的。**t = 0 时，X0 =“I”作为输入给出。t=1 时，X1 =“am”作为输入给出。**第一时间步的状态被记忆，并作为第二时间步的输入，与该时间步的当前输入一起给出。

在前向神经网络中，网络每个样本仅向前传播一次。但是在 RNN，网络正向传播等于每个样本的时间步数。

与 RNN 解决了不同的问题

生成文本

给定一个单词序列，我们希望在给定前一个单词的情况下预测每个单词的概率。

机器翻译

机器翻译类似于语言建模，因为我们的输入是源语言(例如德语)中的单词序列。我们希望输出目标语言(例如英语)的单词序列。

语音识别

给定来自声波的声学信号的输入序列，我们可以预测语音片段序列及其概率。

生成图像描述

与卷积神经网络一起，RNNs 已经被用作模型的一部分来为未标记的图像生成描述。

聊天机器人

聊天机器人可以回答你的问题。当给定一个单词序列作为输入时，将在输出端生成单词序列。

你需要知道的数据科学概念！第一部分

原文：https://towardsdatascience.com/introduction-to-statistics-e9d72d818745?source=collection_archive---------0-----------------------

这是 5 篇系列文章的第一篇，将概述数据科学中的一些核心概念:

值得注意的是，我们将涵盖上图中的统计桶，并试图使其尽可能直观和数学化。

我的背景是南部和东部非洲的大型制药公司、学术界(牛津大学物理和理论化学系的博士)和“为社会公益服务的数据科学”行业。

这一系列帖子旨在介绍并快速发展数据科学和数据分析中的一些核心概念，并特别关注我觉得在其他材料中被忽略或简单处理的领域。因此，这些帖子将适合希望重温统计理论的数据科学家，希望获得全面培训材料的数据分析师，以及希望从数据和员工那里提出更好问题的数据经理。

我们将在本系列中讨论的一些主题包括:

• 功率和样本量计算的非参数方法
• A/B 测试(又名随机试验)的设计和分析
• 如何有效探索数据趋势
• 如何实施、评估和改进线性模型

这些帖子改编自我为国际非政府组织 One Acre Fund 开发的一系列培训材料，One Acre Fund 为东非的小农提供小额贷款。

这篇文章的目的(#1):

• 理解是什么决定了分布的形状。
• 理解非正态性如何影响假设检验和功效计算
• 了解如何使用蒙特卡罗方法进行非参数样本容量/功效计算

这些帖子是用 R-markdown 制作的，它允许在一个地方呈现 R 代码、文本写作和结果。R-markdown 对于编写报告特别有用，因为它使分析透明且可重复。

动机:

统计学本质上是对分布的研究。分布是将一个值的频率与观察到的实际值联系在一起的一种方式。你可以想象一下，例如，如果我们在测量员工的身高，我们可能会发现大多数人的身高在 1.5 到 1.8 米之间，但也会有一些人比这更高或更矮。当我们谈论这样的数据时，我们实际上是在谈论分布:

set.seed(111)
hist(rnorm(1000, mean=1.75, sd=0.5), xlab="Employee height (m)", main="Employee height", breaks=c(0,0.5,1,1.5,2,2.5,3,3.5,4),col="black")

上面的柱状图是我们用图形的方式表述的；我们有沿 X 轴的观察值(这些值被分组为“箱”)和沿 Y 轴的计数或频率(即我们观察该值的次数)。有时我们也可以绘制沿 Y 轴的观察的概率。

阅读上面的直方图，我们可以说最常见的身高是 1.5 到 2 米，在更极端的值下频率会降低(例如，0.5-1 米和 2.5-3 米的人比 1.5-2 米的人少)。更具体地说，通过阅读 X~Y 轴，我们可以说大约有 400 人身高 1.5 到 2 米，大约有 200 人身高 2 到 2.5 米。

每一项分析都必须从理解底层数据及其分布开始。分布的形状将决定:

• 我们使用的分析方法(我们应该使用 T 检验，还是韦尔奇检验？)
• 样本量(我们需要多少个人才能发现一个效应？)
• 我们的误差大小(数据有多分散？)
• 统计显著性(即干预/产品是否成功，其效果有多大)

最终，这种分布将决定我们是否能为我们测试的干预赋予意义和因果关系。

正态分布

上面的直方图来自正态分布，这是你将看到的最常见的分布之一。让我们更深入地研究一下这个分布。在正态分布中(也只有正态分布),我们期望我们的均值和中值位于分布的中心(最高点),我们同样期望分布关于均值大致对称，如下所示:

#lets generate some example normal distributions 
set.seed(111)# by setting a seed we can make sure that R draws random numbers reproducibly
#simulate a random normal distribution
norm1 <- rnorm(20000, mean=0,sd=1)

normdf <- data.frame("val"=norm1, "flag"=0)

#which is a logical, it will evaulate a condition: here "which val is greater than 1sd of norm OR less than -1sd of norm"
# note that | means OR and & means and
sd2 <- which( normdf$val > 1* sd(norm1) | normdf$val <  -1*sd(norm1) )
sd3 <- which( normdf$val > 2*sd(norm1) | normdf$val <  -2*sd(norm1) )
sd1 <- which( normdf$val < 1*sd(norm1) & normdf$val > -1*sd(norm1) )

#now lets use the index made above to reset the values in flag to the following
normdf$flag[sd1] <- "Within 1 SD"
normdf$flag[sd2] <- "Within 2 SD"
normdf$flag[sd3] <- "Within 3 SD"
normdf$flag <- as.factor(normdf$flag)

#lets plot these...
#as a histogram
ggplot(normdf, aes(x=val, fill=flag)) + geom_histogram(binwidth=.5, bins=100, alpha=1, breaks=c(-Inf,min(subset(normdf, flag=="Within 3 SD")$val),min(subset(normdf, flag=="Within 2 SD")$val),min(subset(normdf, flag=="Within 1 SD")$val), max(subset(normdf, flag=="Within 1 SD")$val), max(subset(normdf, flag=="Within 2 SD")$val),max(subset(normdf, flag=="Within 3 SD")$val), Inf)) + geom_vline(aes(xintercept = mean(val), colour="Mean & Median"),color='black',size=2) + geom_text(aes(x=mean(norm1)-0.6,y=15000,label="Mean & \n Median")) + xlab("Value of observation") + ylab("Number of observations")

正态分布由两个指标定义，均值和标准差。对于正态分布，我们预计 68%的值位于 1 个标准差内，95%位于 2 个标准差内，99.7%位于 3 个标准差内(如上)。

请注意，分布的标准偏差与分布的宽度成比例。

当我们展示分布时，我们的目标是使图表尽可能简单，同时仍然忠实于数据。

呈现数据

我们可以用其他格式呈现分布，而不仅仅是直方图:

#
#histo
p1 <- ggplot(normdf, aes(x=val)) +  geom_histogram(binwidth=.5, colour="black", fill="white") + ggtitle("Histogram")

#density
p2 <-ggplot(normdf, aes(x=val)) + geom_density() + ggtitle("Density plot") 

#boxplot
#the theme argument removes the x-axis as we have no information to show here
p3 <- ggplot(normdf, aes(x=1, y=val)) + geom_boxplot() + ggtitle("Boxplot") + xlab("")  + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
#p3 <- ggdraw(add_sub(p3, "Boxplot"))

#finally a beeswarm plot
library(ggbeeswarm)

p4 <- ggplot(normdf) +  geom_beeswarm(aes(x=1,y=val)) + ggtitle("Beeswarm") + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

#this is a funtion defined in the first block which allows us to put more than one graph on the same page
multiplot(p1,p2,p3,p4, cols=2)

我们以不同的方式展示了上述完全相同的数据:

• 我们熟悉的直方图(左上)。请注意，改变直方图的条柱数量(即我们需要多少条柱)会极大地改变我们对数据的理解。
• 盒状图(右上)在 Y 轴上有我们的值，粗中心线代表中间值，上下盒边界代表第一和第三个四分位数(分别)，从盒中伸出的“胡须”代表最大值和最小值，最后点是(假定的)异常值。箱形图信息丰富，但可能隐藏聚类数据。
• 密度图类似于我们的直方图。但是，密度图被平滑，y 轴现在表示总计数的分数(从 0 到 1 ),而不是实际计数。密度图下的面积总和为 1。
• 蜂群图显示了数据框中的每个数据点。图的宽度与观察值的数量成正比，我们可以知道分布的正态性，因为图在 0(我们的平均值)附近较厚，然后在平均值的北部和南部对称地变小。虽然箱线图可以隐藏聚类，密度和直方图依赖于良好的箱宽度，但蜂群可以清楚地显示聚类数据，但看起来有点难看。

为了说明这一点，让我们展示一些聚类数据(也称为二项式分布，注意这不同于正态分布):

#
#lets make two seperate distributions and join them:
set.seed(122)
nor1 <- rnorm(200, mean=0,sd=0.8)
nor2 <- rnorm(50, mean=4,sd=0.7)

clusterdf <- data.frame("val"=c(nor1,nor2), "flag"=0)

# a beeswarm plot

p4 <- ggplot(clusterdf) +  geom_beeswarm(aes(x=1,y=val)) + ggtitle("Beeswarm") + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

#histo
p1 <- ggplot(clusterdf, aes(x=val)) +  geom_histogram(binwidth=.5, colour="black", fill="white") + ggtitle("Histogram")

#density
p2 <-ggplot(clusterdf, aes(x=val)) + geom_density() + ggtitle("Density plot") 

#boxplot
#the theme argument removes the x-axis as we have no information to show here
p3 <- ggplot(clusterdf, aes(x=1, y=val)) + geom_boxplot() + ggtitle("Boxplot") + xlab("")  + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
#p3 <- ggdraw(add_sub(p3, "Boxplot"))

#this is a funtion defined in the first block which allows us to put more than one graph on the same page
multiplot(p4,p1,p2,p3, cols=2)

从蜂群中我们可以看到，在~0 附近有一个很大的数据集群，但在~4 附近也有一个较小的集群。密度图和直方图显示了这一点，但如果我们只检查箱形图，我们将无法清楚地看到这一点。

最后，改变密度图或直方图的条块也可以改变我们对数据的解释:

#density
p1 <-ggplot(clusterdf, aes(x=val)) + geom_density(bw=0.05) + ggtitle("Density plot 1") 
p2 <-ggplot(clusterdf, aes(x=val)) + geom_density(bw=0.2) + ggtitle("Density plot 2") 
p3 <-ggplot(clusterdf, aes(x=val)) + geom_density(bw=2) + ggtitle("Density plot 3") 
p4 <-ggplot(clusterdf, aes(x=val)) + geom_density(bw=6) + ggtitle("Density plot 4") 

multiplot(p1,p3,p2, p4, cols=2)

即使数据是相同的，我们对数据的解释也会根据我们检查的图而发生巨大的变化！我们如何从图 4 中了解我们的数据？图 1 可能太详细，而图 3 和图 4 太不详细(注意图 4 完全隐藏了第二组数据)。

剧情 2 既简单又符合数据。因此，这是这里最好的情节。请注意，找出最佳图的唯一方法是检查基础数据。因此，我建议在制作数据的简化版本之前，使用蜜蜂群或具有小区间宽度(或高区间数)的密度图，首先尽可能详细地检查数据。

有效地比较分布

通常我们会想要比较分布，看看是否有明显的差异，例如当比较治疗组和对照组时。我的首选方法是使用密度图或箱线图。让我们看看下面的一些例子，在这些例子中，我们测量了一种处理方法对卢旺达农民采用化肥的影响:

 #lets make some data
set.seed(122)
nor1 <- data.frame("val"=rnorm(50, mean=0.6,sd=0.15), "Group"="Control")
nor2 <- data.frame("val"=rnorm(50, mean=0.8,sd=0.1), "Group"="Treated")

#bind the data by rows (e.g. stack one ontop of the other)
dat <- rbind(nor1,nor2)

#Density plots with semi-transparent fill
g1 <- ggplot(dat, aes(x=val, fill=Group)) + geom_density(alpha=.3, bw=0.1) + xlab("Fertilizer adoption")

# boxplot and remove redundant legend. 
g2 <- ggplot(dat, aes(x=Group, y=val, fill=Group)) + geom_boxplot() +
    guides(fill=FALSE) + ylab("Fertilizer adoption")

g3 <- ggplot(dat, aes(x=val, fill=Group)) + geom_histogram(alpha=.3, bw=0.1) + xlab("Fertilizer adoption")

g4 <- ggplot(dat, aes(x=val, fill=Group)) + geom_beeswarm(aes(x=0, y=val, colour=Group)) + ylab("Fertilizer adoption") + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

multiplot(g1,g2,g3,g4, cols=2)

我们可以看到，密度图和箱线图对于观察组间差异非常有用。然而，蜂群图和直方图过于繁忙，无法快速得出结论(尽管它们在其他场景中很有用，例如查看异常值)。因此，这些数据应该以密度图或箱线图的形式呈现，以满足我们的标准，即易读但真实的数据。

显著和非显著差异

所以我们现在有了治疗组和对照组的描述图。

ggplot(dat, aes(x=val, fill=Group)) + geom_density(alpha=.3, bw=0.1) + xlab("Fertilizer adoption")

这两组看起来不同，但是我们怎么能更确定，我们怎么能量化我们的不确定性呢？

假设检验

统计测试的存在是因为我们只有一个人口的样本(见词汇表)。例如，我们可以在随机试验中测量东非农民样本的肥料使用情况。但是我们正试图推断出一些关于所有东非农民的情况(即从我们的样本中归纳出全部人口)。我们希望能够了解对实际人群的影响，这样我们就可以对项目干预做出决定；因此我们想要比较人口，但是我们只有样本。假设检验允许我们在给定一些关键假设的情况下检验潜在人群的差异。

假设检验有某些潜在的假设。常用测试的一些示例假设如下:

• T 检验:假设正态分布。该测试是在 20 世纪 50 年代针对小样本量进行的，现在有了更可靠的测试，应该使用
• 配对 T 检验:T 检验(同上)，但用于测量治疗前后的同一受试者(例如，在一次散发传单活动前后测量农民种植依从性)
• 双样本 T 检验:用于两个不同组的 T 检验，例如比较接受治疗的农民和对照农民的玉米产量。
• **Welch 检验:**对 T 检验的改进，对不等方差更稳健(例如，对密度图宽度不同的样本更稳健)。我们可以在 t.test() R 函数中调用带有参数“var.equal=FALSE”的 Welch 测试。请注意，韦尔奇检验仅对正态分布数据有效(如 T 检验)。

我们如何知道何时使用每个测试？我们如何检验这些测试的假设？

常态假设

许多假设检验和样本量计算都假设数据是正态分布的(记住，这意味着数据关于平均值是对称的)。然而，我们经常希望对非正态数据进行假设检验。我们将很快检查非正态数据的假设检验，但是首先，我们如何知道我们的数据是否正态分布？这里有许多方法，我将介绍两种比较常见的方法。

让我们来看看下面的数据集:

#
#lets make some data
set.seed(122)
dat <- data.frame("val"=rnorm(100, mean=0.6,sd=0.15))

#this lets us plot graphs next to each other
par(mfrow = c(1, 2))

#this is a simple histogram
# i call the column named "val" from dataframe named "dat" with the $val addition
hist(dat$val, main="Data A", xlab="values")

hist(log(dat$val), main="Data B",xlab="values")

希望很明显，数据 A 可能是正常的，但数据 B 看起来肯定是非正常的(注意明显的不对称)！

为了更彻底地检验正态性，我们可以绘制一个理想的正态分布，然后将我们的数据放在上面(这被称为 QQ 图)，然后计算关系之间的 R 平方:

#
#plot data against ideal normal
#note that the subset command means we will only take the val column of data with Group=="Control"
# subset(dat, Group=="Treated") would return the whole dataframe where the column "Group" is equal to "Treated", adding the $val allows us to access the val column directly
par(mfrow = c(1, 2))

qqnorm(dat$val, main="QQ plot A")
qqline(dat$val)

qqnorm(log(dat$val), main="QQ plot B")
qqline(log(dat$val))

#get Rsquared
qn=qqnorm(log(dat$val), plot.it=FALSE)
rsq_B <- cor(qn$x,qn$y)

#get Rsquared
qn=qqnorm(dat$val, plot.it=FALSE)
rsq_A <- cor(qn$x,qn$y)

QQ 图将完美的正态分布绘制为一条实对角线，然后将我们的实际数据作为标记(圆圈)放在其上。

从图中我们可以看出，我们的数据 A 很好地遵循了正态(实斜线)分布。然而，检查图 B 中的数据，我们可以看到与正常线的实质性偏差，这些偏差显然不是随机的，而是遵循一个趋势(随机偏差在线的任一侧出现的次数相同)。

为了量化这一点，我们可以提取 R 平方(我们将在后面的课程中了解更多关于 R 平方的内容，但它本质上是在 0 到 1 的范围内测量我们的点与对角线的接近程度)。A 和 B 的 R 平方值分别为 0.997 和 0.974。查看数据就足以认识到曲线 B 是非正态的，但是我们也可以选择一个截止值，例如 R 的平方为 0.975，并将低于这个值的任何值称为非正态，将高于这个值的任何值称为正态。

另一种方法是使用夏皮罗检验，它将对我们的数据是否正常进行假设检验:

shapiro.test(dat$val)

这里的关键部分是 P 值，数据 a 的 P 值为 0.84 值越高，分布越正态。通常我们会使用 0.1 的 P 值临界值来确保正态性。因此，任何 p 值为<0.1 we would class as non-normal. This approach is complementary to the above, and I would recommend both plotting a QQ-plot (as above) 和的数据都要进行夏皮罗检验(并在任何报告中报告这两个值)。

让我们看看一些非正态(特别是对数正态)数据的 QQ 图和夏皮罗测试:

non_norm <- log(rnorm(75, mean=5,sd=3))

par(mfrow = c(1, 2))
hist(non_norm, main="Histogram of log-normal data", xlab="Values")
qqnorm(non_norm, main="QQ plot for log-normal data")
qqline(non_norm)

上述分布的夏皮罗测试 P 值为 2.4e-09，QQ 图显示了与正常线的系统偏差。这强烈表明我们的数据 B 是非正态的。

假设检验正态分布

一旦我们确信我们的数据是正态分布的(我们稍后将处理非正态数据)，我们就可以进行一些假设检验。T-test 通常被用作“go-to”假设检验，但这里我想说服您使用 Welch 检验(记住，要访问 R 中的 Welch 检验，我们调用 t.test()函数，但添加参数 var.equal=FALSE )。韦尔奇检验更适合于我们在现实世界中经常遇到的不等样本量和不等方差的情况。

下面，我运行了一个模拟来演示在方差数据不相等的情况下使用 T-test 优于 Welch tests 的危险性:

#set up vars
nSims <- 20000 #number of simulations
p1 <-c()
p2 <- c()

#create variables for dataframe
catx<-rep("x",38)
caty<-rep("y",22)
condition<- c(catx,caty)

#run simulations
for(i in 1:nSims){ #for each simulated experiment
  sim_a<-rnorm(n = 38, mean = 0, sd = 1.11) #simulate participants condition a
  sim_b<-rnorm(n = 22, mean = 0, sd = 1.84) #simulate participants condition b
  p1[i]<-t.test(sim_a,sim_b, alternative = "two.sided", var.equal = TRUE)$p.value #perform the t-test and store p-value
  p2[i]<-t.test(sim_a,sim_b, alternative = "two.sided", var.equal = FALSE)$p.value #perform the Welch test and store p-value
}

par(mfrow = c(1, 2))
hist(p1, main="Histogram of t-test p-values ", xlab=("Observed p-value"))
hist(p2, main="Histogram of Welch's p-values", xlab=("Observed p-value"))

我们在这里绘制了多次运行模拟后从两次测试中得出的 P 值直方图。我们会记得 p 值是组间差异纯粹是偶然观察到的概率(即组间没有真正的差异)。这里的两个总体具有相同的均值，这意味着我们应该接受零假设，拒绝替代假设(见术语表)。

正如我们所知，人群之间没有真正的差异，我们应该接受零假设，并看到一个平坦的 p 值分布(我们应该看到一个 p 值为<0.05 5% of the time by definition, as a p-value of 0.05 implies a 5% chance of a false result in frequentist statistics). Note that in a real situation we wouldn’t already know if there was or was not a difference between populations!

We see that the right hand histogram of the Welch test results is flat (this is good, it matches what we know about the data). However, the T-test (left hand histogram) performs very poorly here and reports more smaller p-values. This means the chances of us wrongly concluding that 那里 是 )组之间的显著差异要高得多。准确地说，如果我们在这里使用 T 检验而不是韦尔奇检验，我们有 72 %的可能做出不正确的结论！！！这可能导致接受实际上没有效果的程序或产品创新，浪费时间和金钱。

这个故事的寓意是在正确的时间使用正确的测试。如果数据是正常的，那么我们应该总是使用韦尔奇检验(忽略 T 检验)，因为韦尔奇检验比 T 检验有更少的假设。如果数据是非正常的，那么我们需要一个非常不同的测试。

假设检验非正态分布

以上适用于正态分布的数据。但是非正态数据呢？

在现实生活中，我们经常会遇到非正态数据。任何我们的值接近 0 或 100%的情况都可能是非正常的。例如，如果我们观察一家老牌银行的贷款偿还情况，或者如果我们观察一种非常受欢迎(或不受欢迎)的产品的采用情况。

卢旺达农民的非正态分布示例如下(贷款偿还百分比和太阳能灯采用率):

df <- read.csv("./Sample_data.csv")
df$X..Repaid <- df$TotalRepaid/df$TotalCredit * 100
df$solar_A <- df$solar_A *100
d1 <- ggplot(df, aes(x=X..Repaid)) + geom_density(alpha=0.3, fill="red") + ggtitle("% repaid density plot")  + xlab("% repaid") + xlim(75,100)
d2<- ggplot(df, aes(x=solar_A)) + geom_density(alpha=0.3, fill="red") + ggtitle("Solar adoption density plot")  + xlab("% solar adoption") 

multiplot(d1,d2,cols=2)

我们可以从密度图中清楚地看到，我们的数据在两种情况下都是非正态的。但是让我们使用上面概述的诊断来确保:

par(mfrow = c(1, 2)) qqnorm(df$X..Repaid, main="QQ plot for % repaid") qqline(df$X..Repaid) qqnorm(df$solar_A, main="QQ plot for solar adoption") qqline(df$solar_A)#Shapiro test has an upper limit of 5000, so we will randomly take a subset of data to test instead shapiro.test(sample(df$X..Repaid, size=1000))Shapiro-Wilk normality test data: sample(df$X..Repaid, size = 1000) W = 0.54519, p-value < 2.2e-16shapiro.test(sample(df$solar_A, size=1000))Shapiro-Wilk normality test data: sample(df$solar_A, size = 1000) W = 0.87404, p-value < 2.2e-16

哇！这些是一些看起来很糟糕的 QQ 图，和一些非常低的 P 值！数据远非正常。

然而，我们仍然希望对这些数据进行假设检验。

如何检验非正态数据

现在让我们比较一下测试增加 1.5%太阳能利用率的干预措施的方法。我们将使用 T 检验(对非正态数据不正确)、韦尔奇检验(对非正态数据也不正确)和威尔科森检验(对非正态数据的适当检验)。

#load data
norm1 <- df$solar_A
norm2 <- norm1 *1.015

T 检验:

#t test
t.test(norm1,norm2, alternative = "two.sided", var.equal = TRUE)

P 值= 0.1048

韦尔奇试验:

#welch test
t.test(norm1,norm2, alternative = "two.sided", var.equal = FALSE)

P 值= 0.1048

威尔科克森试验:

wilcox.test(norm1 , norm2 ,paired=FALSE, correct=TRUE)

P 值= 1.55e-06

我们可以看到我们得到了非常不同的结果！T 检验和 Welch 检验报告 p 值为 0.105 和 0.105，这表明在典型的 p 值阈值 0.05 下无显著性。而 Wilcoxon 试验表明 p 值为 1.55e-06。

你可以看到，根据我们使用的测试，我们会得出完全不同的结果。T 检验和 Welch 检验会使我们相信干预是没有用的，而 Wilcoxon 检验清楚地表明干预似乎有效果。因此，对正确的数据使用正确的测试非常重要。想象一下，如果我们一直在测试一种对农民有害的干预措施，但使用了错误的测试方法，因此得出干预措施是安全的结论！

我们现在可以将 Wilcoxon 测试添加到我们的库中:

• Wilcoxon-test :对假设很少的非正态数据的良好测试。我们正在测试两个数据集具有不同分布的假设。这里的主要假设是:1)样本随机代表总体，2)样本相互独立，3)值有顺序(如 3 大于 1，但不能说真大于假)。

报告数据

我们现在已经了解了如何计算数据的稳健 P 值。然而，我们报告数据的方式也取决于分布的形状。对于正常数据，我们可以简单地报告平均值和置信区间，例如平均值为每英亩 4.5(CI:1.5–7.5)千克(您也可以报告带有标准误差的平均值，只要您报告的内容清晰明了)。非正态数据有点复杂，我们通常应该报告中位数、第一和第三四分位数，例如中位数 4.5 (Q1 2.5 & Q3 6) 。

我们将记住，正态分布是对称的，正是非正态分布的不对称性阻止了我们报告 CIs、SEMs 或 StDevs(因为这些都是基于对称性的)。

与上述假设检验思想相关的一个概念是统计功效。幂用于计算样本量，也用于解释试验后 RCT 的结果。低功率意味着我们的结果不太可靠，我们可能会对干预做出错误的结论。

我们使用正确的功率计算非常重要，就像我们使用正确的假设检验一样重要。功率计算取决于数据分布，R 中最常用的功率计算器适用于正态分布的数据。这个计算器需要一个“效果大小”，这是一个简单的方法来衡量治疗的效果。效应大小基本上是样本间平均值的差异除以样本的标准差。这意味着单位是标准偏差。因此，效应大小为 1 意味着数据移动了等于 1 个标准差的量。我在下面加入了一个函数，使用科恩的 D 方程计算效应大小:

cohen_d <- function(d1,d2) {  

  m1 <- mean(d1, na.rm=TRUE)
  m2 <- mean(d2, na.rm=TRUE)
  s1 <- sd(d1, na.rm=TRUE)
  s2 <- sd(d2, na.rm=TRUE)
  spo <- sqrt((s1**2 + s2**2)/2)
  d <- (m1 - m2)/spo
  rpb <- d / sqrt((d**2)+4)
  ret <- list("rpb" = rpb, "effectsi" = d)
  return(ret)  }

现在我们有了一种计算效果大小的方法，让我们用它来计算一项有正常数据和 100 名参与者的研究的功效。我们将使用 pwr 库中的函数“pwr.t.test”。该函数需要以下参数:

• h —这是由上面的科恩 D 方程计算出的效应大小。我们用对象后面的“$effectsi”来调用它(详见代码)。
• n——每个样本的样本量
• sig.level —这是我们希望在分析中使用的 p 值临界值。这通常是 0.05 或 0.1。
• power —我们希望将此设置为 NULL，以获得将电源返回给我们的功能。

我们一起将这个库称为“pwr”。

其中“length()”参数返回名为“norm”的列表中元素的数量。

norm <- rnorm(200,40,10)
norm2 <- norm*1.1
inp <- cohen_d(norm, norm2)

这里我们看到我们有一个 0.97 的幂。理想情况下，我们希望功效高于 0.75–0.8，以获得准确检测效果的良好概率。

请注意，我们可以通过填写幂参数并将 n 设置为空来计算 0.8 的幂所需的样本量，如下所示(这是一个略有不同的等式，因为我们假设对照组和治疗组的大小相等):

*pwr.t.test(n= NULL, d= inp$effectsi, power=0.8, alternative="two.sided")*

我们可以循环查看不同样本大小的功效，然后用一条线标出功效为 0.8 的位置。

power_out <- c() # initialise this empty so we can append to it
sample_sizes <- c(2,4,6,10,20,30,40,50,70,90,110,130)
for(i in sample_sizes){
  x <- pwr.t.test(n= i, d= inp$effectsi, power=NULL, alternative="two.sided")$power
  power_out <- c(power_out, x)

}

p <- ggplot() + geom_line(aes(x=sample_sizes, y= power_out)) + geom_hline(yintercept = 0.8, color="red", size=2) + ggtitle("Example power curve")
p

像上面这样的图可以帮助你了解什么样的样本大小和功效是合理的。我强烈建议绘制功效曲线的下限(即你可能看到的最小功效)和上限(即你可能看到的最大功效)。

运行一个完整的例子:让我们想象我们想要运行一个 RCT(或 A/B 测试),检查一个新的利率对肥料使用的影响。我们有两组，一组是对照组，利率和之前一样。以及利率降低 3 个点(例如 15%的利息而不是 18%)的治疗组。我们对客户的了解意味着，我们认为通过这种处理方式，化肥的使用量将增加 15%至 25%。我们绘制如下的功率曲线*(请检查代码并确保您理解了代码)*:

control <- rnorm(100, 0.5,0.2)
treat_lower_estimate <- control *1.15
treat_upper_estimate <- control *1.25

power_lw <- c() # initialise this empty so we can append to it
power_hi <- c() # initialise this empty so we can append to it
sample_sizes <- c(10,20,30,40,50,70,90,100,125,150,175,200,250,300,350,400,450,500,550,600,700)
for(i in sample_sizes){

  lower_cohen <- cohen_d(control, treat_lower_estimate)
  a <- pwr.t.test(d = lower_cohen$effectsi , n=i,  sig.level = 0.05, power = NULL)$power
  power_lw <- c(power_lw, a)

  upper_cohen <- cohen_d(control, treat_upper_estimate)
  b <- pwr.t.test(d = upper_cohen$effectsi , n=i,  sig.level = 0.05, power = NULL)$power
  power_hi <- c(power_hi, b)

}

marker <- pwr.t.test(d = lower_cohen$effectsi , n=NULL,  sig.level = 0.05, power = 0.8)$n
marker2 <- pwr.t.test(d = upper_cohen$effectsi , n=NULL,  sig.level = 0.05, power = 0.8)$n
ggplot() + geom_ribbon(aes(x=sample_sizes, ymin= power_lw, ymax=power_hi), alpha=0.2, colour="blue", fill="blue")  + xlab("Sample size") + ylab("Power") + geom_vline(xintercept = marker, linetype="dotted" ) + geom_hline(yintercept=0.8, linetype="dotted", size=2) + geom_vline(xintercept = marker2 , linetype="dotted") + labs(title="Power curve example", caption="Power curve indicating sample sizes needed for a 0.8 power (dotted lines)") + theme(plot.caption = element_text(hjust = 0.5))

上面的图让我们看到它需要一个大的样本量(n 大约是。200，所以我们在这个试验中需要 400 个参与者——200 个控制组和 200 个治疗组)来检测 15%的采用变化，其功效为 0.8。像这样绘制数据允许更灵活的决策，强烈建议这样做。

非正态分布数据和样本大小/功效

对于正常的功率计算来说，以上都没问题。但是非正常力量呢？

我们对此也有解决方案！

我在下面做了一个函数(“MCpower”)，它将计算非正态分布的功率。它需要样本 1 和 2(处理和对照)的数据(或预期数据)和样本大小(“大小”参数)。

请确保你理解这个代码，并舒适地使用它

MCpower = function(sample1, sample2, size) {

  reps = 1000
    results  <- sapply(1:reps, function(r) {
        resample1 <- sample(sample1, size=size, replace=TRUE) 
        resample2 <- sample(sample2, size=size, replace=TRUE) 
        test <- wilcox.test(resample1, resample2, alternative="two.sided",paired=FALSE, correct=TRUE)
        test$p.value
    })
    sum(results<0.05)/reps
}

#we use it like this
Non_norm_power <- MCpower(control, treat_upper_estimate, 100)

我们可以比较这些方法来理解为什么我们需要对正确的数据使用正确的测试。让我们进行一些模拟实验。下面的数据是非正态，样本之间是有差异的。因此，我们期望看到统计上的显著差异。请注意，虽然正常的功效测试或假设测试将测试均值，但这里的功效计算将测试分布的总体相似性。

让我们运行一个模拟，计算用不同方法测试两个总体之间的差异所需的功效:

#set up vars
nSims <- 500 #number of simulations

#initialise empty vectors for appending to
x1 <- c()
x2 <- c()
p1 <- c()
p2 <- c()

#run simulations
for(i in 1:nSims){ #for each simulated experiment
  set.seed(i)
  sim_a <- c(rnorm(150, mean=130,sd=20) )**2
  sim_b <- c(rnorm(150, mean=120,sd=35) )**2
  sim_b[which(is.nan(sim_b)==TRUE)] <- 1
  sim_a[which(is.nan(sim_a)==TRUE)] <- 1
  inp <- cohen_d(sim_a, sim_b) # effect size
  x1[i] <- pwr.t.test(d= inp$effectsi , n =length(sim_a), sig.level = 0.05, power = NULL)$power
  x2[i] <- MCpower(sim_a, sim_b,  size=150) 
  p1[i] <- t.test(sim_a,sim_b)$p.value
  p2[i] <- wilcox.test(sim_a,sim_b, paired=FALSE, correct=TRUE)$p.value
}

#hist(sim_b, bins=20)
#shapiro.test(sim_b)
par(mfrow = c(1, 2))
hist(x1, main="Histogram of normal powers ", xlab=("Observed power"), xlim =c(0,1))
hist(x2, main="Histogram of MC powers", xlab=("Observed power"), xlim =c(0,1))

数百次 t 检验模拟的平均 p 值为 0.133，MC 检验的平均 p 值为 0.0305(参见上述分布直方图)。我们现在知道 T-检验假设正态分布，这在这里是不正确的，而 MC 方法只需要满足很少的假设。我们可以看到，MC 方法通常具有更高的功效，并且经常正确地识别出这些样本之间存在和显著差异，不像 T 检验表现不佳。

注意，对于正常数据，正常功率计算的结果和我的 MC 方法非常相似。比较效应大小为 1 的两种分布，标准方法和 MC 方法的 P 值分别为 0.012 和 0.014，功效分别为 0.86 和 0.85。

摘要

我们现在已经了解了如何绘制和呈现数据，以使其简单且真实地反映基础数据，我们还了解了如何计算 RCT 中正常和非正常数据的 p 值和乘方，以及当我们的假设有误时会发生什么。我附上了一个简短的流程图来说明需要作出的决定:

分析 RCT 的最佳实践:

• 即使您在 RCT 前计算样本量，计算 RCT 后的功效(与样本量相关)也很重要，以确保我们的研究有足够的功效。请记住，低功率意味着我们不太可能检测到干预的效果，通常会使我们的结果不太可靠。我希望每一次 RCT 分析都从关键变量的功效计算开始
• 我们还必须绘制关键变量的分布图。在 RCT 中，我们通常会观察 2 个或更多组(例如对照组和治疗组)。因此，我建议在分析开始时为每个关键变量绘制密度图或箱线图。R-Markdown(这是本课的内容)在这里非常有用，将成为 AMP 的重点。
• 检验我们的假设。在决定使用什么方法之前，我们应该对我们的数据进行 Shapiro 测试，以确保它是正常的*(例如 T-测试或 Wilcoxon 测试)。*
• 对于正态分布数据，始终报告 p 值的平均值和置信区间，对于非正态数据，报告中位数和第一个及第三个四分位数。
• 确保图表简单，但真实的数据。

统计测试概述

• T 检验:假设正态分布。该测试是在 20 世纪 50 年代针对小样本量进行的，现在有了更可靠的测试，应该使用
• 配对 T-检验:T-检验(如上)，但用于测量治疗前后的同一受试者(例如，测量发传单活动前后的农民种植依从性)
• 双样本 T 检验:用于两个不同组的 T 检验，例如比较接受治疗的农民和对照农民的玉米产量。
• **Welch 检验:**对 T 检验的改进，对不等方差更稳健(例如，对密度图宽度不同的样本更稳健)。我们可以在 t.test() R 函数中调用带有参数“var.equal=FALSE”的 Welch 测试。请注意，韦尔奇检验仅对正态分布数据有效(如 T 检验)。
• Wilcoxon-test :对假设很少的非正态数据的良好测试。我们正在测试两个数据集具有不同分布的假设。这里的主要假设是:1)样本随机代表总体，2)样本相互独立，3)值有顺序(如 3 大于 1，但不能说真大于假)。

资源

• 设置/安装说明 R
• 互动 R 基础教程——强烈推荐
• Base R cheatsheet —强烈推荐，打印出来放在办公桌上！
• 非互动教程(如果那是你的事)
• Youtube 教程
• 备忘单 —仅推荐给更高级的 R 用户
• StackOverFlow —针对具体的技术问题

词汇表

• 观察:在统计学中，这是一个变量在特定时间的值。我们也可以把观察看作数据帧中的数据行
• 变量:可以观察到的特性。这也可以被称为一个特性，可以被认为是数据帧中的一列(例如:性别、糖价或收成)。
• 正态:正态分布是近似对称的，均值一定等于中位数。数据点离平均值越远，发生的可能性就越小。正常数据也称为“钟形数据”。
• 方差和标准差:方差是正态分布的宽度。标准差是方差的平方根。我们可以认为这是数据的分布，高方差意味着大量的可变性，低方差意味着大多数观测值将具有非常相似的值。
• 四分位数:四分位数是另一种数据度量。如果我将一个变量的所有值按顺序排列，那么第一个四分位数就是总长度的四分之一。同样，第三个四分位数将是四分之三长度的值。例如，数据 1，3，6，8，10，14，15，18，20，22，25，26 的第一个四分位数为 6 ，第三个四分位数为 20 。
• 人口:从大量潜在数据中抽取样本。
• 样本:我们有数据的总体的子集。例如，我们可能正在研究一个由 1000 名肯尼亚农民组成的样本，但是想要推广到整个肯尼亚农民人口。
• P 值:观察结果完全随机产生的概率(即群体之间没有真正差异的概率)
• 功效:在给定样本量、分布和误差率的情况下，我们能够检测出两个总体之间差异的概率。请注意，功率也会影响我们的结果的稳健程度。
• 效果大小:2 组之间的标准化差异。它有单位 stdev 。因此，1 的效应大小等于 1 stdev 的差值。群体之间。
• 假设检验:统计检验。每个假设检验有两个部分，一个是零假设，表示两个群体是相同的(即没有差异)，另一个是替代假设，表示两个群体之间有差异。通常，零假设被称为 H0，替代假设被称为 H1。
• 标准形式/科学记数法:我们可以用标准形式写出有很多零的数字。例如 0.001 变成 1e-3，1000 变成 1e+3(意思是 1 *10 的 3 次方)。
• 直方图:矩形图(也称为条柱)，条柱的高度与观察值的数量成比例，宽度与包含的值的范围成比例。

下一篇文章:

• 最小可检测影响
• 计划 A/B 测试
• 集群内效应
• 还有更多！

最初发布于Michael-bar . github . io。

如果你在这些帖子中看到任何错误，请留下评论！同样，如果任何事情需要进一步的解释，然后评论！