python随机数用法实例 python 随机数 原理

前言

还在为生成随机数而头疼吗？还为生成随机数不够随机而苦恼吗？学了二项式分布不知道在实际生活中有啥用处吗？如何用用随机数来观察正态分布？

读下去，下面，就是你要的答案。

NumPy随机数

随机数常用于蒙特卡罗法、随机积分等方面。然而，真正的随机数很难获得，实际中使用的都是伪随机数。大部分情况下，伪随机数就足以满足我们的需求。当然，某些特殊情况除外，如进行高精度的模拟实验时。对于NumPy，与随机数有关的函数都在random子程序包中。

提示：

NumPy核心的随机数发生器是基于梅森旋转算法的，详见
https://en.wikipedia.org/wiki/Mersenne_twister$2。

我们既可以生成连续分布的随机数，也可以生成非连续分布的随机数。分布函数有一个可选的size参数，它能通知NumPy要创建多少个数字。我们可以用整型或者元组来给这个参数赋值，这时会得到相应形状的数组，其值由随机数填充。离散分布包括几何分布、超几何分布和二项式分布。连续分布包括正态分布和对数正态分布。

用二项式分布进行博弈
二项式分布模拟的是在进行整数次独立实验中成功的次数，其中每次实验的成功机会是一定的。

假如我们身在17世纪的赌场，正在对8片币玩法下注。当时流行用9枚硬币来玩。如果人头朝上的硬币少于5枚，那么我们将输掉一个8分币；否则，我们就赢一个8分币。下面，我们开始模拟这个游戏，假设我们手上有一千硬币作为初始赌资。我们将使用random模块提供的binomial()函数进行模拟：

提示：

完整的代码详见本书代码包中的headortail.py文件。

import numpy as np
from matplotlib.pyplot import plot, show

cash = np.zeros(10000)
cash[0] = 1000
outcome = np.random.binomial(9, 0.5, size=len(cash))

for i in range(1, len(cash)):

　　 if outcome[i] < 5:
　　　　　　 cash[i] = cash[i - 1] - 1
　　 elif outcome[i] < 10:
　　　　　　 cash[i] = cash[i - 1] + 1
　　 else:
　　　　　　 raise AssertionError("Unexpected outcome " + outcome)

print outcome.min(), outcome.max()

plot(np.arange(len(cash)), cash)
show()

若要了解binomial()函数，请看下列步骤。

1．用binomial()函数。

将数组初始化为零，即现金余额为零。调用binomial()函数，将size参数设为10000，这就表示我们要掷10000次硬币。

cash = np.zeros(10000)
cash[0] = 1000
outcome = np.random.binomial(9, 0.5, size=len(cash))

2．更新现金余额。

利用抛郑硬币的结果来更新cash数组。显示outcome数组中的最大值和最小值，只要保证没有罕见的异常值即可。

for i in range(1, len(cash)):
　　 if outcome[i] < 5:
　　　　 cash[i] = cash[i - 1] - 1
　　 elif outcome[i] < 10:
　　　　　　cash[i] = cash[i - 1] + 1
　　 else:
　　　　　　raise AssertionError("Unexpected outcome " + outcome)
print outcome.min(), outcome.max()

不出所料，其值在0～9。

0 9

3．用matplotlib绘出cash数组的图像。

plot(np.arange(len(cash)), cash)
show()

从图3-1可以看出，现金余额的曲线类似于随机行走（即不按照固定模式，而是随机游走）。

图3-1

当然，我们每一次执行程序代码，都会得到一个不同的随机游走。若想总是得到相同的结果，需要给NumPy的随机数子程序包中的binomial()函数一个种子值。

正态分布采样
连续分布是通过 概率密度函数（pdf） 进行建模的。在特定区间发生某事件的可能性可以通过概率密度函数的积分运算求出。NumPy的random模块提供了许多表示连续分布的函数，如beta、chisquare、exponential、f、gamma、gumbel、laplace、
lognormal、logistic、multivariate_normal、noncentral_chisquare、
noncentral_f、normal等。

利用NumPy的random子程序包中的normal()函数，可以把正态分布以直观的形式图示出来。下面给出随机产生的数值的钟形曲线和条形图（完整的脚本详见本书代码包中的normaldist.py文件）：

import numpy as np
import matplotlib.pyplot as plt

N=10000

normal_values = np.random.normal(size=N)
dummy, bins, dummy = plt.hist(normal_values, np.sqrt(N), normed=True, 
lw=1)
sigma = 1
mu = 0
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins - 
mu)**2 / (2 * sigma**2) ),lw=2)
plt.show()

随机数可以按照正态分布的要求生成，它们的分布情况可以用条形图来显示。若要绘制正态分布，请执行下列步骤。

1．成数值。

借助于NumPy的random子程序包提供的normal()函数，可以创建指定数量的随机数。

N=100.00
normal_values = np.random.normal(size=N)

2．画出条形图和理论上的pdf。

下面使用matplotlib来绘制条形图和理论上的pdf，这里中心值为0，标准差为1：

dummy, bins, dummy = plt.hist(normal_values,　　
　　np.sqrt(N), normed=True, lw=1)
sigma = 1
mu = 0
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi))　　
　　* np.exp( - (bins - mu)**2 / (2 * sigma**2) ),lw=2)
plt.show()

图3-2展示了著名的钟形曲线。

用SciPy进行正态检验

正态分布被广泛用于科学和统计学领域，按照中心极限定理，随着独立观测的随机样本数量的增加，它们会趋向呈正态分布。正态分布的特性已经为大家所熟知，并且它还非常便于使用，不过，它需要满足许多必要条件，如数据点的数量要足够大，并且还要求这些数据点必须是相互独立的。工作中，我们需要养成检查数据是否符合正态分布的好习惯。正态检验的方法有很多，其中有一些已经在scipy.stats程序包中实现了。本节将用到这些检验方法，这里的样本是取自
https://www.google.org/flutrends/data.txt$2的流感趋势数据。这里我们对原始文件进行了简化处理，仅留下日期和来自阿根廷的数据两列，下面给出示例行：

Date,Argentina
29/12/02,
05/01/03,
12/01/03,
19/01/03,
26/01/03,
02/02/03,136

这些数据还可以从本书代码包中的goog_flutrends.csv文件中找到，和前面一样，这里也按照正态分布对这些数据进行抽样。得到的数组的大小与流感趋势数组相同，如果顺利通过正态检验，就以它为金标准（the golden standard）。

提示：

这里的代码取自代码包中的normality_test.py文件。

import numpy as np
from scipy.stats import shapiro
from scipy.stats import anderson
from scipy.stats import normaltest

flutrends = np.loadtxt("goog_flutrends.csv", delimiter=',', 
usecols=(1,), skiprows=1, converters = {1: lambda s: float(s or 0)}, 
unpack=True)
N = len(flutrends)
normal_values = np.random.normal(size=N)
zero_values = np.zeros(N)

print "Normal Values Shapiro", shapiro(normal_values)
print "Zeroes Shapiro", shapiro(zero_values)
print "Flu Shapiro", shapiro(flutrends)
print

print "Normal Values Anderson", anderson(normal_values)
print "Zeroes Anderson", anderson(zero_values)
print "Flu Anderson", anderson(flutrends)
print

print "Normal Values normaltest", normaltest(normal_values)
print "Zeroes normaltest", normaltest(zero_values)
print "Flu normaltest", normaltest(flutrends)

作为一个反面例子，下面使用的这个数组跟前面提到的填满零的数组具有同样的大小。实际上，当处理一些小概率事件（如世界性疫情的爆发）时，就可能遇到这种数据。

在这个数据文件中，一些单元是空的；当然，这种问题经常会碰到，所以必须习惯于首先清洗数据。我们假定正确的数值是0，下面使用一个转换器来填上这些0值：

flutrends = np.loadtxt("goog_flutrends.csv", delimiter=',',　　
usecols=(1,), skiprows=1, converters = {1: lambda s:　　
float(s or 0)}, unpack=True)

夏皮罗-威尔克检验法可以对正态性进行检验。相应的SciPy函数将返回一个元组，其中第一个元素是一个检验统计量，第二个数值是 p 值。需要注意的是，这种填满了零的数组会引发警告。事实上，本例中用到的3个函数在处理这个数组方面确实有困难，同时还引发警告。下面是得到的结果：

Normal Values Shapiro (0.9967482686042786,　　
0.2774980068206787)
Zeroes Shapiro (1.0, 1.0)
Flu Shapiro (0.9351990818977356, 2.2945883254311397e-15)

这种用零填充的数组看起来有点怪，虽然会收到相关的警告，但是大可以忽略不计。这里得到的p值类似于本例中后面的第三个检验的结果。分析结果基本上是一样的。

Anderson-Darling检验可以用来检验正态分布以及其他分布，如指数分布、对数分布和冈贝尔（Gumbel）分布等。相关的SciPy函数涉及一个检验统计量和一个数组，该数组存放了百分之15、10、5、2.5和1等显著性水平所对应的临界值。如果该统计量大于显著性水平的临界值，就可以断定它不具有正态性。我们将得到下列数值：

Normal Values Anderson (0.31201465602225653, array([ 0.572,　　 
0.652,　　0.782,　　0.912,　　1.085]), array([ 15. ,　　10. ,　　 5.　　
,　　 2.5,　　 1. ]))
Zeroes Anderson (nan, array([ 0.572,　　0.652,　　0.782,　　 
0.912,　　1.085]), array([ 15. ,　　10. ,　　 5. ,　　 2.5,　　 1.　　
]))
Flu Anderson (8.258614154768793, array([ 0.572,　　0.652,　　 
0.782,　　0.912,　　1.085]), array([ 15. ,　　10. ,　　 5. ,　　 2.5,　　　　
1. ]))

对于这种用零填充的数组，我们没有什么好说的，因为返回的统计量根本就不是一个数值。就像我们期望的那样，这个金标准数组是符合正态分布的。然而，流感趋势数据返回的这个统计量大于所有的有关临界值，因此我们可以肯定它不符合正态分布。这3个测试函数中，这个看起来是最容易使用的一个。

SciPy的normaltest()函数还实现了D’Agostino检验和Pearson检验功能。该函数返回的元组和shapiro()函数一样，也包括一个统计量和p值。这里的p值是双边卡方概率（two-sided Chi-squared probability），卡方分布是另一种著名的分布，这种检验本身是基于偏度和峰度检验的z分数的。偏度系数用来表示分布的对称程度的，这样，由于正态分布是对称的，所以偏度系数为零。峰度系数描述的是分布的形状（尖峰，肥尾）。这种正态分布的峰度系数为3，超额峰度的系数为0。以下是本次检验获得的数值：

Normal Values normaltest (3.102791866779639, 0.21195189649335339)
Zeroes normaltest (1.0095473240349975, 0.60364218712103535)
Flu normaltest (99.643733363569538, 2.3048264115368721e-22)

因为这里处理的是p值的概率，所以它越大越好，最好接近1。对于这种用零填充的数组，会得到很奇怪的结果，不过，由于我们已经得到提示，所以对于这个特殊数组来说，结果是不可信的。此外，如果p值大于等于0.5，则认为它具有正态性。对于金标准数组，我们会得到一个更小的数值，这说明我们需要进行更多的观察，这一点留作作业，请读者自行练习。