Schnelle beliebige Verteilung zufällige Stichprobe

Question

Das random Modul (http://docs.python.org/2/library/random.html) hat mehrere feste Funktionen, um zufällig Stichproben aus. Zum Beispiel wird random.gauss zufälligen Punkt von einer Normalverteilung mit einem gegebenen Mittelwert und Sigma-Werten abtasten.

Ich bin auf der Suche nach einer Möglichkeit, eine Reihe N von Stichproben zwischen einem bestimmten Intervall mit meiner eigenen Verteilung so schnell wie möglich in python zu extrahieren. Das ist, was ich meine:

def my_dist(x): 
    # Some distribution, assume c1,c2,c3 and c4 are known. 
    f = c1*exp(-((x-c2)**c3)/c4) 
    return f 

# Draw N random samples from my distribution between given limits a,b. 
N = 1000 
N_rand_samples = ran_func_sample(my_dist, a, b, N) 

wo ran_func_sample ist, was ich nach und a, b sind die Grenzen, um die Proben zu ziehen. Gibt es so etwas in python?

Answer 1

Sie müssen Inverse Transformation Stichprobenverfahren Methode verwenden, um zufällige Werte nach einem Gesetz verteilt erhalten Sie wollen. Mit dieser Methode können Sie nur invertierte Funktion zu Zufallszahlen mit Standardgleichmäßige Verteilung im Intervall [0,1] anwenden können.

Nachdem Sie die umgekehrte Funktion finden, erhalten Sie 1000 Nummern verteilt entsprechend der benötigten Verteilung diese offensichtliche Art und Weise:

[inverted_function(random.random()) for x in range(1000)] 

Mehr auf inverse Transformation Sampling:

• http://en.wikipedia.org/wiki/Inverse_transform_sampling

Auch gibt es eine gute Frage zu StackOverflow Relat ed zum Thema:

• Pythonic way to select list elements with different probability

Answer 2

Dieser Code implementiert die Probennahme von n-d diskreten Wahrscheinlichkeitsverteilungen. Indem ein Flag auf das Objekt gesetzt wird, kann es auch als eine stückweise konstante Wahrscheinlichkeitsverteilung verwendet werden, die dann verwendet werden kann, um beliebige PDFs zu approximieren. Nun, beliebige PDFs mit kompakter Unterstützung; Wenn Sie extrem lange Tails effizient testen möchten, ist eine nicht einheitliche Beschreibung der PDF erforderlich. Aber auch für Dinge wie luftig-spreizende Funktionen (für die ich sie ursprünglich geschaffen habe) ist das immer noch effizient. Die interne Sortierung von Werten ist dort absolut kritisch, um Genauigkeit zu erhalten; die vielen kleinen Werte in den Schwänzen sollten wesentlich beitragen, aber sie werden in fp Genauigkeit ohne Sortierung übertönt werden. relevantes Beispiel

class Distribution(object): 
    """ 
    draws samples from a one dimensional probability distribution, 
    by means of inversion of a discrete inverstion of a cumulative density function 

    the pdf can be sorted first to prevent numerical error in the cumulative sum 
    this is set as default; for big density functions with high contrast, 
    it is absolutely necessary, and for small density functions, 
    the overhead is minimal 

    a call to this distibution object returns indices into density array 
    """ 
    def __init__(self, pdf, sort = True, interpolation = True, transform = lambda x: x): 
     self.shape   = pdf.shape 
     self.pdf   = pdf.ravel() 
     self.sort   = sort 
     self.interpolation = interpolation 
     self.transform  = transform 

     #a pdf can not be negative 
     assert(np.all(pdf>=0)) 

     #sort the pdf by magnitude 
     if self.sort: 
      self.sortindex = np.argsort(self.pdf, axis=None) 
      self.pdf = self.pdf[self.sortindex] 
     #construct the cumulative distribution function 
     self.cdf = np.cumsum(self.pdf) 
    @property 
    def ndim(self): 
     return len(self.shape) 
    @property 
    def sum(self): 
     """cached sum of all pdf values; the pdf need not sum to one, and is imlpicitly normalized""" 
     return self.cdf[-1] 
    def __call__(self, N): 
     """draw """ 
     #pick numbers which are uniformly random over the cumulative distribution function 
     choice = np.random.uniform(high = self.sum, size = N) 
     #find the indices corresponding to this point on the CDF 
     index = np.searchsorted(self.cdf, choice) 
     #if necessary, map the indices back to their original ordering 
     if self.sort: 
      index = self.sortindex[index] 
     #map back to multi-dimensional indexing 
     index = np.unravel_index(index, self.shape) 
     index = np.vstack(index) 
     #is this a discrete or piecewise continuous distribution? 
     if self.interpolation: 
      index = index + np.random.uniform(size=index.shape) 
     return self.transform(index) 


if __name__=='__main__': 
    shape = 3,3 
    pdf = np.ones(shape) 
    pdf[1]=0 
    dist = Distribution(pdf, transform=lambda i:i-1.5) 
    print dist(10) 
    import matplotlib.pyplot as pp 
    pp.scatter(*dist(1000)) 
    pp.show() 

Und als mehr der realen Welt:

x = np.linspace(-100, 100, 512) 
p = np.exp(-x**2) 
pdf = p[:,None]*p[None,:]  #2d gaussian 
dist = Distribution(pdf, transform=lambda i:i-256) 
print dist(1000000).mean(axis=1) #should be in the 1/sqrt(1e6) range 
import matplotlib.pyplot as pp 
pp.scatter(*dist(1000)) 
pp.show() 

Answer 3

import numpy as np 
import scipy.interpolate as interpolate 

def inverse_transform_sampling(data, n_bins, n_samples): 
    hist, bin_edges = np.histogram(data, bins=n_bins, density=True) 
    cum_values = np.zeros(bin_edges.shape) 
    cum_values[1:] = np.cumsum(hist*np.diff(bin_edges)) 
    inv_cdf = interpolate.interp1d(cum_values, bin_edges) 
    r = np.random.rand(n_samples) 
    return inv_cdf(r) 

Also, wenn wir unsere Datenprobe geben, die eine bestimmte Verteilung hat, wird die inverse_transform_sampling Funktion einen Datensatz mit genau Rückkehr der gleiche Verteilung. Hier ist der Vorteil ist, dass wir unsere eigene Probengröße indem sie sie in der n_samples Variable Angabe bekommen.