mle2公式调用的自定义密度函数定义错误 [英] Error with custom density function definition for mle2 formula call

问题描述

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

I want to define my own density function to be used in the formula call to mle2 from R's bbmle package. The parameters of the model are estimated but I cannot apply functions like residuals or predict on the returned mle2 object.

这是一个示例，其中我为简单的泊松模型定义了一个函数.

This is an example in which I define a function for a simple Poisson model.

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上调用残差时，出现错误，提示:

In myfit, the lambda is estimated correctly but when I call residuals on myfit, I get an error which says:

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

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

On the other hand, if I simply fit the model as follows using R's built-in dpois function, the residuals are computed:

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

任何人都可以告诉我myf的函数定义中我做错了什么吗?

Could anyone please tell me what I am doing wrong in the function definition of myf?

谢谢

推荐答案

文档中没有很清楚地说明，但是使用自定义密度函数有一些先决条件:

It's not very clearly explained in the documents, but there are a few prerequisites for using custom density functions:

• 函数名称必须以d 开头，必须具有第一个自变量x，并且必须具有命名自变量log. ( log参数必须做一些明智的事情:尤其是mle2将使用log=TRUE调用该函数，并且该函数最好返回对数似然！)通常，尽管它不是必需的，直接计算对数似然然后在log=FALSE处取幂，而不是计算可能性并在log=TRUE处记录它，在数字上更明智(在某些情况下，例如零膨胀模型，真的不可行).例如，将我的dmyf()定义与OP代码中的myf()定义进行比较...
• 为了使用诸如predict之类的其他方法，您必须定义一个名称以s开头的附加函数.它返回指定参数的弯矩列表，摘要统计信息等-参见下面的示例，该示例是从bbmle::spois复制而来的.

• the function's name must start with d, must have first argument x, and must have a named argument log. (The log argument has to do something sensible: in particular, mle2 will call the function with log=TRUE, and the function had better return the log-likelihood!) In general, although it's not required, it's more numerically sensible to compute the log-likelihood directly and then exponentiate if log=FALSE, rather than computing the likelihood and logging it if log=TRUE (there are cases, such as zero-inflated models, where this isn't really feasible). For example, compare my dmyf() definition with the myf() definition in the OP's code ...
• in order to use additional methods such as predict you have to define an additional function whose name starts with s; it returns a list of moments, summary statistics, etc. for a specified parameter -- see example below, which is copied from 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)

这篇关于mle2公式调用的自定义密度函数定义错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！