R 函数 loglik() 返回 -inf？

问题描述

在 R 中模拟 SIR 模型。我有一个数据集，我正在尝试使用该模型准确绘图。我现在正在使用粒子过滤器功能，然后想对结果使用相应的logLik方法。当我这样做时，结果是“[1] -Inf”。我在文档中找不到为什么会这样以及如何避免它。我的模型参数不够准确吗？还有什么不对的吗？

我的函数如下所示： Sirsim %>% pfilter(Np=5000) -> pf logLik(pf)

来自名为 Likelihood for POMPS https://kingaa.github.io/sbied/pfilter/ 的在线课程，这是该课程的 R 脚本。但是，代码在这里有效......我不确定如何用它重现我的具体问题，不幸的是无法共享我正在使用的数据集或代码，因为它用于学术研究。

library(tidyverse)
library(pomp)
options(stringsAsFactors=FALSE)
stopifnot(packageVersion("pomp")>="3.0")
set.seed(1350254336)

library(tidyverse)
library(pomp)

sir_step <- Csnippet("
double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
double dN_IR = rbinom(I,1-exp(-mu_IR*dt));
S -= dN_SI;
I += dN_SI - dN_IR;
R += dN_IR;
H += dN_IR;
")

sir_init <- Csnippet("
S = nearbyint(eta*N);
I = 1;
R = nearbyint((1-eta)*N);
H = 0;
")

dmeas <- Csnippet("
lik = dbinom(reports,H,rho,give_log);
")

rmeas <- Csnippet("
reports = rbinom(H,rho);
")

read_csv("https://kingaa.github.io/sbied/pfilter/Measles_Consett_1948.csv") 
%>%
  select(week,reports=cases) %>%
  filter(week<=42) %>%
  pomp(
    times="week",t0=0,rprocess=euler(sir_step,delta.t=1/7),rinit=sir_init,rmeasure=rmeas,dmeasure=dmeas,accumvars="H",statenames=c("S","I","R","H"),paramnames=c("Beta","mu_IR","eta","rho","N"),params=c(Beta=15,mu_IR=0.5,rho=0.5,eta=0.06,N=38000)
  ) -> measSIR

measSIR %>%
  pfilter(Np=5000) -> pf
logLik(pf)

library(doParallel)
library(doRNG)
registerDoParallel()
registerDoRNG(652643293)
foreach (i=1:10,.combine=c) %dopar% {
  measSIR %>% pfilter(Np=5000)
} -> pf
logLik(pf) -> ll
logmeanexp(ll,se=TRUE)

解决方法

如果我在上面的代码中设置 Beta=100，我可以获得负无穷对数似然。

用这个替换测量误差片段：

dmeas <- Csnippet("
      double ll = dbinom(reports,H,rho,give_log);
      lik =  (!isfinite(ll) ? -1000 : ll );
")

似乎“解决”了问题，尽管您应该小心一点；像这样掩盖数字裂缝有时是可以的，但可以想象，以后可能会以某种方式回来咬你。如果您只需要足够长的时间避免非有限值以进入合理的参数范围，这可能没问题...

关于为什么会发生这种情况的一些猜测：

不知何故，当潜在的真实感染数量为零时，您会遇到“不可能”的情况，例如报告的病例数为正数。
有时，当非常小的正概率下溢到零时，会出现非有限对数似然。这里的等价物很可能是感染概率 1-exp(-Beta*I/N*dt) 达到 1.0；那么任何观察到的结果是，低于 100% 的人口被感染是不可能的。

您可以通过查看过滤后的轨迹实际是什么样子并将其与数据进行比较，或者通过在代码中添加调试语句来尝试诊断情况。如果有一种方法可以使用您的参数值运行确定性模拟，这可能会很快告诉您出了什么问题。

一种更简单/更直接的调试方法是用 R 函数替换您用于 dmeas 的 Csnippet：这会更慢但更容易使用（特别是如果您不熟悉C 编码）。如果您取消注释下面的 browser() 语句，当您遇到糟糕的情况时，代码将进入调试模式......

dmeas  <- function(reports,log,...) {
    lik <- dbinom(reports,size=H,prob=rho,log=log)
    if (!is.finite(lik)) {
        lik <- -1000
        ## browser()
    }
    return(lik)
}

例如：

(t = 3,reports = 2,S = 2280,I = 0,R = 35721,H = 0,Beta = 100,mu_IR = 0.5,rho = 0.5,eta = 0.06,N = 38000,log = TRUE)
Browse[1]> debug at /tmp/SO65554258.R!ZlSILG#7: return(lik)
Browse[2]> reports
[1] 2
Browse[2]> H
[1] 0
Browse[2]> rho
[1] 0.5

这表明问题确实在于，当感染为零时，您报告的病例数为正数……R 正在尝试计算观察 reports 个病例的二项式概率，当有 {{ 1}} 种可能可报告的感染，每种感染的报告概率为 H。当二项式概率 rho 中的试验次数 N 为零时，唯一可能的结果是零个“成功”（报告的案例），概率为 1。所有其他结果的概率为 0（和 log-概率-Inf).

function function modeling r r simulation