python随机抽取数据 python 随机取样

一、什么是采样

• 在信号系统和数字信号处理中，采样是每隔一定的时间测量一次声音信号的幅值，把时间连续的模拟信号转换成时间离散幅值的采样信号。
• 如果采样的时间间隔相等，这种采样称为均匀采样。
• 在计算机系统中有一个重要的问题叫随机采样：就是给定一个概率分布$p(x)$,我们如何按照该分布产生样本
• 在机器学习或统计学习中,我们是给定一堆样本数据,通过参数估计的方法求出这堆样本所符合的概率分布$p(x)$;而随机采样刚相反：给定一个概率分布p(x)，如何生成满足条件的样本？

二、均匀分布的采样

• 首先需要明确,计算机程序都是确定性的,因此并不能产生真正意义上的随机数,只能产生伪随机数.
• 另外由于计算机的存储和计算单元只能处理离散状态值,因此也不能产生连续均匀分布的随机数,只能通过离散分布来逼近连续分布.
• 均匀分布式是一种最简单的分布。在计算机中使rand()函数生成的伪随机数序列，就可以看做是一种均匀分布。
• 一般可采用线性同余法来生成离散均匀分布伪随机数,计算公式为:$x_{t+1}=(ax_t+c) mod m$
• 根据当前生成的随机数$x_t$来进行适当变换,进而产生下一次的随机数$x_{t+1}$.初始值$x_0$称为随机种子.上式可得到区间$[0,m-1]$上的伪随机数.
• 想要真正的随机数可以通过自然界的物理产生,比如放射性物质的衰变,温度,气流的扰动等.有一些网站可以提供基于大自然的随机现象的随机数,有兴趣的读者可以尝试一下.
• 代码：

import time
start = int(time.time()); # 以当前时间作为种子
m = (1 << 31) -1; # 2^31 -1
a = 1103515245;
b = 12345;
R=start
for i in range(5):
    R = ( a *R + b ) % m;
    print(R/m)#产生[0,1]的随机数

三、离散分布采样

• 离散分布是比均匀分布稍微复杂一点的情况。下面用一个例子来详细说明：
• 令概率分布为：$p(x)=[0.1,0.2,0.3,0.4]$，$x=["hello", “java”, “python”, "scala"]$,我们可以使用将每个概率视为区间[0,1]中的一个小线段，线段长度为概率的大小，然后用上面均匀分布产生随机数，随机所在区间就是当前要采的样本。
• 代码：

import numpy as np
from collections import defaultdict

dic = defaultdict(int)#key是x的取值，value是x被采样的次数
def sample():
    u = np.random.rand()#在[0,1]上按照均匀分布产生一个随机数
    if u <= 0.1:#[0,0.1]=>x=hello
        dic["hello"] += 1
    elif u <= 0.3:#(0.1,0.3]=>x=java
        dic["java"] += 1
    elif u <= 0.6:#(0.3,0.6]=>x=python
        dic["python"] += 1
    else:#(0.6,1]=>x=scala
        dic["scala"] += 1


for i in range(10000):#进行100词采样
    sample()
for k,v in dic.items():
    print(k,v)

输出：
python 2991
scala 3940
hello 1037
java 2032

四、逆变换采样法

• 对于连续型随机变量$X$，一点有两个函数：

• 概率密度函数：$f(x)$
• 分布函数：$F_X (x)=P(X≤x)=∫_{-∞}^xf(t)dt$ :表示随后变量取值在区间$(-∞,x]$的概率
• 两个函数的图如下所示：



• 
  
• 并且有下面性质：

• $F_X (x)$是$f(x)$的一个原函数
• 并且$0≤F_X (x)≤1$.
• $F_X (x)$是一个单独递增的函数，并存在反函数$F_X^{-1}$.

• 重要定理： 若联系型随机变量$Y$与连续型随机变量$X$是关系为：$Y=F_X (X)$,即他们的关系函数是$X$是分布函数。则随机变量$Y$是$0,1]$上的均匀分布,即：$Y$~$U(0.1)$.下面证明一下：$\because Y=F_X (X)~~~ \therefore X=F_X^{-1} (Y)$ $\because F_X (x)是单调递增的~~~ \therefore F_X^{-1} (Y)也是单调递增的$ $\because F_X (x)=P(X≤x)=P(F_X^{-1} (Y)≤x)=P(Y≤y$$\therefore P(Y≤F_X (x))=F_X (x),且F_X (x)∈[0,1],$ 这恰恰是$[0,1]$上均匀分布的分布函数,所以：$Y$~$U(0,1)$
• 若Y~U(0,1)则(只是为了说明上面的证明)：

• 概率密度函数为：$f(y) = \left\{ \begin{array}{lr} 1 & 0≤y≤1\\ 0 & 其他 \end{array} \right.$
• 分布函数：：$F_Y(y) = \left\{ \begin{array}{lr} 0 & 0<y\\ y & 0≤y≤1\\ 1 & y>1 \end{array} \right.$

• 有了上面这个定理，我们可以：(逆变换采样法)

• 首先对于随机变量Y进行均匀采样抽取一个数$u$
• 由于：$X=F_X^{-1} (Y)$
• 我们可以得到一个随机变量X的样本$x$：$x=F_X^{-1} (u)$
• 过程可以如下图所示：



• 
  

• 例子：

• 假设随机变量X的概率密度函数为：$f(x)=3x^2,(0<x<1)$
• 则其相应的分布函数为：$F_X (x)=x^3, (0<x<1)$
• 在区间[0,1]上进行均匀分布的采样,得到序列:$[u_1,u_2,…,u_n ]$
• 将序列带入$F_X^{-1} (u)=u^{1/3}$,得到随机变量X的样本$[u_1^{1/3},u_2^{1/3},…,u_n^{1/3}]=[x_1,x_2,…,x_n ]$

五、Box-Muller算法

• Box-Muller算法实质是逆变换采样法的一个特例，即随机变量$X$服从正态分布的情况。
• 假设$X$~$N(0,1)$,其概率密度函数为：$f(x)=\frac{1}{\sqrt{2π}} e^{-\frac{x^2}{2}}$
• 但是注意到直接对$f(x)$进行定积分得到其分布函数$F_X (x)$是行不通。
• 假设：$X$ ~ $N(0,1)$和$Y$~$N(0,1)$并且$X$和$Y$相互独立，则其联合概率密度函数为：$f(x,y)=\frac{1}{2π} e^{-\frac{x^2+y^2}{2}}$
• 进行下面随机变量变换：$\left\{ \begin{array}{lr} X=Rcos(θ)\\ Y=Rsin(θ) \end{array} \right.$其中$R$和$θ$都为随机变量
• 则有：$∫_{-∞}^{+∞}∫_{-∞}^{+∞} \frac{1}{2π} e^{-\frac{x^2+y^2}{2}} dxdy=∫_{0}^{2π}∫_{0}^{+∞} \frac{1}{2π} e^{-\frac{r^2}{2}} rdrdθ=1$
• 则随机变量$R$和$θ$的边缘概率分布函数为：$F_R (r)=∫_{0}^{2π}∫_{0}^{r} \frac{1}{2π} e^{-\frac{r^2}{2}} rdrdθ=1-e^\frac{-r^2}{2}$ $F_θ (∅)=∫_{0}^{∅}∫_{0}^{+∞} \frac{1}{2π} e^{-\frac{r^2}{2}} rdrdθdxdy=\frac{∅}{2π}$
• 由$F_θ (∅)$表达式可知：$θ$~$U(0,2π)$
• 令随机变量$Z$与随机变量$R$的关系为：$Z=F_R (R)=1-e^\frac{-r^2}{2}$$R=F_R^{-1} (Z)=\sqrt{-2ln⁡(1-Z)}$$则Z$~$U(0,1)$
• 令$2πU_1=θ$~$U(0,2π)$ 和$1-U_2=Z$~$U(0,1)$ 则：$U_1$~$U(0,1); U_2$~$U(0,1)$
• 则$U_1$和$U_2$可以通过按$[0,1]$上的均匀分布采样得到。并且有：

• $θ=2πU_1$
• $Z=1-U_2$
• $R=\sqrt{-2ln⁡(1-Z)}=\sqrt{-2n⁡U_2 }$
• $X=Rcos(θ)$
• $Y=Rsin(θ)$

• 代码：

import numpy as np
from pylab import *

def sample():
    U1 = np.random.rand()  # 两个均匀分布分别为x, y
    U2 = np.random.rand()
    theta = 2 * np.pi * U1
    R = np.sqrt(-2 * np.log(U2))
    X= R * np.cos(theta)
    Y= R * np.sin(theta)
    return X


def sampleNtimes():
    list = []
    n = 100000
    for i in range(n):
        x = sample()
        list.append(x)
    y=np.array(list).reshape(n,1)
    hist(y, 20,normed=1, fc='c')  # 采样直方图

    x = arange(-5, 5, 0.1)
    plot(x, 1 / np.sqrt(2 * np.pi) * np.exp(-0.5 * x ** 2), 'g', lw=6)  # 标准正态分布
    xlabel('x', fontsize = 24)
    ylabel('p(x)', fontsize = 24)
    show()
sampleNtimes()

输出：

五、接受/拒绝采样

• 很多实际问题中，对于目标分布$p(z)$直接采样比较难, 因此选取一个容易采样的参考分布$q(z)$比如高斯分布,使得对于任意$z$都有$p(z)≤kq(z)$,则可以按如下过程采用:

• 从参考分布$q(z)$中随机抽取一个样本$z_i$
• 从均匀分布$U(0,1$)产生一个随机数$u_i$
• 如果$u_i≤\frac{P(z_i )}{kq(z_i}$ ,则接受样本$z_i$;否则拒绝重新进行步骤1~3,直到新产生的样本 被接
• $p(z)$和$q(z)$的关系图：



• 
  

• 现在我们来分析一下接受/拒绝采样在$z_i$采样数量的决定因素：

• 第一个因素是$P(z_i )$ ：因为如果$P(z_i )$大，那么$q(z_i )$就会比较大，而我们采样是按照$q(z_i )$来进行采样的则$z_i$被采样的概率就大，
• 第一个因素是$\frac{P(z_i )}{kq(z_i )}$ ：$0≤\frac{P(z_i )}{kq(z_i )} ≤1$ 越大而$u_i$是$[0,1]$上的均匀分布，则$\frac{P(z_i )}{kq(z_i )} >u_i$
• 当然最主要的因素是第一个因素。

• 拒绝采用的关键是为目标分布$p(z)$选取一个合适的包络网络函数$kq(z)$:包络函数越紧,每次采样时样本被接受的概率越大,采样效率越高.
• 代码：

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from pylab import *

#标准正态分布
def qsample():
    U1 = np.random.rand()  # 两个均匀分布分别为x, y
    U2 = np.random.rand()
    theta = 2 * np.pi * U1
    R = np.sqrt(-2 * np.log(U2))
    X= R * np.cos(theta)
    Y= R * np.sin(theta)
    return X

#正态分布函数
def q(x):
    return 1 / np.sqrt(2 * np.pi) * np.exp(-0.5 * x ** 2)

#目标分布
def p(x):
    x=x+1.2#左移一点点让正态分布正好覆盖
    return 0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)


def rejection(nsamples):
    K=3  #  k值
    samples = np.zeros(nsamples, dtype=float)
    for i in range(nsamples):
        accept = False
        while not accept:
            x = qsample()
            u = np.random.rand()
            if u < p(x)/(K*q(x)):
                accept = True
                samples[i] = x
    return samples


x = np.arange(-4, 4, 0.01)
x2 = np.arange(-4.5, 4.5, 0.1)
realdata = p(x)  #目标概率函数p(x)图
stardline=3*q(x)    #正态分布函数q(x)
box = np.ones(len(x2))*0.8
plt.plot(x, realdata, 'g', lw=3)
plt.plot(x, stardline, 'r', lw=3)
samples = rejection(100000)

plt.hist(samples, 100, normed=1, fc='c')
plt.xlabel('x', fontsize=24)
plt.ylabel('p(x)', fontsize=24)
plt.show()

输出：

参考连接：

• https://drivingc.com/p/5c46be544b0f2b64f2179aea
• https://zhuanlan.zhihu.com/p/55653004
• https://zhuanlan.zhihu.com/p/38638710