用于 mle2 公式调用的自定义密度函数定义错误

Question

我想定义我自己的密度函数，用于从R 的bbmle 包中对mle2 的公式调用。模型的参数是估计的，但我不能在返回的 mle2 对象上应用 residuals 或 predict 之类的函数。

这是我为简单泊松模型定义函数的示例。

library(bbmle)

set.seed(1)
hpoisson <- rpois(1000, 10)

myf <- function(x, lambda, log = FALSE) {
  pmf <- (lambda^x)*exp(-lambda)/factorial(x)
  if (log)
    log(pmf)
  else
    pmf
}

myfit <- mle2(hpoisson ~ myf(lambda), start = list(lambda=9), data=data.frame(hpoisson))
residuals(myfit)

在myfit 中，lambda 估计正确，但是当我在myfit 上调用残差时，我收到一条错误消息：

Error in myf(9.77598906811668) : 
  argument "lambda" is missing, with no default

另一方面，如果我简单地使用R 的内置dpois 函数按如下方式拟合模型，则会计算残差：

myfit <- mle2(hpoisson ~ dpois(lambda), start = list(lambda=9), data=data.frame(hpoisson))
    residuals(myfit)

谁能告诉我我在myf的函数定义中做错了什么？

谢谢

Answer 1

文档中解释的不是很清楚，但是使用自定义密度函数有几个先决条件：

• 函数的名称必须以d开头，必须有第一个参数x，并且必须有一个命名参数log。 （log 参数必须做一些明智的事情：特别是，mle2 将使用 log=TRUE 调用函数，并且函数最好返回对数似然！）一般来说，虽然这不是必需的，直接计算对数似然，然后在 log=FALSE 时取幂，而不是计算似然度并在 log=TRUE 时记录它（在某些情况下，例如零膨胀模型，这不是'不是真的可行）。例如，将我的 dmyf() 定义与 OP 代码中的 myf() 定义进行比较...
• 为了使用附加方法，例如predict，您必须定义一个名称以s 开头的附加函数；它返回指定参数的时刻列表、汇总统计信息等 - 请参见下面的示例，该示例复制自 bbmle::spois。

library("bbmle")
set.seed(1)
hpoisson <- rpois(1000, 10)

dmyf <- function(x, lambda, log = FALSE) {
    logpmf <- x*log(lambda)-lambda-lfactorial(x)
    if (log) return(logpmf)  else return(exp(logpmf))
}
smyf <- function(lambda) {
    list(title = "modified Poisson",
         lambda = lambda, mean = lambda,
         median = qpois(0.5, lambda),
         mode = NA, variance = lambda, sd = sqrt(lambda))
}
myfit <- mle2(hpoisson ~ dmyf(lambda),
              start = list(lambda=9), data=data.frame(hpoisson))
residuals(myfit)

Answer 2

不是真正的答案，但需要更多帮助：

我用它来尝试制作一个“自定义”beta-binomial 函数来模仿 bbmle 小插图的第一位中的那个。

set.seed(1001)
x1 <- rbetabinom(n=1000, prob=0.1, size=50, theta=10)
dmybetabinom <- function(x, N, theta, p, log=FALSE) {
    (choose(N,x)*beta(N-x+theta*(1-p),x+theta*p))/beta(theta*(1-p),theta*p)
}

该函数的工作原理类似于 dbetabinom：

dbetabinom(0:9,size=9,theta=4, prob=0.5) [1] 0.04545455 0.08181818 0.10909091 0.12727273 0.13636364 0.13636364 0.12727273 0.10909091 [9] 0.08181818 0.04545455 dmybetabinom(0:9,N=9,theta=4, p=0.5) [1] 0.04545455 0.08181818 0.10909091 0.12727273 0.13636364 0.13636364 0.12727273 0.10909091 [9] 0.08181818 0.04545455

但是当我尝试对其使用 mle2 功能时，我遇到了这个错误：

m0fa <- mle2(x1~dmybetabinom( N=50, theta, p), start=list(p=0.2, theta=9), data=data.frame(x1) )`

Error in optim(par = c(0.2, 9), fn = function (p)  :   non-finite finite-difference value [1] `