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

Question

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

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

使用 Box-Müller 方法：

def sample(self):
  x = np.random.uniform(0,1,[2])
  z = np.sqrt(-2*np.log(x[0]))*np.cos(2*np.pi*x[1])
  return z * self.sigma_sq + self.mu

Answer 2

我最终使用了@sascha 的建议。我查看了实现函数 rk_gauss（实现 Box Muller 变换）、rk_double 和 rk_random（实现模拟的 Mersenne Twister 随机数生成器）的 this wikipedia article and the Numpy source and found this randomkit.c 文件Box Muller 变换所需的均匀分布随机变量）。

然后我从 here and implemented the Box Muller Transform to simulate a gaussian (more information about Random Twister Generator here).

改编了 Mersenne 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)

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

How do I implement the Probability density function of a Gaussian Distribution

python

numpy

scipy

probability-density

probability-distribution