2

我想定义我自己的密度函数，以在对mle2from包R的公式调用中使用。bbmle模型的参数是估计的，但我不能在返回的对象上应用residuals或应用函数。predictmle2

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

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什么？

谢谢

4

2 回答 2

3

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

  • 函数的名称必须以 开头d，必须有第一个参数x，并且必须有一个命名参数log。（参数必须做一些有意义的log事情：特别是，mle2将调用函数log=TRUE，并且函数最好返回对数似然！）一般来说，虽然不是必需的，但直接计算对数似然在数值上更合理，并且然后取幂 if log=FALSE，而不是计算可能性并记录它 if log=TRUE（在某些情况下，例如零膨胀模型，这实际上并不可行）。例如，将我dmyf()的定义与myf()OP 代码中的定义进行比较...
  • 为了使用其他方法，例如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)
于 2015-02-07T16:32:41.930 回答
0

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

我用它来尝试制作一个“自定义”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] `
于 2016-02-15T19:20:35.100 回答