更新

在我的问题中，泊松强度的估计取决于对单个事件的观察。这是一个重要因素，因为它不会破坏泊松通常的 MLE 并不明显。

我运行了一个模拟，其中绘制了随机强度，然后从这些强度中绘制了随机频率。我们只关注等于 18 年的频率。最后，我们使用等于观察频率的强度估计器，就像在泊松的通常 MLE 中一样，并与真实强度进行比较。

我怀疑对于在一个样本周期内精确观察到一次事件的罕见事件，泊松强度的估计可能比观察到的频率更好。然而，这个模拟显示了通常泊松 MLE 的零偏差。

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

ns = 100000
T = 10
lam = np.random.gamma(1,1/T,ns) # random intensities ~1/18 (per year)
obsf = np.random.poisson(lam*T,ns) # observe frequencies for 18 years

plt.hist(obsf)
plt.title('random intensities')
plt.show()

sidx = np.where(obsf==1)[0] # we only look at trials where 1 event happened in 18 years
err = (obsf[sidx]/T - lam[sidx]) # error when intensity estimator is set to observed frequency
print(stats.describe(err))

plt.hist(err)
plt.title('errors')
plt.show()

输出：

DescribeResult(nobs=25146, minmax=(-0.30917753984994045, 0.055486021558327742), mean=-0.00044121203771743959, variance=0.0015661042110224508, skewness=-1.434031238643681, kurtosis=2.993433694497429)

更新 2

我使用 Gamma 作为模拟参数分布。这是不公平的，因为伽马是泊松之前的。我将使用其他一些发行版。

2个回答

我会将其视为一次观察，假设您从 2000 年 12 月 31 日午夜开始，4 年，而另一次观察在 12.956 年被审查。在指数分布的情况下，这是一个很容易解决的问题：

p (x_{1}, x_{2} | λ) = (\frac{1}{λ} e^{- x_{1} / λ}) (e^{- c_{2} / λ})

$p(x_1, x_2 | \lambda) = \left(\frac{1}{\lambda}e^{-x_1/\lambda}\right)\left(e^{-c_2/\lambda}\right)$

其中第一项是观察的概率 $x_1$ 第二个是观察的概率 $x_2 > c_2$ ，和 $c_2$ 是 12.956 年的审查点。取对数并将导数设置为零给我们：

\hat{λ} = \frac{1}{x_{1} + c_{2}}

$\hat{\lambda} = \frac{1}{x_1 + c_2}$

这给了我们大致 $1/18$ 年，因为我们在 2017 年底而不是开始。

这只是通过一般结果的特定情况来处理，即删失指数分布的 MLE 等于观察到的总时间的总和，无论它们是否被删失，除以观察到的到达总数。

结论：这两种方法将给出相同的结果。

根据评论进行编辑：

基于泊松和指数的估计器将始终给出相同的答案。在泊松的情况下，你有 $n$ 在一段时间内观察到的事件 $T$ 给出一个估计 $T/n$ . 在指数的情况下，您观察到的总时间是 $T$ ，你有 $n+1$ 观察，最后一个在观察期结束时被审查，所以估计是 $n/T$ .

更仔细地观察你的模拟，我可以更清楚地看到你在做什么。您已经设置了一个数据生成过程，其中随机强度为每个伽马强度采用伽马分布，有一个泊松实现，然后您转向估计具有非零计数的子样本中的泊松率。

这与我对日期问题的设想有些不同，因为您有来自样本的数据，其中一些被故意排除在外，并且该样本具有重要的异质性。

这里有几件事要解决：如果您没有对泊松过程的随机输入，则观察到的数据的正确概率模型是零截断泊松分布。https://en.wikipedia.org/wiki/Zero-truncated_Poisson_distribution。请注意，强度的矩估计方法可以数值求解（没有易于处理的解析解）

0 = \bar{Y} - λ \exp (λ) / (\exp (λ) - 1)

$\begin{equation} 0 = \bar{Y} - \lambda \exp(\lambda) / (\exp(\lambda) - 1) \end{equation}$

set.seed(123)
ns <- 10000
zcpois <- function(y) uniroot(function(x) mean(y) - x*exp(x)/(exp(x)-1), c(1e-8, max(y)))
y <- rpois(ns, 0.7)
zcpois(y[y!=0]) ## very precise MoM

> zcpois(y[y!=0])$root
[1] 0.6950045

您还可以使用最大似然来回忆截断 RV 的可能性： $P(Y=y | Y>0) = P(Y=y)/P(Y>0)$ .

y0 <- y[y!=0]
negloglik <- function(lambda, y) {
  -sum(dpois(x=y, lambda=lambda, log=T) - ppois(q=0, lambda=lambda, lower.tail=F, log=T))
}

lseq <- seq(from=0.001, to=2, by=0.001)
nll <- sapply(lseq, negloglik, y=y0)
plot(lseq, nll, type='l')
i <- which.min(nll)
abline(v=lseq[i])
text(lseq[i], max(nll), paste('MLE:', lseq[i]), pos=4, xpd=T)

但是，重要的是要考虑泊松/伽玛混合以及对所得概率模型的影响。混合模型不遵循零截断泊松分布。潜在的泊松过程（没有遗漏 0 个计数）将遵循过度分散的泊松过程。

lambdas <- rgamma(ns, 1, 18) ## rate parametrization
counts <- rpois(ns, lambdas)
nzcounts <- counts[counts!=0]

> zcpois(nzcounts)$root
[1] 0.1311875
> mean(lambdas)
[1] 0.05573894

一个贝叶斯模型，其似然性遵循零截断泊松模型，其强度遵循具有非信息性形状和尺度参数的伽马分布，使您能够估计后验分布，其后验模式将为您提供对输入参数的良好估计伽马分布。一种常客方法会调用 EM 算法。

其它你可能感兴趣的问题

上一篇Softmax 层怎么不能使用硬目标收敛下一篇异常值过多时如何处理？