Transformierte Dichtevorgabe (TDR)

• Erforderlich: T-konkave PDF, dPDF

• Optional: Modus, Zentrum

• Geschwindigkeit

  • Einrichtung: langsam

  • Sampling: schnell

TDR ist eine Annahme/Zurückweisungs-Methode, die die Konkavität einer transformierten Dichte verwendet, um eine Hut-Funktion zu konstruieren und automatisch zu verengen. Solche PDFs werden als T-konkav bezeichnet. Derzeit sind die folgenden Transformationen implementiert

\[\begin{split}c = 0.: T(x) &= \log(x)\\ c = -0.5: T(x) &= \frac{1}{\sqrt{x}} \text{ (Standard)}\end{split}\]

Zusätzlich zur PDF wird auch die Ableitung der PDF bezüglich x (d. h. der Zufallsvariablen) benötigt. Diese Funktionen müssen als Methoden eines Python-Objekts vorhanden sein, das dann den Generatoren übergeben werden kann, um deren Objekt zu instanziieren. Die implementierte Variante verwendet Verengungen proportional zur Hut-Funktion ([1]).

Ein Beispiel für die Verwendung dieser Methode ist unten gezeigt

>>> import numpy as np
>>> from scipy.stats.sampling import TransformedDensityRejection
>>> from scipy.stats import norm
>>>
>>> class StandardNormal:
...     def pdf(self, x):
...         # note that the normalization constant is not required
...         return np.exp(-0.5 * x*x)
...     def dpdf(self, x):
...         return -x * np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>>
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rng.rvs()
-1.526829048388144

Im obigen Beispiel haben wir die TDR-Methode verwendet, um Stichproben aus der Standardnormalverteilung zu ziehen. Beachten Sie, dass wir die Normierungskonstante bei der Berechnung der PDF weglassen können. Dies beschleunigt normalerweise die Sampling-Phase. Beachten Sie auch, dass die PDF nicht vektorisiert sein muss. Sie sollte einen Skalar akzeptieren und zurückgeben.

Es ist auch möglich, die Umkehrfunktion der CDF der Hut-Verteilung direkt über die Methode ppf_hat auszuwerten.

>>> rng.ppf_hat(0.5)
-0.00018050266342362759
>>> norm.ppf(0.5)
0.0
>>> u = np.linspace(0, 1, num=10)
>>> rng.ppf_hat(u)
array([       -inf, -1.22227372, -0.7656556 , -0.43135731, -0.14002921,
        0.13966423,  0.43096141,  0.76517113,  1.22185606,         inf])
>>> norm.ppf(u)
array([       -inf, -1.22064035, -0.76470967, -0.4307273 , -0.1397103 ,
        0.1397103 ,  0.4307273 ,  0.76470967,  1.22064035,         inf])

Neben der PPF-Methode können auch andere Attribute abgerufen werden, um zu sehen, wie gut der Generator der gegebenen Verteilung entspricht. Dies sind

• „squeeze_hat_ratio“: (Fläche unter der Verengung) / (Fläche unter dem Hut) für den Generator. Es ist eine Zahl zwischen 0 und 1. Näher an 1 bedeutet, dass die Hut- und Verengungsfunktionen die Verteilung eng umschließen und weniger PDF-Auswertungen zur Erzeugung von Stichproben erforderlich sind. Die erwartete Anzahl von Auswertungen der Dichte ist begrenzt durch (1/squeeze_hat_ratio) - 1 pro Stichprobe. Standardmäßig wird sie über 0,99 gehalten, dies kann jedoch durch Übergabe eines Parameters max_squeeze_hat_ratio geändert werden.

• „hat_area“: Fläche unter dem Hut für den Generator.

• „squeeze_area“: Fläche unter der Verengung für den Generator.

  >>> rng.squeeze_hat_ratio
  0.9947024204884917
  >>> rng.hat_area
  2.510253139791547
  >>> rng.squeeze_area
  2.4969548741894876
  >>> rng.squeeze_hat_ratio == rng.squeeze_area / rng.hat_area
  True
  

Die Verteilung kann durch Übergabe eines Domain-Parameters gekürzt werden

>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, domain=[0, 1], random_state=urng)
>>> rng.rvs(10)
array([0.05452512, 0.97251362, 0.49955877, 0.82789729, 0.33048885,
       0.55558548, 0.23168323, 0.13423275, 0.73176575, 0.35739799])

Wenn die Domain nicht angegeben ist, wird die Methode support des Objekts dist verwendet, um die Domain zu bestimmen

>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...     def dpdf(self, x):
...         return -x * np.exp(-0.5 * x*x)
...     def support(self):
...         return -np.inf, np.inf
...
>>> dist = StandardNormal()
>>>
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rng.rvs(10)
array([-1.52682905,  2.06206883,  0.15205036,  1.11587367, -0.30775562,
       0.29879802, -0.61858268, -1.01049115,  0.78853694, -0.23060766])

Wenn das Objekt dist keine support-Methode bereitstellt, wird die Domain als (-np.inf, np.inf) angenommen.

Um das squeeze_hat_ratio zu erhöhen, übergeben Sie max_squeeze_hat_ratio

>>> dist = StandardNormal()
>>> rng = TransformedDensityRejection(dist, max_squeeze_hat_ratio=0.999,
...                                   random_state=urng)
>>> rng.squeeze_hat_ratio
0.999364900465214

Sehen wir uns an, wie sich dies auf die Rückrufe an die PDF-Methode der Verteilung auswirkt

>>> from copy import copy
>>> class StandardNormal:
...     def __init__(self):
...         self.callbacks = 0
...     def pdf(self, x):
...         self.callbacks += 1
...         return np.exp(-0.5 * x*x)
...     def dpdf(self, x):
...         return -x * np.exp(-0.5 * x*x)
...
>>> dist1 = StandardNormal()
>>> urng1 = np.random.default_rng()
>>> urng2 = copy(urng1)
>>> rng1 = TransformedDensityRejection(dist1, random_state=urng1)
>>> dist1.callbacks  # evaluations during setup
139
>>> dist1.callbacks = 0  # don't consider evaluations during setup
>>> rvs = rng1.rvs(100000)
>>> dist1.callbacks  # evaluations during sampling
527
>>> dist2 = StandardNormal()
>>> # use the same stream of uniform random numbers
>>> rng2 = TransformedDensityRejection(dist2, max_squeeze_hat_ratio=0.999,
...                                    random_state=urng2)
>>> dist2.callbacks  # evaluations during setup
467
>>> dist2.callbacks = 0  # don't consider evaluations during setup
>>> rvs = rng2.rvs(100000)
>>> dist2.callbacks  # evaluations during sampling
84      # may vary

Wie wir sehen können, sind bei Erhöhung des squeeze_hat_ratio während des Samplings deutlich weniger PDF-Auswertungen erforderlich. Die PPF-Hut-Funktion ist ebenfalls genauer

>>> abs(norm.ppf(0.975) - rng1.ppf_hat(0.975))
0.0027054565421578136
>>> abs(norm.ppf(0.975) - rng2.ppf_hat(0.975))
0.00047824084476300044

Beachten Sie jedoch, dass dies auf Kosten erhöhter PDF-Auswertungen während der Einrichtung geht.

Für Dichten mit Modi, die nicht nahe 0 liegen, wird empfohlen, entweder den Modus oder das Zentrum der Verteilung festzulegen, indem die Parameter mode oder center übergeben werden. Letzteres ist die ungefähre Lage des Modus oder des Mittelwerts der Verteilung. Diese Lage liefert einige Informationen über den Hauptteil der PDF und wird verwendet, um numerische Probleme zu vermeiden.

>>> # mode = 0 for our distribution
>>> # if exact mode is not available, pass 'center' parameter instead
>>> rng = TransformedDensityRejection(dist, mode=0.)

Standardmäßig verwendet die Methode 30 Konstruktionspunkte, um den Hut zu konstruieren. Dies kann durch Übergabe eines Parameters construction_points geändert werden, der entweder ein Array von Konstruktionspunkten oder eine Ganzzahl sein kann, die die Anzahl der zu verwendenden Konstruktionspunkte darstellt.

>>> rng = TransformedDensityRejection(dist,
...                                   construction_points=[-5, 0, 5])

Diese Methode akzeptiert viele weitere Einrichtungsparameter. Eine vollständige Liste finden Sie in der Dokumentation. Weitere Informationen zu den Parametern und der Methode finden Sie in Abschnitt 5.3.16 des UNU.RAN-Benutzerhandbuchs.

Weitere Einzelheiten zu dieser Methode finden Sie unter [1] und [2].

Referenzen

[1] (1,2)

UNU.RAN Referenzhandbuch, Abschnitt 5.3.16, „TDR - Transformed Density Rejection“, http://statmath.wu.ac.at/software/unuran/doc/unuran.html#TDR

[2]

Hörmann, Wolfgang. „A rejection technique for sampling from T-concave distributions.“ ACM Transactions on Mathematical Software (TOMS) 21.2 (1995): 182-193