Comment implémenter la fonction de densité de probabilité d'une distribution gaussienne

Question

Je dois implémenter une classe en Python, qui représente une distribution normale univariée (pour l'instant). Ce que j'ai à l'esprit est comme suit

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() 

Les échantillons générés doivent être distribués en fonction de la fonction de densité de probabilité suivante

Je sais que scipy.stats et Numpy fournissent des fonctions pour ce faire, mais J'ai besoin de comprendre comment ces fonctions sont implémentées. Toute aide serait appréciée, merci :)

Answer 1

J'ai fini utilisant les conseils par @sascha. J'ai regardé l'article wikipedia this et la source Numpy et trouvé ce fichier qui implémente les fonctions rk_gauss (qui implémente la transformée Box Muller), rk_double et rk_random (qui implémente le générateur de nombres aléatoires Mersenne Twister qui simule une variable aléatoire uniformément distribuée, requis par la transformée Box Muller). J'ai ensuite adapté le générateur Mersenne Twister de here et mis en œuvre la boîte Muller Transform pour simuler un gaussien (plus d'informations sur Random Twister Generator here).

Voici le code que je fini par écrire:

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 

Cela fonctionne parfaitement, et génère une assez bonne gaussienne

Norm().plot(10000) 

Answer 2

En utilisant la méthode 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