随机采样和随机模拟

Posted 2022-08-13 -柚子皮-

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了随机采样和随机模拟相关的知识，希望对你有一定的参考价值。

http://blog.csdn.net/pipisorry/article/details/50615652

随机采样方法

模拟方法:是一种基于“随机数”的计算方法，基于数值采样的近似推断方法,也被称为蒙特卡罗( MonteCarlo )方法、随机模拟方法。

通常均匀分布Uniform(0,1) 的样本，即我们熟悉的类rand()函数，可以由线性同余发生器生成；而其他的随机分布都可以在均匀分布的基础上，通过某种函数变换得到，比如，正态分布可以通过Box-Muller变换得到。然而，这种变换依赖于计算目标分布的积分的反函数，当目标分布的形式很复杂，或者是高维分布时，很难简单变换得到。

当一个问题无法用分析的方法来求精确解，此时通常只能去推断该问题的近似解，而随机模拟(MCMC)就是求解近似解的一种强有力的方法。

随机模拟的核心就是对一个分布进行抽样(Sampling)。

Monte Carlo 方法有这些

产生独立样本
基本方法:概率积分变换(第一部分已讲)
接受—拒绝(舍选)采样
重要性采样
产生相关样本:Markov Chain Monte Carlo
Metropolis-Hastings算法
Gibbs Sampler

采样的目的：为什么要采样

对于大多数实际应用中的概率模型来说,精确推断是不可行的,因此我们不得不借助与某种形式的近似。有基于确定性近似的推断方法,它包括诸如变分贝叶斯方法以及期望传播。这里,我们考虑蒙特卡罗方法。
虽然对于一些应用来说,我们感兴趣的是非观测变量上的后验概率分布本身,但是在大部分情况下,后验概率分布的主要用途是计算期望,例如在做预测的情形下就是这样。因此,我们希望解决的基本的问题涉及到关于一个概率分布 p(z) 寻找某个函数 f (z) 的期望。这里, z 的元素可能是离散变量、连续变量或者二者的组合。因此,在连续变量的情形下,我们希望计算下面的期望
∫E[f ] =f (z)p(z) d z
在离散变量的情形下,积分被替换为求和。图11.1图形化地说明了单一连续变量的情形。我们假设,使用解析的方法精确地求出这种期望是十分复杂的。

采样方法背后的一般思想是得到从概率分布 p(z) 中独立抽取的一组变量 z (l) , 其中 l = 1, . . . , L 。这使得期望可以通过有限和的方式计算。
只要样本 z (l) 是从概率分布 p(z) 中抽取的,那么 E[ f ] = E[f ] ,因此估计 f 具有正确的均值。值得强调的一点是,估计的精度不依赖于 z 的维度,并且原则上,对于数量相对较少的样本 z (l) ,可能会达到较高的精度。在实际应用中,10个或者20个独立的样本就能够以足够高的精度对期望做出估计。

然而,问题在于样本 z (l) 可能不是独立的,因此有效样本大小可能远远小于表面上的样本大小。并且,回到图11.1,我们注意到如果 f (z) 在 p(z) 较大的区域中的值较小,反之亦然,那么期望可能由小概率的区域控制,表明为了达到足够的精度,需要相对较大的样本大小。

在由无向图定义的概率分布的情形中,如果我们希望从没有观测变量的先验概率分布中采样,那么不存在一遍采样的方法。相反,我们必须使用计算量更大的方法,例如吉布斯采样。
除了从条件概率分布中采样之外,我们可能也需要从边缘概率分布中采样。如果我们已经有了一种从联合概率分布 p(x, v) 中采样的方法,那么得到从边缘概率分布 p(u) 中的样本是很容易的,只需忽略每个样本中的 v 的值即可。

应用采样进行积分近似的实例

所以对于蒙特卡罗积分来说，最重要的当然就是怎么采样了，采样出来了，使用均值进行近似就ok了。

皮皮blog

这里只写一下直接采样、接受-拒绝采样和吉布斯采样，其它具体的，太懒了，给个网上的链接，自己看看吧。。。

直接采样

通过对均匀分布采样，实现对任意分布采样。

一般而言均匀分布 Uniform(0,1)的样本是相对容易生成的。通过线性同余发生器可以生成伪随机数，我们用确定性算法生成[0,1]之间的伪随机数序列后，这些序列的各种统计指标和均匀分布 Uniform(0,1) 的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

而我们常见的概率分布，无论是连续的还是离散的帆布，都可以基于Uniform(0,1)的样本生成。

直接采样步骤：

从Uniform(0,1)随机产生一个样本z
另z=h(y)，其中h(y)为y的累积概率分布CDF
计算y=h−1(z)
结果y为对p(y)的采样

Note：需要知道累积概率分布的解析表达式，且累积概率分布函数存在反函数。但是如果h(x)不能确定或者没有无法解析求逆则直接抽样不再合适。对于复杂的现实模型其实是不常用的。

皮皮blog

接受-拒绝采样（Acceptance-Rejection Sampling）

简称拒绝采样，基本思想：假设我们需要对一个分布f(x)进行采样，但是却很难直接进行采样，所以我们想通过另外一个容易采样的分布g(x)的样本，用某种机制去除掉一些样本，从而使得剩下的样本就是来自与所求分布f(x)的样本。

采样过程

几何解释及数学证明

Note: 也就是随机产生了Mq(x)下的点（采样点x值的数目与Mq(x)的pdf成正相关），我们只接受pi(x)内部的点，这样采样点的就与pi(x)的pdf成正相关了。每个x点值对应的pi(x)/Mq(x)是一个[0,1]均匀分布的随机变量，几何上Mq(x)与pi(x)越接近，当然应该更容易接受（如图中接近a的第一个峰值对应的地方），而pi(x)/Mq(x)接近1，产生的均匀[0,1]变量U更容易<=pi(x)/Mq(x)，更容易接受。

接受率K(M)及其计算

一般写作K或者M

Note：可能是高维下更不容易找到接近的q(x)，这样两者之间空间大，容易被拒绝。

对截断正态分布采样实例1

Note: 这里接受率M的计算过程如下：

[概率论：正态分布]

对截断正态分布采样实例2

吉布斯采样Gibbs Sampling

[随机采样和随机模拟：吉布斯采样Gibbs Sampling]

其它随机采样方法参考

[Bishop [Pattern Recognition and Machine Learning. Springer, 2006.] Chapter 11 for discussion of rejection sampling, adaptive rejection sampling, adaptive rejection Metropolis sampling, importance sampling,sampling-importance-sampling,...]

[随机采样方法整理与讲解（MCMC、Gibbs Sampling等）]

[随机采样方法整理（MCMC、Gibbs Sampling等）]

[随机模拟的基本思想和常用采样方法（sampling）]

[Deep learning：十八(关于随机采样)]

正交采样

基于正交实验设计？比如有30个参数，每个参数有不同数目的可选值，随机采样出多组30个参数值作为一次实验，这些组正交？

[测试基础---测试用例之正交试验]

[发布了一个正交表计算的程序 (Python 调用的 pyd) c++ 写的]

正交示例 [正交表测试用例设计 ]

[http://blog.csdn.net/vincetest]

类似正交采样的拉丁超立方抽样latin hypercube sampling。

皮皮blog

随机采样算法的实现

接受-拒绝采样的python实现

对截断正态分布[Truncated normal distribution]采样

k的求解

不要求截断面积时的k值计算：

Note: 这里计算应该是错的，因为考虑的并不是截断正态分布，也就是正态分布函数相除时，少除了一个截断的面积，推导可参考上面的接受率M的计算过程。当然在python代码中不用这么复杂的计算了嘛，计算最小值计算机做就可以了。

接受-拒绝采样python代码

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
__title__ = '接受拒绝采样'
__author__ = '皮'
__mtime__ = '4/8/2016-008'
__email__ = 'pipisorry@126.com'
# code is far away from bugs with the god animal protecting
    I love animals. They taste delicious.
"""
from scipy import optimize
from scipy.stats import norm, uniform
import numpy as np
import matplotlib.pyplot as plt

trunc = [0, 4]  # 实际分布截断坐标点
p = [1, 1]  # 实际分布参数（均值，标准差）
q = [0, 2]  # 建议分布参数（均值，标准差）


def k(p, q, trunc):
    '''
    求k = max_xp(x)/q(x)
    '''
    trunc_factor = (norm.cdf(trunc[1], p[0], p[1]) - norm.cdf(trunc[0], p[0], p[1])) / p[1]
    print(trunc_factor)
    exp_k = lambda x: norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1]) / trunc_factor
    # exp_k = lambda x: norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1])#未截断的
    max_x = optimize.minimize_scalar(lambda x: -norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1]), bounds=trunc)['x']
    k = exp_k(max_x)
    print("max_x = \\nk = \\n".format(max_x, k))
    return k, max_x


def show_k(k, max_x, p, q, trunc):
    '''
    求出k后绘制建议分布概率密度和实际分布概率密度图，看p(x)和k*q(x)是否相切
    '''
    # x = np.linspace(norm.ppf(0.01, loc=p[0], scale=p[1]), norm.ppf(0.99, loc=p[0], scale=p[1]), N)
    x = np.linspace(trunc[0], trunc[1], 100)
    q = k * norm.pdf(x, loc=q[0], scale=q[1])  # 建议分布概率密度
    p = norm.pdf(x, loc=p[0], scale=p[1]) / (
        norm.cdf(trunc[1], p[0], p[1]) - norm.cdf(trunc[0], p[0], p[1]))  # 实际分布概率密度
    plt.plot(x, q, 'r')
    plt.plot(x, p, 'g')
    plt.axvline(max_x, color='b', label=max_x)  # 相切点
    plt.text(max_x, 0, str(round(max_x, 2)))
    plt.show()


def acc_rej_sample(k, p, q, trunc, N):
    '''
    接受拒绝采样
    :param N: 采样数
    '''
    z = norm.rvs(loc=q[0], scale=q[1], size=N)  # 从建议分布采样
    mu = uniform.rvs(size=N)  # 从均匀分布采样
    z = z[(mu <= norm.pdf(z, p[0], p[1]) / (k * norm.pdf(z, q[0], q[1])))]  # 接受-拒绝采样
    z = z[z >= trunc[0]]
    z = z[z <= trunc[1]]
    # print("sampled z = \\n\\n".format(z))
    return z


def show_z(z, p, trunc):
    '''
    采样得到采样样本z后看是否采样得到实际正态分布的近似
    '''
    # 采样分布概率密度图
    cnts, bins = np.histogram(z, bins=500, normed=True)
    bins = (bins[:-1] + bins[1:]) / 2
    plt.plot(bins, cnts, label='sampling dist')
    # plt.hist(z, bins=500, normed=True)

    # 实际分布概率密度图（截断后的）
    x = np.linspace(trunc[0], trunc[1], 100)
    trunc_factor = (norm.cdf(trunc[1], p[0], p[1]) - norm.cdf(trunc[0], p[0], p[1])) / p[1]
    plt.plot(x, norm.pdf(x, loc=p[0], scale=p[1]) / trunc_factor, 'r', label='real dist')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    k, max_x = k(p, q, trunc)
    # show_k(k, max_x, p, q, trunc)
    z = acc_rej_sample(k, p, q, trunc, N=10000000)
    show_z(z, p, trunc)

截断面积trunc_factor：0.839994848037

k最大时的x点max_x = 1.333333333395553

求解出的k值：k = 2.812780139369929

根据求解出的k值绘制的实际截断分布和建议分布相切图：