如何实现高斯分布的概率密度函数

Question

我需要在 Python 中实现一个类，它代表单变量（目前）正态分布。我的想法如下

class Norm():
    def __init__(self, mu=0, sigma_sq=1):
        self.mu = mu
        self.sigma_sq = sigma_sq
        # some initialization if necessary

    def sample(self):
        # generate a sample, where the probability of the value 
        # of the sample being generated is distributed according
        # a normal distribution with a particular mean and variance
        pass

N = Norm()
N.sample()

生成的样本应根据以下概率密度函数分布

我知道这一点scipy.stats并Numpy提供了执行此操作的函数，但我需要了解这些函数是如何实现的。任何帮助将不胜感激，谢谢:)

Answer 1

我最终采用了@sascha 的建议。我查看了这篇维基百科文章和 Numpy 源代码，发现这个randomkit.c文件实现了函数rk_gauss（实现了 Box Muller 变换），rk_double并且rk_random（实现了模拟均匀分布随机变量的 Mersenne Twister 随机数生成器，需要通过 Box Muller 变换）。

然后，我从此处改编了 Mersenne Twister Generator ，并实现了 Box Muller 变换来模拟高斯（有关随机 Twister Generator 的更多信息，请参见此处）。

这是我最终编写的代码：

import numpy as np
from datetime import datetime
import matplotlib.pyplot as plt

class Distribution():
    def __init__(self):
        pass

    def plot(self, number_of_samples=100000):
        # the histogram of the data
        n, bins, patches = plt.hist([self.sample() for i in range(number_of_samples)], 100, normed=1, facecolor='g', alpha=0.75)
        plt.show()

    def sample(self):
        # dummy sample function (to be overridden)
        return 1

class Uniform_distribution(Distribution):
    # Create a length 624 list to store the state of the generator
    MT = [0 for i in xrange(624)]
    index = 0

    # To get last 32 bits
    bitmask_1 = (2 ** 32) - 1

    # To get 32. bit
    bitmask_2 = 2 ** 31

    # To get last 31 bits
    bitmask_3 = (2 ** 31) - 1

    def __init__(self, seed):
        self.initialize_generator(seed)

    def initialize_generator(self, seed):
        "Initialize the generator from a seed"
        global MT
        global bitmask_1
        MT[0] = seed
        for i in xrange(1,624):
            MT[i] = ((1812433253 * MT[i-1]) ^ ((MT[i-1] >> 30) + i)) & bitmask_1

    def generate_numbers(self):
        "Generate an array of 624 untempered numbers"
        global MT
        for i in xrange(624):
            y = (MT[i] & bitmask_2) + (MT[(i + 1 ) % 624] & bitmask_3)
            MT[i] = MT[(i + 397) % 624] ^ (y >> 1)
            if y % 2 != 0:
                MT[i] ^= 2567483615

    def sample(self):
        """
        Extract a tempered pseudorandom number based on the index-th value,
        calling generate_numbers() every 624 numbers
        """
        global index
        global MT
        if index == 0:
            self.generate_numbers()
        y = MT[index]
        y ^= y >> 11
        y ^= (y << 7) & 2636928640
        y ^= (y << 15) & 4022730752
        y ^= y >> 18

        index = (index + 1) % 624
        # divide by 4294967296, which is the largest 32 bit number 
        # to normalize the output value to the range [0,1]
        return y*1.0/4294967296

class Norm(Distribution):
    def __init__(self, mu=0, sigma_sq=1):
        self.mu = mu
        self.sigma_sq = sigma_sq
        self.uniform_distribution_1 = Uniform_distribution(datetime.now().microsecond)
        self.uniform_distribution_2 = Uniform_distribution(datetime.now().microsecond)
        # some initialization if necessary

    def sample(self):
        # generate a sample, where the value of the sample being generated 
        # is distributed according a normal distribution with a particular 
        # mean and variance
        u = self.uniform_distribution_1.sample()
        v = self.uniform_distribution_2.sample()
        return ((self.sigma_sq**0.5)*((-2*np.log(u))**0.5)*np.cos(2*np.pi*v)) + self.mu

这非常有效，并且生成了非常好的高斯分布

Norm().plot(10000)