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

最后我使用了@sascha的建议。我查看了this维基百科文章和Numpy源代码，发现了这个randomkit.c文件，它实现了函数rk_gauss（它实现了Box-Muller变换）、rk_double和{}（它实现了模拟均匀分布随机变量的Mersenne Twister随机数生成器，由Box-Muller变换所需）。在

然后，我从here改编了Mersenne扭曲器生成器，并实现了Box-Muller变换来模拟高斯（关于随机扭曲生成器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

这是完美的工作，并产生一个相当好的高斯

^{pr2}$