如何實現高斯分佈

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使用的建議。我查看了this維基百科文章和Numpy源文件，發現這個randomkit.c文件實現了函數rk_gauss（實現Box Muller變換），rk_double和rk_random（它實現了模擬均勻分佈隨機變量的Mersenne Twister Random Number Generator， Box Muller變換所需的）。

我然後從here適於Mersenne扭曲發生器和實現的框穆勒變換到模擬高斯（約隨機倍捻機發生器here更多信息）。

這是我最後寫代碼：

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) 

Answer 2

使用箱穆勒方法：

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