numpy下 随机抽样

numpy.random模块可以用来生成多种概率分布的样本值的函数，是对python内置的random函数的补充。

numpy.ranom.seed(seed=None) 生成随机数

使用相同的seed值，可以保证每次都生成相同的随机数。

离散型随机变量

二项分布

numpy.random.binomial(n, p, size=None)

表示对一个二项分布进行采样，n为n重伯努利实验次数，p为成功的概率，size为采样的次数，函数返回n中成功的次数。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy import stats

# 野外正在进行9（n=9）口石油勘探井的发掘工作，每一口井能够开发出油的概率是0.1（p=0.1）。
# 请问，最终所有的勘探井都勘探失败的概率？

# 生成一个随机种子，固定时间使每次生成的随机数相同
np.random.seed(20200605)
n = 9
p = 0.1
size = 50000
# n：n重伯努利试验，p：成功的概率，size：采样的次数（每次采样做n次试验）
# 函数返回n次中成功的次数
x = np.random.binomial(n, p, size)

'''
或者使用binom.rvs(n, p, size=1)函数模拟一个二项随机变量,可视化地表现概率
y = stats.binom.rvs(n, p, size=size)#返回一个numpy.ndarray
'''
# mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体：解决plot不能显示中文问题
# mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题

print(np.sum(x == 0) / size) # n次实验均失败的概率

plt.hist(x) # 绘制直方图（histogram）
plt.xlabel('随机变量：成功次数')
plt.ylabel('样本中出现的次数')
plt.show()

# 它返回一个列表，列表中每个元素表示随机变量中对应值（n次中成功k次）的概率
s = stats.binom.pmf(range(10), n, p)
print(np.around(s, 3))
# [0.387 0.387 0.172 0.045 0.007 0.001 0. 0. 0. 0. ]

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 模拟投硬币，投2次，两次都为正面的概率？

np.random.seed(20200605)
n = 2
p = 0.5 # 投硬币为正面的概率
size = 50000
x = np.random.binomial(n, p, size)

print(np.sum(x == 0) / size) # 0.25154
print(np.sum(x == 1) / size) # 0.49874
print(np.sum(x == 2) / size) # 0.24972

# density为true，则返回的元组的第一个参数frequency将为频率而非默认的频数。
plt.hist(x, density=True) 
plt.xlabel('随机变量：硬币为正面次数')
plt.ylabel('50000个样本中出现的次数')
plt.show()

# 它返回一个列表，列表中每个元素表示随机变量中对应值的概率
s = stats.binom.pmf(range(n + 1), n, p)
print(np.around(s, 3))
# [0.25 0.5 0.25]

# 计算期望和方差 m为期望 n为方差
# 期望：E(x) = np
# 方差：Var(x) = np(1‐p)
m,v = stats.binom.stats(n, p, loc=0, moments='mv')
print(m,v)

泊松分布

numpy.random.poisson(lam=1.0, size=None)

表示对一个泊松分布进行采样，lam表示单位时间内发生事件的平均值，函数返回单位时间内事件发生次数。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 假定某航空公司预定票处平均每小时接到42次订票电话，那么10分钟内恰好接到6次电话的概率是多少

np.random.seed(20200605)
lam = 42 / 6 # 平均值：平均每十分钟接到42/6次订票电话
size = 50000
# lam 表示一个单位内发生事件的平均值，函数的返回值表示一个单位内事件发生的次数
x = np.random.poisson(lam, size)
'''或者
# 模拟服从泊松分布的50000个随机变量
x = stats.poisson.rvs(lam,size=size)
'''
print(np.sum(x == 6) / size) # 0.14988

plt.hist(x)
plt.xlabel('随机变量：每十分钟接到订票电话的次数')
plt.ylabel('50000个样本中出现的次数')
plt.show()

# 用poisson.pmf(k, mu)求对应分布的概率:概率质量函数 (PMF)
x = stats.poisson.pmf(6, lam)
print(x) # 0.14900277967433773

超几何分布

numpy.random.hypergeometric(ngood, nbad, nsample, size=None)

表示对一个超几何分布进行采样，ngood表示成功标志元素个数，nbad表示没有成功标志元素个数，nsample表示抽样次数，函数返回抽取nsample个元素中具有成功标志的元素个数。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 一共20只动物里有7只是狗，抽取12只有3只狗的概率（无放回抽样）。

np.random.seed(20200605)
size = 500000
# ngood 表示总体中具有成功标志的元素个数，nbad 表示总体中不具有成功标志的元素个数，ngood+nbad 表示总体样本容量
# nsample 表示抽取元素的次数（小于或等于总体样本容量），函数的返回值表示抽取nsample 个元素中具有成功标识的元素个数。
x = np.random.hypergeometric(ngood=7, nbad=13, nsample=12, size=size)
'''或者
# 用rvs(M, n, N, loc=0, size=1, random_state=None)模拟
# M 为总体容量，n 为总体中具有成功标志的元素的个数，N、k 表示抽取N个元素有k个是成功元素
x = stats.hypergeom.rvs(M=20,n=7,N=12,size=size)
'''
print(np.sum(x == 3) / size) # 0.198664

plt.hist(x, bins=8) # bins为柱子个数
plt.xlabel('狗的数量')
plt.ylabel('50000个样本中出现的次数')
plt.title('超几何分布',fontsize=20)
plt.show()

x = range(8)
#用hypergeom.pmf(k, M, n, N, loc)来计算k次成功的概率
s = stats.hypergeom.pmf(k=x, M=20, n=7, N=12)
print(np.round(s, 3))
# [0. 0.004 0.048 0.199 0.358 0.286 0.095 0.01 ]

# 用stats(M, n, N, loc=0, moments='mv')计算均值和方差
# 均值E(x) = N(n/M)
# 方差Var(x) = N(n/M)(1‐n/M)((M‐N)/(M‐1))
m,v = stats.hypergeom.stats(20,7,12,moments='mv')
print(m,v)

连续型随机变量

均匀分布

numpy.random.uniform(low=0.0, high=1.0, size=None)

表示在均匀分布[low,high)中抽取样本

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 在low到high范围内，创建大小为size的均匀分布的随机数。

np.random.seed(20200614)
a = 0 # low
b = 100 # high
size = 50000
x = np.random.uniform(a, b, size=size) # 返回值：ndarray类型
print(x)

print(np.all(x >= 0)) # True
print(np.all(x < 100)) # True
y = (np.sum(x < 50) - np.sum(x < 10)) / size
print(y) # 0.40144

plt.hist(x, bins=20)
plt.title('均匀分布',fontsize=20)
plt.show()

a = stats.uniform.cdf(10, 0, 100) # 累计概率密度
b = stats.uniform.cdf(50, 0, 100)
print(b - a) # 0.4

numpy.random.rand(d0, d1, …, dn)

可以得到均匀分布在[0,1)之间的随机数

import numpy as np

# 根据指定大小产生[0,1)之间均匀分布的随机数。

np.random.seed(20200614)
print(np.random.rand())
# 0.7594819171852776

print(np.random.rand(5))
# [0.75165827 0.16552651 0.0538581 0.46671446 0.89076925]

print(np.random.rand(4, 3))
# [[0.10073292 0.14624784 0.40273923]
# [0.21844459 0.22226682 0.37246217]
# [0.50334257 0.01714939 0.47780388]
# [0.08755349 0.86500477 0.70566398]]

np.random.seed(20200614)
print(np.random.uniform()) # 0.7594819171852776
print(np.random.uniform(size=5))
# [0.75165827 0.16552651 0.0538581 0.46671446 0.89076925]

print(np.random.uniform(size=(4, 3)))
# [[0.10073292 0.14624784 0.40273923]
# [0.21844459 0.22226682 0.37246217]
# [0.50334257 0.01714939 0.47780388]
# [0.08755349 0.86500477 0.70566398]]

numpy.random.randint(low, high=None, size=None, dtype=‘l’)

可以得到均匀分布在[low,high)之间的随机整数

import numpy as np

# 若high 不为None 时，取[low,high)之间随机整数，否则取值[0,low)之间随机整数。

np.random.seed(20200614)
x = np.random.randint(2, size=10) # 取[0,2)之间整数
print(x)
# [0 0 0 1 0 1 0 0 0 0]

x = np.random.randint(1, size=10) # [0,1)
print(x)
# [0 0 0 0 0 0 0 0 0 0]

x = np.random.randint(5, size=(2, 4)) # [0,5)
print(x)
# [[3 3 0 1]
# [1 1 0 1]]

x = np.random.randint(1, 10, [3, 4]) # [1,10)
print(x)
# [[2 1 7 7]
# [7 2 4 6]
# [8 7 2 8]]

正态分布

numpy.random.randn(d0, d1, …, dn)

得到标准正态分布的样本

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 根据指定大小产生满足标准正态分布的数组（均值为0，标准差为1）。

np.random.seed(20200614)
size = 50000
x = np.random.randn(size)
print(x)
y1 = (np.sum(x < 1) - np.sum(x < -1)) / size
y2 = (np.sum(x < 2) - np.sum(x < -2)) / size
y3 = (np.sum(x < 3) - np.sum(x < -3)) / size
print(y1) # 0.68596
print(y2) # 0.95456
print(y3) # 0.99744

plt.hist(x, bins=20)
plt.title('正态分布',fontsize=20)
plt.show()

y1 = stats.norm.cdf(1) - stats.norm.cdf(-1) # cdf 累计概率密度
y2 = stats.norm.cdf(2) - stats.norm.cdf(-2)
y3 = stats.norm.cdf(3) - stats.norm.cdf(-3)
print(y1) # 0.6826894921370859
print(y2) # 0.9544997361036416
print(y3) # 0.9973002039367398

numpy.random.normal(loc=0.0, scale=1.0, size=None)

得到给定均值和标准差的高斯分布（正态分布）样本，其中loc为均值，scale为标准差。

# numpy.random.normal(loc=0.0, scale=1.0, size=None)
import numpy as np
import matplotlib.pyplot as plt

# 生成均值为5，标准差为0.5的高斯分布（正态分布）

np.random.seed(20200614)
x = 0.5 * np.random.randn(2, 4) + 5 # 2×4维，使用标准正态分布生成高斯分布（×sigma+mu）
'''或者
# 模拟10000个随机变量
x = 0.5*stats.norm.rvs(size=(2,4))+5
'''
print(x)
# [[5.39654234 5.4088702 5.49104652 4.95817289]
# [4.31977933 4.76502391 4.70720327 4.36239023]]

np.random.seed(20200614)
mu = 5 # 平均值
sigma = 0.5 # 标准差
# normal()为创建均值为 loc（mu），标准差为 scale（sigma），大小为 size 的数组。
x = np.random.normal(mu, sigma, (2, 4))
print(x)
# [[5.39654234 5.4088702 5.49104652 4.95817289]
# [4.31977933 4.76502391 4.70720327 4.36239023]]

size = 50000
x = np.random.normal(mu, sigma, size)
print(np.mean(x)) # 4.996403463175092
print(np.std(x, ddof=1)) # 0.4986846716715106（样本标准差） ddof=1使得正确调用std


plt.hist(x, bins=20)
plt.title('高斯分布',fontsize=20)
plt.show()

指数分布

numpy.random.exponential(scale=1.0, size=None)

得到指数分布的样本

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 指数分布 scale = 1/lambda 标准差

np.random.seed(20200614)
lam = 7
size = 50000
x = np.random.exponential(1 / lam, size)
# 或者rvs(loc=0, scale=1/lam, size=size, random_state=None)

y1 = (np.sum(x < 1 / 7)) / size
y2 = (np.sum(x < 2 / 7)) / size
y3 = (np.sum(x < 3 / 7)) / size
print(y1) # 0.63218
print(y2) # 0.86518
print(y3) # 0.95056

plt.hist(x, bins=20)
plt.title('指数分布',fontsize=20)
plt.show()

y1 = stats.expon.cdf(1 / 7, scale=1 / lam) # F(1/7) F(x)=1-e^(-lam*x)
y2 = stats.expon.cdf(2 / 7, scale=1 / lam)
y3 = stats.expon.cdf(3 / 7, scale=1 / lam)
print(y1) # 0.6321205588285577
print(y2) # 0.8646647167633873
print(y3) # 0.950212931632136

其他随机函数

随机从序列中获取元素

numpy.random.choice(a, size=None, replace=True, p=None)

若a为一整数，则从range(a)中获取元素；若a为数组，则从数组元素中获取。replace可指定获取的元素是否重复。p可指定选取a中各元素的概率。

import numpy as np

# 随机从序列中获取元素

np.random.seed(20200614)
x = np.random.choice(10, 3) # 从range(10)中随机获取
print(x) # [2 0 1]

x = np.random.choice(10, 3, p=[0.05, 0, 0.05, 0.9, 0, 0, 0, 0, 0, 0]) # 以一定概率从range(10)中随机获取
print(x) # [3 2 3]

x = np.random.choice(10, 3, replace=False, p=[0.05, 0, 0.05, 0.9, 0, 0, 0, 0, 0, 0]) # 以一定概率不重复从range(10)中随机获取
print(x) # [3 0 2]

aa_milne_arr = ['pooh', 'rabbit', 'piglet', 'Christopher']
x = np.random.choice(aa_milne_arr, 5, p=[0.5, 0.1, 0.1, 0.3]) # 从数组中获取
print(x) # ['pooh' 'rabbit' 'pooh' 'pooh' 'pooh']

np.random.seed(20200614)
x = np.random.randint(0, 10, 3) # 返回0到10之间的任意整数
print(x) # [2 0 1]

对数据集进行洗牌操作

机器学习中很多算法要求数据之间相互独立，所以要先对数据集进行洗牌操作（改变元素序列）。

numpy.random.shuffle(x)

shuffle函数在原数组上进行操作，改变自身序列，无返回值。

import numpy as np

# 洗牌，改变自身内容，打乱顺序。

np.random.seed(20200614)
x = np.arange(10)
print(x)
np.random.shuffle(x) # 洗牌
print(x)
# [6 8 7 5 3 9 1 4 0 2]

print(np.random.shuffle([1, 4, 9, 12, 15])) # 无返回值
# None

x = np.arange(20).reshape((5, 4))
print(x)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]
# [12 13 14 15]
# [16 17 18 19]]

np.random.shuffle(x) # 在原数组上进行
print(x)
# [[ 8  9 10 11]
#  [ 0  1  2  3]
#  [12 13 14 15]
#  [16 17 18 19]
#  [ 4  5  6  7]]

numpy.random.permutation(x)

permutation函数不在原数组上进行操作，不改变自身数组，返回新数组。

import numpy as np

# 洗牌

np.random.seed(20200614)
x = np.arange(10)
y = np.random.permutation(x) 
print(y)
# [6 8 7 5 3 9 1 4 0 2]

print(np.random.permutation([1, 4, 9, 12, 15])) # 返回新的数组
# [ 4 1 9 15 12]

x = np.arange(20).reshape((5, 4))
print(x)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]
# [12 13 14 15]
# [16 17 18 19]]

y = np.random.permutation(x) # 在新的数组的上进行
print(y)
# [[ 8 9 10 11]
# [ 0 1 2 3]
# [12 13 14 15]
# [16 17 18 19]
# [ 4 5 6 7]]
print(x)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]
# [12 13 14 15]
# [16 17 18 19]]