R函数loglik（）返回-inf？

Question

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

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

来自题为 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)

Answer 1

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

将测量错误片段替换为：

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

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

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

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

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

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

dmeas  <- function(reports,H,rho,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 病例的二项式概率H 种可能需要报告的感染，每种感染的报告概率为 rho。当二项式概率 Binom(N,p) 中的试验次数 N 为零时，唯一可能的结果为零 'successes'（报告的病例），概率为 1。所有其他结果的概率为 0（并且对数概率 -Inf).

R函数loglik（）返回-inf？

R function loglik() returning -inf?

simulation

modeling

r

function