快速任意分布随机采样（逆变换采样）

Question

random 模块 (http://docs.python.org/2/library/random.html) 有几个固定 函数可以从中随机采样。例如random.gauss 将从具有给定均值和 sigma 值的正态分布中随机抽取点。

我正在寻找一种方法，使用我自己在python 中的分布尽可能快 提取给定间隔之间的随机样本数N。这就是我的意思：

def my_dist(x):
    # Some distribution, assume c1,c2,c3 and c4 are known.
    f = c1*exp(-((x-c2)**c3)/c4)
    return f

# Draw N random samples from my distribution between given limits a,b.
N = 1000
N_rand_samples = ran_func_sample(my_dist, a, b, N)

ran_func_sample 是我所追求的，a, b 是从中抽取样本的限制。 python 里有这种东西吗？

Answer 1

您需要使用逆变换采样方法来获得根据您想要的规律分布的随机值。使用这种方法你可以应用 反函数 在区间 [0,1] 内具有标准均匀分布的随机数。

找到反转函数后，你会得到 1000 个根据需要分布的数字，这种方式很明显：

[inverted_function(random.random()) for x in range(1000)]

更多关于逆变换采样：

• http://en.wikipedia.org/wiki/Inverse_transform_sampling

另外，StackOverflow 上有一个与主题相关的好问题：

• Pythonic way to select list elements with different probability

Answer 2

此代码实现了 n-d 离散概率分布的采样。通过在对象上设置一个标志，它也可以用作分段常数概率分布，然后可以用来逼近任意 pdf。好吧，具有紧凑支持的任意 pdf；如果您想有效地对极长的尾巴进行采样，则需要对 pdf 进行不统一的描述。但是即使对于像airy-point-spread函数（我最初创建它）这样的东西，这仍然是有效的。值的内部排序对于获得准确性是绝对关键的；尾部的许多小值应该有很大的贡献，但它们会在没有排序的情况下被 fp 精度淹没。

class Distribution(object):
    """
    draws samples from a one dimensional probability distribution,
    by means of inversion of a discrete inverstion of a cumulative density function

    the pdf can be sorted first to prevent numerical error in the cumulative sum
    this is set as default; for big density functions with high contrast,
    it is absolutely necessary, and for small density functions,
    the overhead is minimal

    a call to this distibution object returns indices into density array
    """
    def __init__(self, pdf, sort = True, interpolation = True, transform = lambda x: x):
        self.shape          = pdf.shape
        self.pdf            = pdf.ravel()
        self.sort           = sort
        self.interpolation  = interpolation
        self.transform      = transform

        #a pdf can not be negative
        assert(np.all(pdf>=0))

        #sort the pdf by magnitude
        if self.sort:
            self.sortindex = np.argsort(self.pdf, axis=None)
            self.pdf = self.pdf[self.sortindex]
        #construct the cumulative distribution function
        self.cdf = np.cumsum(self.pdf)
    @property
    def ndim(self):
        return len(self.shape)
    @property
    def sum(self):
        """cached sum of all pdf values; the pdf need not sum to one, and is imlpicitly normalized"""
        return self.cdf[-1]
    def __call__(self, N):
        """draw """
        #pick numbers which are uniformly random over the cumulative distribution function
        choice = np.random.uniform(high = self.sum, size = N)
        #find the indices corresponding to this point on the CDF
        index = np.searchsorted(self.cdf, choice)
        #if necessary, map the indices back to their original ordering
        if self.sort:
            index = self.sortindex[index]
        #map back to multi-dimensional indexing
        index = np.unravel_index(index, self.shape)
        index = np.vstack(index)
        #is this a discrete or piecewise continuous distribution?
        if self.interpolation:
            index = index + np.random.uniform(size=index.shape)
        return self.transform(index)


if __name__=='__main__':
    shape = 3,3
    pdf = np.ones(shape)
    pdf[1]=0
    dist = Distribution(pdf, transform=lambda i:i-1.5)
    print dist(10)
    import matplotlib.pyplot as pp
    pp.scatter(*dist(1000))
    pp.show()

作为一个更现实的例子：

x = np.linspace(-100, 100, 512)
p = np.exp(-x**2)
pdf = p[:,None]*p[None,:]     #2d gaussian
dist = Distribution(pdf, transform=lambda i:i-256)
print dist(1000000).mean(axis=1)    #should be in the 1/sqrt(1e6) range
import matplotlib.pyplot as pp
pp.scatter(*dist(1000))
pp.show()

Answer 3

import numpy as np
import scipy.interpolate as interpolate

def inverse_transform_sampling(data, n_bins, n_samples):
    hist, bin_edges = np.histogram(data, bins=n_bins, density=True)
    cum_values = np.zeros(bin_edges.shape)
    cum_values[1:] = np.cumsum(hist*np.diff(bin_edges))
    inv_cdf = interpolate.interp1d(cum_values, bin_edges)
    r = np.random.rand(n_samples)
    return inv_cdf(r)

因此，如果我们提供具有特定分布的数据样本，inverse_transform_sampling 函数将返回具有完全相同分布的数据集。这里的优点是我们可以通过在n_samples 变量中指定我们自己的样本大小。

Answer 4

我遇到了类似的情况，但我想从多元分布中采样，所以我实现了 Metropolis-Hastings 的基本版本（这是一种 MCMC 方法）。

def metropolis_hastings(target_density, size=500000):
    burnin_size = 10000
    size += burnin_size
    x0 = np.array([[0, 0]])
    xt = x0
    samples = []
    for i in range(size):
        xt_candidate = np.array([np.random.multivariate_normal(xt[0], np.eye(2))])
        accept_prob = (target_density(xt_candidate))/(target_density(xt))
        if np.random.uniform(0, 1) < accept_prob:
            xt = xt_candidate
        samples.append(xt)
    samples = np.array(samples[burnin_size:])
    samples = np.reshape(samples, [samples.shape[0], 2])
    return samples

这个函数需要一个函数target_density，它接收一个数据点并计算它的概率。

详情请查看我的detailed answer。

Answer 5

这是使用装饰器执行inverse transform sampling 的一种相当不错的方式。

import numpy as np
from scipy.interpolate import interp1d

def inverse_sample_decorator(dist):
    
    def wrapper(pnts, x_min=-100, x_max=100, n=1e5, **kwargs):
        
        x = np.linspace(x_min, x_max, int(n))
        cumulative = np.cumsum(dist(x, **kwargs))
        cumulative -= cumulative.min()
        f = interp1d(cumulative/cumulative.max(), x)
        return f(np.random.random(pnts))
    
    return wrapper

在高斯分布上使用这个装饰器，例如：

@inverse_sample_decorator
def gauss(x, amp=1.0, mean=0.0, std=0.2):
    return amp*np.exp(-(x-mean)**2/std**2/2.0)

然后您可以通过调用该函数从分布中生成样本点。关键字参数x_min 和x_max 是原始分布的限制，可以与参数化分布的其他关键字参数一起作为参数传递给gauss。

samples = gauss(5000, mean=20, std=0.8, x_min=19, x_max=21)

或者，这可以作为将分布作为参数的函数来完成（如您的原始问题中所示），

def inverse_sample_function(dist, pnts, x_min=-100, x_max=100, n=1e5, 
                            **kwargs):
        
    x = np.linspace(x_min, x_max, int(n))
    cumulative = np.cumsum(dist(x, **kwargs))
    cumulative -= cumulative.min()
    f = interp1d(cumulative/cumulative.max(), x)
        
    return f(np.random.random(pnts))