随机采样和随机模拟

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)的样本生成。

直接采样步骤:

  1. 从Uniform(0,1)随机产生一个样本z
  2. 另z=h(y),其中h(y)为y的累积概率分布CDF
  3. 计算y=h−1(z)
  4. 结果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_x{p(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值绘制的实际截断分布和建议分布相切图:

采样数据:

[ 1.47996773 -0.55490135  2.63931978 ...,  1.43872433 0.797398    0.53202651]

采样结果同真实结果的比较:

接受-拒绝采样的matlab实现

Note:这个代码也是没有考虑截断面积的!

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

ref: prml chap11

 

相关推荐
©️2020 CSDN 皮肤主题: 编程工作室 设计师:CSDN官方博客 返回首页