采样算法

2019-05-12 Sunday     linux , misc

经常会在各类的算法中看到 “蒙特卡罗” (Monte Carlo) 的字样,例如 Markov Chain Monte Carlo, MCMC、以及 AlphaGo 使用的蒙特卡洛搜索树。

在使用 Monte-Carlo 时,除了所谓的采样之外,在解决贝叶斯问题的时候,很多地方还使用了积分,这时就是通过 Monte-Carlo 的样本和去近似积分。

而所谓的采样,实际上就是根据某种分布去生成一些数据点,最简单的例如抛硬币、掷骰子等,服从均匀分布。

蒙特卡罗方法

Monte Carlo Method 也就是通过重复随机采样模拟对象的概率与统计的问题,在物理、化学、经济学和信息技术领域均具有广泛应用。

它诞生于上个世纪 40 年代美国的 “曼哈顿计划”,主要是由于计算机的发展,其名字来源于赌城蒙特卡罗,象征概率。

就是通过大量的随机样本,去了解一个系统,进而得到所要计算的值。对于许多问题来说,它往往是最简单的计算方法,有时甚至是唯一可行的方法。

简单来说,就是通过随机采样的方式,以频率估计概率。

其实,”蒙特卡罗” 并非是一个算法,而是一个思想或者方法的统称。

PI 计算

假设要估计圆周率 Pi 的值,选取一个边长为 1 的正方形,在正方形内作一个内切圆,那么可以计算出圆的面积与正方形面积之比为 Pi/4

pi example

在正方形内随机生成大量的点,那么圆形区域内点的个数与所有点的个数之比,就可以认为近似等于 Pi/4

其中采样值越大,越接近真实值 3.141592653589793

import math
import random

inside = 0
total = 1000

for i in range(0, total):
        # generate random x, y in [0, 1].
        x = random.random()
        y = random.random()

        # increment if inside unit circle.
        if math.sqrt(x**2 + y**2) <= 1.0:
                inside += 1
# inside / total = pi / 4
pi = (float(inside) / float(total)) * 4
print(pi)

简单采样

对于简单分布 (高斯、指数、均匀、Gamma等) 的采样,实际上在计算机中都已经实现,可以直接通过 numpy.random 模块生成,但是对于一些复杂的分布,尤其在贝叶斯、概率编程里面,则需要对这些复杂分布进行采样。

例如,对于均匀分布来说,就可以通过随机数发生器来实现,包括了伪随机数发生器,然后再以此为基础实现一些比较复杂的采样。

实际上这些库也基本上是这么做的,不同的分布会有不同的算法。

离散分布

可以把概率分布看做一个区间段,然后判断随机数 u 落在哪个区间段内,而区间段的长度和概率成正比,这样采样完全符合原来的分布。

例如有样本空间 {A, B, C, D} ,对应的概率为 p(x)=[0.1, 0.5, 0.3, 0.1]

import numpy as np

def sample_discrete(vec):
	u = np.random.rand()
	start = 0
	for i, num in enumerate(vec):
		if u > start:
			start += num
		else:
			return i-1
	return i

sample_space = ["A", "B", "C", "D"]
sample = dict([(w, 0) for w in sample_space])
for i in range(10000):
    s = sample_discrete([0.1, 0.5, 0.3, 0.1])
    sample[sample_space[s]] += 1
for k in sample:
    print("%s: %d" % (k, sample[k]))

最终输出的内容类似如下,基本上满足上述的概率分布。

A: 991
C: 3060
B: 4954
D: 995

正态分布

一般来说会使用 Box-Muller 变换,也就是通过服从均匀分布的随机变量,来构建服从高斯分布的随机变量的一种方法。

简单来说,选取两个服从 [0, 1] 上均匀分布的随机变量 U1U2 ,满足如下公示,那么 XY 满足均值为 0 ,方差为 1 的高斯分布。

公示推导

对于正态分布函数 (PDF) 有如下公式,其中 $\mu$ 是平均值,$\sigma$ 是标准差。

假设有两个随机变量

import numpy as np
import matplotlib.pyplot as plt

def norm_sample(size, mu, sigma):
    u1 = np.random.rand(size)
    u2 = np.random.rand(size)

    x = np.sqrt(-2*np.log(u1))*np.cos(2*np.pi*u2)
    y = np.sqrt(-2*np.log(u1))*np.sin(2*np.pi*u2)

    x1 = mu + x * sigma
    y1 = mu + y * sigma

    return x1, y1

def norm_pdf(x, mu, sigma):
    return np.exp(-((x - mu)**2)/(2*(sigma**2)))/(sigma*np.sqrt(2*np.pi))

size = int(1e+05)
mu, sigma = 0, 1
x = np.arange(-4., 4., 0.1)

(y1, y2) = norm_sample(size, mu, sigma)

f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
ax1.plot(x, norm_pdf(x, mu, sigma), 'r', lw=2)
ax1.hist(y1, bins=200, density=True)

ax2.plot(x, norm_pdf(x, mu, sigma), 'r', lw=2)
ax2.hist(y2, bins=200, density=True)

plt.show()

输出的内容如下,在采样之后,两者基本都可以满足正态分布。

sample normal distribution

参考



如果喜欢这里的文章,而且又不差钱的话,欢迎打赏个早餐 ^_^


About This Blog

Recent Posts

Categories

Related Links

  • RTEMS
    RTEMS
  • GNU
  • Linux Kernel
  • Arduino

Search


This Site was built by Jin Yang, generated with Jekyll, and hosted on GitHub Pages
©2013-2019 – Jin Yang