scipy.stats.sampling.

RatioUniforms#

class scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)[Quelle]#

Generiere Zufallsstichproben aus einer Wahrscheinlichkeitsdichtefunktion mittels der Ratio-of-Uniforms-Methode.

Parameter:
pdfcallable

Eine Funktion mit der Signatur pdf(x), die proportional zur Wahrscheinlichkeitsdichtefunktion der Verteilung ist.

umaxfloat

Die obere Grenze des umschließenden Rechtecks in u-Richtung.

vminfloat

Die untere Grenze des umschließenden Rechtecks in v-Richtung.

vmaxfloat

Die obere Grenze des umschließenden Rechtecks in v-Richtung.

cfloat, optional.

Verschiebungsparameter der Ratio-of-Uniforms-Methode, siehe Hinweise. Standardwert ist 0.

random_state{None, int, numpy.random.Generator,

Wenn seed None ist (oder np.random), wird die Singleton-Instanz von numpy.random.RandomState verwendet. Wenn seed eine Ganzzahl ist, wird eine neue RandomState-Instanz mit seed initialisiert. Wenn seed bereits eine Generator- oder RandomState-Instanz ist, wird diese Instanz verwendet.

Methoden

rvs([size])

Stichprobengenerierung von Zufallsvariablen

Hinweise

Gegeben eine univariaten Wahrscheinlichkeitsdichtefunktion pdf und eine Konstante c, definiere die Menge A = {(u, v) : 0 < u <= sqrt(pdf(v/u + c))}. Wenn (U, V) ein Zufallsvektor ist, der gleichmäßig über A verteilt ist, dann folgt V/U + c einer Verteilung gemäß pdf.

Das obige Ergebnis (siehe [1], [2]) kann verwendet werden, um Zufallsvariablen nur anhand der PDF zu generieren, d.h. es ist keine Inversion der CDF erforderlich. Typische Wahlmöglichkeiten für c sind Null oder der Modus von pdf. Die Menge A ist eine Teilmenge des Rechtecks R = [0, umax] x [vmin, vmax], wobei

  • umax = sup sqrt(pdf(x))

  • vmin = inf (x - c) sqrt(pdf(x))

  • vmax = sup (x - c) sqrt(pdf(x))

Insbesondere sind diese Werte endlich, wenn pdf beschränkt ist und x**2 * pdf(x) beschränkt ist (d.h. subquadratische Schwänze). Man kann (U, V) gleichmäßig auf R generieren und V/U + c zurückgeben, wenn (U, V) auch in A liegt, was direkt überprüft werden kann.

Der Algorithmus wird nicht verändert, wenn man pdf durch k * pdf für eine beliebige Konstante k > 0 ersetzt. Daher ist es oft praktisch, mit einer Funktion zu arbeiten, die proportional zur Wahrscheinlichkeitsdichtefunktion ist, indem unnötige Normierungsfaktoren weggelassen werden.

Intuitiv funktioniert die Methode gut, wenn A den größten Teil des umschließenden Rechtecks ausfüllt, so dass die Wahrscheinlichkeit hoch ist, dass (U, V) in A liegt, wann immer es in R liegt, da die Anzahl der erforderlichen Iterationen sonst zu groß wird. Genauer gesagt ist die erwartete Anzahl von Iterationen, um (U, V) gleichmäßig auf R zu ziehen, so dass (U, V) auch in A liegt, gegeben durch das Verhältnis area(R) / area(A) = 2 * umax * (vmax - vmin) / area(pdf), wobei area(pdf) das Integral von pdf ist (was gleich eins ist, wenn die Wahrscheinlichkeitsdichtefunktion verwendet wird, aber andere Werte annehmen kann, wenn eine Funktion proportional zur Dichte verwendet wird). Die Gleichheit gilt, da die Fläche von A gleich 0.5 * area(pdf) ist (Theorem 7.1 in [1]). Wenn die Stichprobengenerierung nach 50000 Iterationen keine einzige Zufallsvariable generiert (d.h. kein einziger Zug liegt in A), wird eine Ausnahme ausgelöst.

Wenn das umschließende Rechteck nicht korrekt spezifiziert ist (d.h. wenn es A nicht enthält), generiert der Algorithmus Stichproben aus einer anderen Verteilung als der von pdf gegebenen. Es wird daher empfohlen, einen Test wie kstest zur Überprüfung durchzuführen.

Referenzen

[1] (1,2)

L. Devroye, „Non-Uniform Random Variate Generation“, Springer-Verlag, 1986.

[2]

W. Hoermann und J. Leydold, „Generating generalized inverse Gaussian random variates“, Statistics and Computing, 24(4), S. 547–557, 2014.

[3]

A.J. Kinderman und J.F. Monahan, „Computer Generation of Random Variables Using the Ratio of Uniform Deviates“, ACM Transactions on Mathematical Software, 3(3), S. 257–260, 1977.

Beispiele

>>> import numpy as np
>>> from scipy import stats

>>> from scipy.stats.sampling import RatioUniforms
>>> rng = np.random.default_rng()

Simulieren von normalverteilten Zufallsvariablen. Es ist in diesem Fall einfach, das umschließende Rechteck explizit zu berechnen. Der Einfachheit halber lassen wir den Normierungsfaktor der Dichte weg.

>>> f = lambda x: np.exp(-x**2 / 2)
>>> v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
>>> umax = np.sqrt(f(0))
>>> gen = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=rng)
>>> r = gen.rvs(size=2500)

Der K-S-Test bestätigt, dass die Zufallsvariablen tatsächlich normalverteilt sind (Normalität wird auf dem 5%-Signifikanzniveau nicht verworfen)

>>> stats.kstest(r, 'norm')[1]
0.250634764150542

Die Exponentialverteilung bietet ein weiteres Beispiel, bei dem das umschließende Rechteck explizit bestimmt werden kann.

>>> gen = RatioUniforms(lambda x: np.exp(-x), umax=1, vmin=0,
...                     vmax=2*np.exp(-1), random_state=rng)
>>> r = gen.rvs(1000)
>>> stats.kstest(r, 'expon')[1]
0.21121052054580314
Auf dieser Seite