Können negative Binomialparameter wie Poisson behandelt werden?

Sep 14 2020

Ich habe einen Zählprozess, den ich mit einem Poisson-Prozess modellieren möchte. Die Daten werden alle 30 Minuten gemessen, und mit einer Poisson-Verteilung kann ich leicht die Wahrscheinlichkeit messen, dass eine bestimmte Anzahl von Ereignissen in verschiedenen Zeiträumen anomal ist, indem ich einen angepassten Lambda-Wert verwende, dh "ist die Anzahl der Ereignisse, die wir in der 30 min anomal? Was ist mit der letzten Stunde? Ist die Anzahl der Ereignisse, die wir in den letzten 1,5 Stunden gesehen haben, anomal? "Usw.

Das Problem ist, dass meine Daten überstreut sind und definitiv durch eine negative Binomialverteilung gut beschrieben werden. Ich wähle die Parameter zu verwenden$(\mu, \alpha)$ denn genau das nutzt PyMC3, wo $\mu$ entspricht Lambda aus der Poisson-Verteilung.

Gibt es eine Möglichkeit, die negativen Binomialparameter auf dieselbe Weise wie den Poisson-Ratenparameter zu verwenden, bei dem ich sehe, ob eine Ereignisanzahl in einem bestimmten Zeitraum t anomal ist (wobei ich t auf verschiedene Zeiträume erweitern kann)?

Antworten

7 PedroSebe Sep 19 2020 at 20:38

Ich habe Code für diese Aufgabe in PyMC3 erstellt, da Sie ihn in der Frage erwähnt haben. Der erste Teil, mit dem Sie anscheinend bereits vertraut sind, besteht darin, das Modell anzupassen, um eine posteriore Verteilung der Parameter zu erhalten:

import pymc3 as pm
import numpy as np

# generating simulated data data for a week
data = pm.NegativeBinomial.dist(mu=3, alpha=1).random(size=7*24*2)

# defining the model and sampling (MCMC)
with pm.Model() as model:
    alpha = pm.Exponential("alpha", 2.0)
    mean = pm.Exponential("mean", 0.2)
    obs_data = pm.NegativeBinomial("obs_data", mu=mean, alpha=alpha, observed=data)
    trace = pm.sample()

# plotting the posterior
pm.traceplot(trace)
pm.plot_posterior(trace)

Jetzt kommen wir zu dem Teil, mit dem Sie zu kämpfen scheinen. Wir können diese nette Eigenschaft verwenden: Wenn zwei Zufallsvariablen,$X$ und $Y$ haben dann negative Binomialverteilungen mit dem gleichen Überdispersionsparameter $X+Y$ hat auch eine negative Binomialverteilung mit Mittelwert $\mathbb E[X]+\mathbb E[Y]$ und der gleiche Überdispersionsparameter wie $X$ und $Y$. Beweise für diese Eigenschaft finden Sie hier .

Angenommen, die negativen Binomialparameter sind fest (formal vorausgesetzt, Ihr stochastischer Prozess gehört zur Klasse der Lévy-Prozesse , in denen Poisson-Prozesse enthalten sind), bedeutet dies, dass Sie die Verteilung für die Anzahl der Ereignisse in einer ganzen Stunde kennen möchten oder einen ganzen Tag müssen Sie nur den Mittelwert anpassen, wie Sie es bei einem Poisson-Prozess tun würden.

Um beispielsweise herauszufinden, wie untypisch es wäre, mehr als 200 Ereignisse an einem Tag zu finden, könnten wir Folgendes verwenden:

np.mean(pm.NegativeBinomial.dist(mu=48*trace["mean"], alpha=trace["alpha"]).random(10**4)>200)

Lassen Sie uns diese Codezeile ein wenig aufschlüsseln. Wenn wir verwenden pm.NegativeBinomial.dist(mu=..., alpha=...), rufen wir die PyMC3-Implementierung des negativen Binomials mit einem bestimmten Satz von Parametern auf (wir könnten auch die Numpy-Implementierung verwenden, aber sie sind unterschiedlich parametrisiert, sodass es weniger fehleranfällig ist, sich an PyMC3 zu halten).

Wir verwenden dann die Parameter, die wir vom Seitenzahn abgetastet haben: alpha=trace["alpha"]für die Überdispersion und mu=48*trace["mean"]für den Mittelwert (wir multiplizieren mit 48, um diesen Mittelwert so anzupassen, dass er 24 Stunden statt einer halben Stunde widerspiegelt).

Schließlich werden viele Instanzen aus dieser Verteilung abgetastet und mit dem Wert verglichen, an dem wir interessiert sind ( .random(10**4)>200). Anschließend wird die Wahrscheinlichkeit ermittelt, dass neue Stichproben aus unserem negativen Binomialprozess diese überschreiten (indem wir sie np.meanauf das resultierende Array von Booleschen Werten anwenden). Das Ergebnis ist die Wahrscheinlichkeit, dass Ihr Modell einen Tag mit 200 oder mehr Ereignissen generiert.

Ein paar Einschränkungen hier:

Wenn Ihr Modell zulässt, dass sich die Überdispersion mit der Zeit ändert, funktioniert nichts davon
Wenn Ihr Modell zulässt, dass sich die Poisson-Rate als Funktion mit der Zeit ändert $\lambda(t)$muss dies etwas angepasst werden. Anstatt die Rate mit einer Zahl zu multiplizieren, müssen Sie integrieren$\lambda(t)$, was dies etwas komplizierter macht.

BEARBEITEN: Ich bearbeite, um den Kommentar von @J zu adressieren. Fragen Sie nach Wochentagseffekten. Lassen Sie uns zunächst einige Daten mit starken Wochentagseffekten generieren:

# how many weeks of data are available?
WEEKS = 5
# how many observations are available per day?
OBS_PER_DAY = 24*2

data = pm.NegativeBinomial.dist(mu=[2,3,1,2,5,9,7]*5, alpha=1).random(size=OBS_PER_DAY).T.flatten()

Eine Möglichkeit, dies zu umgehen, besteht darin, 7 verschiedene Mittel anstelle eines einzigen zu verwenden. Das PyMC3-Modell kann wie folgt geschrieben werden:

with pm.Model() as model:
    alpha = pm.Exponential("alpha", 2.0)
    mean = pm.Exponential("mean", 0.2, shape=7)
    day = np.arange(WEEKS*7*OBS_PER_DAY)//OBS_PER_DAY%7
    obs_data = pm.NegativeBinomial("obs_data", mu=mean[day], alpha=alpha,
        observed=data)
    trace = pm.sample()

Die Variable dayhier ordnet jede Beobachtung dem Wochentag zu, von dem sie stammt. Jetzt haben wir ein Modell, das Wochentagseffekte berücksichtigt. Wie können wir überprüfen, ob mehr als 500 Veranstaltungen an einem Freitag untypisch sind? Das Verfahren ähnelt dem homogenen Fall:

friday = 4 # assuming the week starts on monday
np.mean(pm.NegativeBinomial.dist(mu=48*trace["mean"][:,friday], alpha=trace["alpha"]).random(10**4)>500)

OK, was ist, wenn wir überprüfen möchten, ob 3000 Ereignisse in einer Woche ein atypisches Ereignis sind? Die erwartete Anzahl von Ereignissen für eine Woche ist 48*sum(mean), also machen wir das:

np.mean(pm.NegativeBinomial.dist(mu=48*trace["mean"].sum(axis=1), alpha=trace["alpha"]).random(10**4)>3000)

Beachten Sie, dass wir keine ausgefallene Integration brauchten, da dieser Wochentagseffekt bewirkt $\lambda(t)$eine stückweise konstante Funktion . (Hurra!). Sie müssen die Poisson-Rate nicht integrieren, wenn ihre funktionale Form etwas komplizierter ist: zum Beispiel, wenn$\lambda(t)$ ist ein Polynom, ein Exponential, eine Funktion, die aus einem Gaußschen Prozess usw. entnommen wurde. Leider scheint es schwierig zu sein, Ressourcen zu diesem bestimmten Thema im Web zu finden ... Vielleicht werde ich dieser Antwort etwas hinzufügen, das dieses Problem behebt, wenn ich finde die Zeit.

Hoffe ich war hilfreich!

3 SextusEmpiricus Sep 19 2020 at 18:47

Das negative Binom kann wie Poisson behandelt werden, es ist jedoch nicht eindeutig, wie es zu behandeln ist. Dies hängt von dem zugrunde liegenden Prozess ab, der die Überdispersion verursacht. Dies kann auf verschiedene Arten geschehen.

Im Folgenden werde ich zwei Möglichkeiten beschreiben:

Das negative Binomial tritt als Poisson-Verteilung auf, die mit einer Gamma-Verteilung verbunden ist

In diesem Fall ist der Erfolgswahrscheinlichkeitsparameter $p$ Änderungen.
Das negative Binomial tritt als Zählvorgang auf, bei dem das Intervall / die Wartezeit zwischen Ereignissen geometrisch verteilt ist .

In diesem Fall ist die $r$ Parameteränderungen.

1. Zusammengesetzte Verteilung

Sie können die negative Binomialverteilung als Poisson-Verteilung anzeigen, die mit einer Gammaverteilung zusammengesetzt ist.

Wenn

$$Y \sim Poisson(\lambda=X)$$ wo $$X \sim Gamma(\alpha,\beta)$$

Dann $$Y \sim NB(r=\alpha, p = (\beta+1)^{-1})$$

Wenn Sie bei einem Poisson-Prozess ein größeres Zeitintervall berücksichtigen, bezieht sich die Verteilung der Anzahl der Ereignisse auf eine verteilte Poisson-Variable mit einem größeren Ratenkoeffizienten.

Beispielsweise wird die Poisson-Rate in der Verbundverteilung mit einem Faktor skaliert $c$.

$$Y_c \sim Poisson(\lambda=cX)$$

Dies ähnelt der Skalierung der Rate der Gammaverteilung.

$$cX \sim Gamma(\alpha,\beta/c)$$

So wird die zusammengesetzte Verteilung

$$Y_c \sim NB(r=\alpha, p = (\beta/c+1)^{-1})$$

2. Zählvorgang

Sie können die negative Binomialverteilung als in einem Zählvorgang auftretend betrachten, bei dem die Wartezeit zwischen Ereignissen in geometrischer Verteilung erfolgt.

Wenn Sie die geordnete Abfolge von Ereignissen berücksichtigen $1,2,...,k,k+1,...$ wobei die Zeit zwischen Ereignissen einer geometrischen Verteilung folgt:

$$t_k-t_{k-1} \sim Geom(p)$$

Dann die Anzahl der Ereignisse innerhalb eines Intervalls der Länge $t$ folgt einer negativen Binomialverteilung mit $r=\lfloor t \rfloor$ und $p=p$

$$N_{\text{events within $t$}} \sim NB(\lfloor t \rfloor, p)$$

In diesem Fall die Verlängerung des Zeitraums $t$ über die der Zählvorgang durchgeführt wird, entspricht einer Erhöhung des Parameters $r$ in der negativen Binomialverteilung.

Dieser Fall entspricht der Antwort von PedroSebe.

Es hängt also davon ab, welche Art von Prozess die negative Binomialverteilung der Zählungen erzeugt.

1 dariober Sep 14 2020 at 13:59

So würde ich es in R machen. Wenn es richtig ist, sollte es einfach sein, es in Python zu übersetzen.

Schätzen Sie zunächst die Parameter der negativen Binomialverteilung, die am besten zu einem bestimmten Trainingsdatensatz passen. Ordnen Sie dann die neuen Daten der Verteilungsfunktion mit diesen Parametern zu.