Polynomial interpolation based INVersion of CDF (PINV)

• Erforderlich: PDF

• Optional: CDF, modus, zentrum

• Geschwindigkeit

  • Einrichtung: (sehr) langsam

  • Stichprobenziehung: (sehr) schnell

Polynomial interpolation based INVersion of CDF (PINV) ist eine Inversionsmethode, die nur die Dichtefunktion benötigt, um Stichproben aus einer Verteilung zu ziehen. Sie basiert auf der Polynominterpolation der PPF und der Gauss-Lobatto-Integration der PDF. Sie bietet Kontrolle über den Interpolationsfehler und den Integrationsfehler. Ihr Hauptzweck ist die Bereitstellung sehr schneller Stichproben, die für jede gegebene Verteilung nahezu gleich sind, auf Kosten einer moderaten bis langsamen Einrichtungszeit. Sie ist die schnellste bekannte Inversionsmethode für den Fall fester Parameter.

Die Inversionsmethode ist die einfachste und flexibelste Methode, um nicht-uniforme Zufallsvariaten zu ziehen. Für eine Zielverteilung mit CDF \(F\) und einer aus \(\text{Uniform}(0, 1)\) gezogenen uniformen Zufallsvariaten \(U\) wird eine Zufallsvariate X generiert, indem die uniforme Zufallsvariate \(U\) mithilfe der PPF (inverse CDF) der Verteilung transformiert wird

\[X = F^{-1}(U)\]

Diese Methode eignet sich aufgrund ihrer Vorteile für stochastische Simulationen. Einige der attraktivsten sind

• Sie bewahrt die strukturellen Eigenschaften des uniformen Zufallszahlengenerators.

• Transformiert eine uniforme Zufallsvariate \(U\) eins zu eins in nicht-uniforme Zufallsvariaten \(X\).

• Einfache und effiziente Stichprobenziehung aus abgeschnittenen Verteilungen.

Leider ist die PPF selten in geschlossener Form verfügbar oder zu langsam, wenn sie verfügbar ist. Für viele Verteilungen ist auch die CDF nicht einfach zu erhalten. Diese Methode adressiert beide Mängel. Der Benutzer muss nur die PDF und optional einen Punkt nahe dem Modus (genannt „Zentrum“) zusammen mit der Größe des maximal zulässigen Fehlers angeben. Sie verwendet eine Kombination aus adaptiver und einfacher Gauss-Lobatto-Quadratur, um die CDF zu erhalten, und Newtons Interpolation, um die PPF zu erhalten. Die Methode ist nicht exakt, da sie nur Zufallsvariaten der approximierten Verteilung erzeugt. Dennoch kann der maximal tolerierte Approximationsfehler nahe der Maschinenpräzision eingestellt werden. Das Konzept des u-Fehlers wird zur Messung und Kontrolle des Fehlers verwendet. Es ist definiert als

\[\epsilon_{u}(u) = | u - F\left(F^{-1}_{a}(u)\right) |\]

wobei \(u \in (0, 1)\) ein Quantil ist, bei dem wir den Fehler messen wollen, und \(F^{-1}_a\) die approximierte PPF der gegebenen Verteilung ist.

Der maximale u-Fehler ist das Kriterium für Approximationsfehler bei der numerischen Berechnung der CDF und PPF. Der maximal tolerierte u-Fehler eines Algorithmus wird als u-Auflösung des Algorithmus bezeichnet und mit \(\epsilon_{u}\) notiert

\[\sup_{u \in (0,1)} | u - F\left(F^{-1}_{a}(u)\right) | \le \epsilon_{u}\]

Der Hauptvorteil des u-Fehlers ist, dass er leicht berechnet werden kann, wenn die CDF verfügbar ist. Für eine detailliertere Diskussion verweisen wir auf [1].

Außerdem funktioniert die Methode nur für beschränkte Verteilungen. Bei unendlichen Schwänzen werden die Enden der Schwänze so abgeschnitten, dass die Fläche darunter kleiner oder gleich \(0.05\epsilon_{u}\) ist.

Es gibt einige Einschränkungen für die gegebene Verteilung

• Der Träger der Verteilung (d. h. der Bereich, in dem die PDF strikt positiv ist) muss zusammenhängend sein. In der Praxis bedeutet dies, dass der Bereich, in dem die PDF „nicht zu klein“ ist, zusammenhängend sein muss. Unimodale Dichten erfüllen diese Bedingung. Wenn diese Bedingung verletzt wird, kann die Domäne der Verteilung abgeschnitten werden.

• Bei der numerischen Integration der PDF muss die gegebene PDF kontinuierlich und glatt sein.

• Die PDF muss beschränkt sein.

• Der Algorithmus hat Probleme, wenn die Verteilung schwere Schwänze hat (da die inverse CDF bei 0 oder 1 sehr steil wird) und die angeforderte u-Auflösung sehr klein ist. Zum Beispiel wird die Cauchy-Verteilung wahrscheinlich dieses Problem zeigen, wenn die angeforderte u-Auflösung kleiner als 1.e-12 ist.

Die folgenden vier Schritte werden vom Algorithmus während der Einrichtung durchgeführt

• Berechnung der Endpunkte der Verteilung: Wenn eine endliche Unterstützung gegeben ist, wird dieser Schritt übersprungen. Andernfalls werden die Enden der Schwänze so abgeschnitten, dass die Fläche darunter kleiner oder gleich \(0.05\epsilon_{u}\) ist.

• Die Domäne wird in Teilintervalle unterteilt, um die CDF und PPF zu berechnen.

• Die CDF wird mittels Gauss-Lobatto-Quadratur berechnet, sodass der Integrationsfehler höchstens \(0.05I_{0}\epsilon_{u}\) beträgt, wobei \(I_{0}\) ungefähr die Gesamtfläche unter der PDF ist.

• Die PPF wird mithilfe der Newtonschen Interpolationsformel mit einem maximalen Interpolationsfehler von \(0.9\epsilon_{u}\) berechnet.

Um den Generator für die Stichprobenziehung aus einer Standardnormalverteilung zu initialisieren, tun Sie Folgendes

>>> import numpy as np
>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)

Der Generator wurde eingerichtet und wir können mit der Stichprobenziehung aus unserer Verteilung beginnen

>>> rng.rvs((5, 3))
array([[-1.52449963,  1.31933688,  2.05884468],
       [ 0.48883931,  0.15207903, -0.02150773],
       [ 1.11486463,  1.95449597, -0.30724928],
       [ 0.9854643 ,  0.29867424,  0.7560304 ],
       [-0.61776203,  0.16033378, -1.00933003]])

Wir können das Histogramm der Zufallsvariaten betrachten, um zu überprüfen, wie gut sie zu unserer Verteilung passen

>>> import matplotlib.pyplot as plt
>>> from scipy.stats import norm
>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)
>>> rvs = rng.rvs(10000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, num=10000)
>>> fx = norm.pdf(x)
>>> plt.plot(x, fx, "r-", label="pdf")
>>> plt.hist(rvs, bins=50, density=True, alpha=0.8, label="rvs")
>>> plt.xlabel("x")
>>> plt.ylabel("PDF(x)")
>>> plt.title("Samples drawn using PINV method.")
>>> plt.legend()
>>> plt.show()

Der maximal tolerierte Fehler (d. h. u-Auflösung) kann geändert werden, indem das Schlüsselwort u_resolution während der Initialisierung übergeben wird

>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
...                                  random_state=urng)

Dies führt zu einer genaueren Approximation der PPF, und die generierten RVs folgen der exakten Verteilung genauer. Beachten Sie jedoch, dass dies auf Kosten einer teuren Einrichtung geht.

Die Einrichtungszeit hängt hauptsächlich davon ab, wie oft die PDF ausgewertet wird. Sie ist kostspieliger für PDFs, die schwer auszuwerten sind. Beachten Sie, dass wir die Normierungskonstante ignorieren können, um die Auswertungen der PDF zu beschleunigen. PDF-Auswertungen steigen während der Einrichtung für kleine Werte von u_resolution an.

>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def __init__(self):
...         self.callbacks = 0
...     def pdf(self, x):
...         self.callbacks += 1
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> # u_resolution = 10^-8
>>> # => less PDF evaluations required
>>> # => faster setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-8,
...                                  random_state=urng)
>>> dist.callbacks
4095        # may vary
>>> dist.callbacks = 0  # reset the number of callbacks
>>> # u_resolution = 10^-10 (default)
>>> # => more PDF evaluations required
>>> # => slow setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-10,
...                                  random_state=urng)
>>> dist.callbacks
11454       # may vary
>>> dist.callbacks = 0  # reset the number of callbacks
>>> # u_resolution = 10^-12
>>> # => lots of PDF evaluations required
>>> # => very slow setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
...                                  random_state=urng)
13902     # may vary

Wie wir sehen können, ist die Anzahl der erforderlichen PDF-Auswertungen sehr hoch, und eine schnelle PDF ist entscheidend für den Algorithmus. Dies hilft zwar, die Anzahl der Teilintervalle zu reduzieren, die zur Erreichung des Fehlerziels erforderlich sind, was Speicher spart und die Stichprobenziehung schnell macht. NumericalInverseHermite ist eine ähnliche Inversionsmethode, die die CDF basierend auf Hermite-Interpolation invertiert und über die maximale tolerierte Fehlersumme mittels u-Auflösung steuert. Sie verwendet jedoch viel mehr Intervalle als NumericalInversePolynomial

>>> from scipy.stats.sampling import NumericalInverseHermite
>>> # NumericalInverseHermite accepts a `tol` parameter to set the
>>> # u-resolution of the generator.
>>> rng_hermite = NumericalInverseHermite(norm(), tol=1e-12)
>>> rng_hermite.intervals
3000
>>> rng_poly = NumericalInversePolynomial(norm(), u_resolution=1e-12)
>>> rng_poly.intervals
252

Wenn die exakte CDF einer Verteilung verfügbar ist, kann man den vom Algorithmus erreichten u-Fehler schätzen, indem man die Methode u_error aufruft

>>> from scipy.special import ndtr
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...     def cdf(self, x):
...         return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)
>>> rng.u_error(sample_size=100_000)
UError(max_error=8.785949745515609e-11, mean_absolute_error=2.9307548109436816e-11)

u_error führt eine Monte-Carlo-Simulation mit einer gegebenen Anzahl von Stichproben durch, um den u-Fehler zu schätzen. Im obigen Beispiel werden 100.000 Stichproben von der Simulation verwendet, um den u-Fehler zu approximieren. Sie gibt den maximalen u-Fehler (max_error) und den mittleren absoluten u-Fehler (mean_absolute_error) in einem UError benannten Tupel zurück. Wie wir sehen, liegt der max_error unter der Standard-u_resolution (1e-10).

Es ist auch möglich, die PPF der gegebenen Verteilung auszuwerten, sobald der Generator initialisiert wurde

>>> rng.ppf(0.975)
1.959963985701268
>>> norm.ppf(0.975)
1.959963984540054

Dies können wir beispielsweise verwenden, um den maximalen und den mittleren absoluten u-Fehler zu überprüfen

>>> u = np.linspace(0.001, 0.999, num=1_000_000)
>>> u_errors = np.abs(u - dist.cdf(rng.ppf(u)))
>>> u_errors.max()
8.78600525666684e-11
>>> u_errors.mean()
2.9321444940323206e-11

Die vom Generator bereitgestellte approximierte PPF-Methode ist viel schneller auszuwerten als die exakte PPF der Verteilung.

Während der Einrichtung wird eine Tabelle von CDF-Punkten gespeichert, die zur Approximation der CDF verwendet werden kann, sobald der Generator erstellt wurde

>>> rng.cdf(1.959963984540054)
0.9750000000042454
>>> norm.cdf(1.959963984540054)
0.975

Wir können dies verwenden, um zu überprüfen, ob der Integrationsfehler bei der Berechnung der CDF \(0.05I_{0}\epsilon_{u}\) überschreitet. Hier ist \(I_0\) \(\sqrt{2\pi}\) (die Normierungskonstante für die Standardnormalverteilung)

>>> x = np.linspace(-10, 10, num=100_000)
>>> x_error = np.abs(dist.cdf(x) - rng.cdf(x))
>>> x_error.max()
4.506062190046123e-12
>>> I0 = np.sqrt(2*np.pi)
>>> max_integration_error = 0.05 * I0 * 1e-10
>>> x_error.max() <= max_integration_error
True

Die während der Einrichtung berechnete CDF-Tabelle wird zur Auswertung der CDF verwendet, und es ist nur eine weitere Feinabstimmung erforderlich. Dies reduziert die Aufrufe der PDF, aber da der Feinabstimmungsschritt die einfache Gauss-Lobatto-Quadratur verwendet, wird die PDF mehrmals aufgerufen, was die Berechnung verlangsamt.

Referenzen

[1]

Derflinger, Gerhard, Wolfgang Hörmann, und Josef Leydold. „Random variate generation by numerical inversion when only the density is known.“ ACM Transactions on Modeling and Computer Simulation (TOMACS) 20.4 (2010): 1-25.