Fast willkürliche Verteilung (random sampling

Dieser code implementiert die Probenahme des n-d-diskrete Wahrscheinlichkeitsverteilungen. Indem Sie ein flag auf das Objekt, es kann auch verwendet werden, als eine stückweise Konstante Wahrscheinlichkeitsverteilung, die dann verwendet werden können, um annähernd beliebige pdf-Dateien. Gut, beliebige PDF-Dateien mit compact-Unterstützung; wenn Sie effizient probieren will extrem langen Schwanz, eine nicht-einheitliche Beschreibung der pdf-Datei erforderlich wäre. Aber das ist immer noch wirksam, auch für Dinge wie luftig-point-spread-Funktionen (die ich erstellt es für Sie, zunächst). Die interne Sortierung der Werte ist absolut wichtig, es zu bekommen Genauigkeit; die vielen kleinen Werte in den Schwanz, der dazu beitragen sollte deutlich, aber Sie bekommen Sie ertrunken im fp-Genauigkeit ohne Sortierung.

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 relevant Beispiel:

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

InformationsquelleAutor Eelco Hoogendoorn

War ich in einer ähnlichen situation, aber ich wollte Stichprobe aus einer multivariaten Verteilung, also, ich implementiert eine rudimentäre version des Metropolis-Hastings (MCMC-Methode).

def metropolis_hastings(target_density, size=500000):
    burnin_size = 10000
    size += burnin_size
    x0 = np.array([[0, 0]])
    xt = x0
    samples = []
    for i in range(size):
        xt_candidate = np.array([np.random.multivariate_normal(xt[0], np.eye(2))])
        accept_prob = (target_density(xt_candidate))/(target_density(xt))
        if np.random.uniform(0, 1) < accept_prob:
            xt = xt_candidate
        samples.append(xt)
    samples = np.array(samples[burnin_size:])
    samples = np.reshape(samples, [samples.shape[0], 2])
    return samples

Erfordert auch diese Funktion eine Funktion target_density die in einem Daten-Punkt und berechnet die Wahrscheinlichkeit.

Details check-out dieser ausführliche Antwort von mir.

InformationsquelleAutor Abdul Fatir