逆变换采样 (inverse transform sampling) 的原理

前文介绍了，对随机变量做函数变换 $Y=f(X)$ 后的概率密度函数 (PDF) 之间的变化：
$$P_Y(y)=P_X(f^{-1}(y))\left|\frac{d{f}^{-1}(y)}{dy}\right|=P_X(x)\left|\frac{dx}{dy}\right|$$

在这先声明一下记法：

• $Pr(\cdot )$ 表示概率；
• $P_X(x)$ 表示随机变量 $X$ 的概率密度函数；
• $F_X(x)$ 表示随机变量 $X$ 的累计分布函数。

借助该引理，可以得到对任意分布采样的逆变换采样方法。

假设变量变换函数 $Y=f(X)$ 正好是变量 $X$ 的累积分布函数
$$f(x)=F_X(x)=Pr(X\le x)={\int }_{-\infty }^x{P}_X(s)ds$$

那么变换后 $Y=f(X)$ 的概率密度函数为：
$$P_Y(y)=P_X(x)\left|\frac{dx}{dy}\right|=P_X(x)|f^{\prime }(x){|}^{-1}=\frac{P_X(x)}{|P_X(x)|}=1$$
且 $Y$ 的值域为 $[0,1]$.

所以以累积分布函数作变换后得到的随机变量 $Y$ 服从 $[0,1]$ 上的均匀分布！！！

反之，如果先从 $[0,1]$ 上的均匀分布上随机采样 $y_i$，再做某个累计分布函数的逆变换 $x_i=F_X^{-1}(y_i)$，得到的 $x_i$ 的累积分布函数正好是 $F_X(x)$ ！！！

这就是实现了对任意分布 $F_X(x)$ 采样。

其实逆变换采样和轮盘赌采样是一回事~

code

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


class TreeSampling(object):
    """
    Sampling from a large population
    Construct in O(N) time
    Sample and update in O(log(N)) time
    """

    def __init__(self, dimension, weights=None):
        self.dimension = dimension
        self.layers = int(np.ceil(np.log2(dimension)))
        self.F = [np.array([])] * self.layers

        self.initialize(weights)

    def initialize(self, weights=None):
        """
        initialize F+ tree with uniform weights
        """
        # initialzie last layer with weights
        if weights is None:
            weight = 1.0 / self.dimension
            self.F[-1] = np.ones((self.dimension,)) * weight
        else:
            self.F[-1] = weights

        for l in range(self.layers - 2, -1, -1):
            length = int(np.ceil(self.F[l + 1].shape[0] / 2.0))
            self.F[l] = np.ones((length,))
            if len(self.F[l + 1]) % 2 != 0:
                self.F[l][:-1] = self.F[l + 1][:-1].reshape((-1, 2)).sum(axis=1)
                self.F[l][-1] = self.F[l + 1][-1]
            else:
                self.F[l] = self.F[l + 1].reshape((-1, 2)).sum(axis=1)

    def print_graph(self):
        if self.dimension > 1000:
            print("Are you crazy?")
            return
        for fl in self.F:
            for prob in fl:
                print(prob, end=" ")
            print("||")

    def total_weight(self):
        """
        return the total weight sum
        """
        return self.F[0][0] + self.F[0][1]

    def get_weight(self, indices):
        """
        return the weight of given indices
        """
        return self.F[-1][indices]

    def sample_batch(self, batch_size):
        """
        sample a batch without replacement
        """
        indices = np.zeros((batch_size,), dtype=np.int)
        weights = np.zeros((batch_size,), dtype=np.float)
        for i in range(batch_size):
            indices[i] = self.__sample()
            weights[i] = self.F[-1][indices[i]]
            self.__update(indices[i], 0)  # wighout replacement
        self.update_batch(indices, weights)  # resume their original weights
        return indices

    def update_batch(self, indices, probs):
        """
        update weights of a given batch
        """
        for i, p in zip(indices, probs):
            self.__update(i, p)

    def __sample(self):
        """
        sample a single node, in log(N) time
        """
        u = np.random.sample() * self.total_weight()
        i = 0
        for fl in self.F:
            # i_left = 2*i
            # i_right = 2*i +1
            if u > fl[2 * i] and fl.shape[0] >= 2 * (i + 1):  # then chose i_right
                u -= fl[2 * i]
                i = 2 * i + 1
            else:
                i = 2 * i
        return i

    def __update(self, idx, prob):
        """
        update weight of a single node, in log(N) time
        """
        delta = prob - self.F[-1][idx]

        for l in range(self.layers - 1, -1, -1):
            self.F[l][idx] += delta
            idx = idx // 2

N = 10000
idx = np.array([i for i in range(N)])
weights = np.square(np.sin(0.001*idx))
plt.plot(weights)
f = TreeSampling(N, weights)  

samples = []
for i in range(100):
    samples += list(f.sample_batch(batch_size=100))

_ = plt.hist(samples, bins=100)

相关知识点还可以参考概率积分变换（Probability Integral transform）