np.random.choice 未生成预期的直方图

np.random.choice not producing expected histogram

我希望生成 1 和 0 之间的 random normally distributed 数字,但是随着 mean 接近 1 或 0,右侧或左侧分别被“压扁”。

修改正态分布并在 geogebra 中使用滑块后,我得出以下结论:

接下来我需要在 python 中创建一个方法,该方法将生成根据此 PDF 分发的随机样本。

最初我认为做到这一点的唯一方法是尝试推导一个新的方程来生成随机数,如 Box-Muller 证明中所示(我通过遵循 this 教程获得).

不过,我认为使用 numpy 库的 np.random.choice() 方法可能有更简单的方法。

毕竟,我应该能够以非常小的步长对 PDF 进行积分,并获得所述步骤的各种概率(当然是大约)。

因此我编写了以下脚本:

# Standard libs
import math

# Third party libs
import numpy as np

from alive_progress import alive_bar
from matplotlib import pyplot as plt

class RandomNumberGenerator:
    def __init__(self):
        pass

    def clamped_normal_distribution(self, mu: float, 
            stddev: float, x: float):
        """ Computes a value from the clamped normal distribution """
        divideByZeroAvoider = 1e-5
        if x < 0 or x > 1:
            return 0
        elif x >= 0 and x <= mu:
            return math.exp(-0.5*( (x - mu) / (stddev)  )**2 \
                    * (1/(x**2 + divideByZeroAvoider)))
        elif x <= 1 and x > mu:
            return math.exp(-0.5*( (x - mu) / (stddev)  )**2 \
                    * (1/((1-x)**2 + divideByZeroAvoider))) 
        else:
            print("This shouldn't happen!: {}".format(x))
            return 0

if __name__ == '__main__':
    rng = RandomNumberGenerator()

    mu = 0.7
    stddev = 1
    stepSize = 1e-3
    x = np.linspace(stepSize,1, int(1/stepSize) - 1)

    # Determine the total area under the curve
    samples = []
    print("Generating samples...")
    with alive_bar(len(x.tolist())) as bar:
        for i in x:
            samples.append(rng.clamped_normal_distribution(
                    mu, stddev, i))
            bar()

    area = np.trapz(samples, dx=stepSize)
    print("Area = {}".format(area))

    # Determine the probability of x falling in a specific interval
    probabilities = []

    print("Generating probabilties...")
    with alive_bar(len(x.tolist())) as bar:
        for i in x:
            lead = rng.clamped_normal_distribution(mu, 
                    stddev, i)
            lag = rng.clamped_normal_distribution(mu, 
                    stddev, i - stepSize)
            probability = np.trapz(
                    np.array([lag, lead]), 
                    dx=stepSize)
            
            # Divide by the area because this isn't a standard normal
            probabilities.append(probability / area)
            bar()
    
    # Should be approximately 1
    print("Probability: {}".format(sum(probabilities)))

    plt.plot(x, probabilities)
    plt.show()

    y = []
    print("Performing distribution test...")
    testSize = int(10e3)
    with alive_bar(testSize) as bar:
        for _ in range(testSize):
            randSamp = np.random.choice(samples, p=probabilities)
            y.append(randSamp)
            bar()

    plt.hist(y,300)
    plt.show()

线性间隔样本的第一个概率图看起来很有希望,给出了下图:

但是,如果我们使用这些样本作为具有给定概率的选择,我们会得到以下直方图:

我不知道为什么这不能正常工作。

我尝试了其他(较小的)示例,例如 numpy website 中列出的示例,它们根据给定的概率数组生成了直方图。

如果可能的话,我真的很感激 advice/intuition:).

调用 np.random.choice(samples, p=probabilities) 中的第一个参数似乎有问题。第一个参数应该是 x,而不是 samples.

作者添加:

原因是 samples 是曲线的值(即 y-axis 而不是 x-axis)。

因此,概率最高的值(即平均值附近的样本)的值都约为 1,这就是为什么我们在值 1 附近看到如此巨大的尖峰。

将其更改为 x 会得到以下图表(针对 10e3 个样本):

按预期工作,非常好。