我在 PyMC3 中放置了一些代码来执行此任务，因为您在问题中提到了它。您似乎已经熟悉的第一部分是拟合模型以获得参数的后验分布：

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)

现在我们来看看你似乎在苦苦挣扎的部分。我们可以利用这个很好的性质：当两个随机变量和具有相同的过离散参数的负二项式分布时，那么也具有负二项式分布，均值为和相同的过度分散参数。您可以在此处找到此属性的证明。$X$$Y$$X+Y$$\mathbb E[X]+\mathbb E[Y]$$X$$Y$

假设负二项式参数是固定的（正式地，假设您的随机过程属于Lévy 过程的类，其中包括泊松过程），这意味着如果您想知道一小时内事件数的分布或者一整天，你只需要调整平均值，就像你对泊松过程所做的那样。

例如，要找出在一天内找到 200 多个事件的非典型情况，我们可以使用以下命令：

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

让我们稍微分解一下这行代码。当我们使用 时pm.NegativeBinomial.dist(mu=..., alpha=...)，我们使用一组特定的参数调用负二项式的 PyMC3 实现（我们也可以使用 Numpy 实现，但它们的参数化方式不同，因此坚持使用 PyMC3 不太容易出错）。

然后，我们使用从后验中采样的参数：alpha=trace["alpha"]过度离散和mu=48*trace["mean"]均值（我们乘以 48 以调整该均值以反映 24 小时而不是半小时）。

最后，我们从这个分布中采样许多实例，并将它们与我们感兴趣的值 ( .random(10**4)>200) 进行比较，然后从我们的负二项式过程中找到新样本超过它的概率（通过应用于np.mean结果布尔数组）。结果是您的模型每天生成 200 个或更多事件的概率。

这里有一些警告：

• 如果您的模型允许过度分散随时间变化，那么这些都不起作用
• 随时间变化，则必须对其进行一些调整。您不必将速率乘以某个数字，而是必须积分，这会使这变得更复杂一些。$\lambda(t)$$\lambda(t)$

编辑： 我正在编辑以解决@J 的评论，询问有关星期几的影响。因此，让我们首先生成一些具有强烈星期几影响的数据：

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

现在，我们可以解决这个问题的一种方法是使用 7 种不同的方法，而不是单一的方法。PyMC3 模型可以写成：

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

这里的变量将day每个观察结果与它来自的星期几相关联。现在，我们有一个模型可以考虑星期几的效果。我们如何检查星期五有超过 500 个事件是否非典型？该过程类似于同质情况：

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)

好的，现在如果我们要检查一周内的 3000 个事件是否是非典型事件怎么办？一周的预期事件数是48*sum(mean)，所以我们这样做：

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

请注意，我们不需要任何花哨的积分，因为这种星期几的效果使成为分段常数函数。（万岁！）。当函数形式有点复杂时，您将不需要积分泊松率：例如，如果是多项式、指数、从高斯过程中采样的函数等。不幸的是，它似乎是很难在网上找到关于这个特定主题的资源......也许我会在找到时间时在这个答案中添加一些解决这个问题的东西。$\lambda(t)$$\lambda(t)$

希望我有帮助！

负二项式可以像泊松一样处理，但是如何处理它是模棱两可的。这将取决于导致过度分散的潜在过程。这可以以不同的方式发生。

下面我将介绍两种方式：

• 负二项式以泊松分布和伽马分布的形式出现

  在这种情况下，成功概率参数会发生变化。$p$

• 负二项式作为计数过程发生，其中事件之间的间隔/等待时间是几何分布的。

  在这种情况下，参数会发生变化。$r$

---

1. 复合分布

您可以将负二项分布视为泊松分布与伽马分布的复合。

如果

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

那么

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

对于泊松过程，如果考虑更大的时间间隔，则事件数量的分布与具有更大速率系数的泊松分布变量有关。

例如，复合分布中的泊松率用因子缩放。$c$

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

这类似于缩放伽马分布的速率。

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

所以复合分布变为

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

---

2.计数过程

您可以将负二项分布视为发生在计数过程中，其中事件之间的等待时间是 iid 几何分布的。

如果您考虑事件的有序序列，其中事件之间的时间遵循几何分布：$1,2,...,k,k+1,...$

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

的区间内的事件数服从负二项分布，其中和$t$$r=\lfloor t \rfloor$$p=p$

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

在这种情况下，执行计数处理的时间段的增加对应于负二项分布$t$$r$

这种情况对应于 PedroSebe 的回答。

---

所以这将取决于你有什么样的过程会产生负二项式分布的计数。

这就是我在 R 中的做法。如果正确，它应该很容易翻译成 python。

首先估计最适合给定训练数据集的负二项分布的参数。然后使用这些参数将新数据映射到分布函数。

library(MASS)

set.seed(1234)
data_stream <- rnbinom(n= 1000, size= 1, mu= 10)

params <- fitdistr(x= data_stream, densfun= 'negative binomial', lower= c(1e-9, 0))
params
          size           mu     
   0.96289937   10.02900002 
 ( 0.04719405) ( 0.33835666)

new_time_point <- 30

pnbinom(new_time_point, size= params$estimate[1], mu= params$estimate[2])
0.94562 # This is how extreme the new data is