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 = {} k = {} ".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 = {} ".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