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

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

这是一个例子，我定义了一个简单泊松模型的函数。

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，拉姆达正确地估计，但是当我呼吁myfit残差，我得到它说的错误：

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

在另一方面，如果我只是拟合模型如下使用内置的dpois功能R的，残差计算：

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

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

感谢

来源

2015-02-07 Iliya Markov

这不是在文件中解释的很清楚，但也有使用自定义的密度函数的几个前提条件：

函数名必须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)

来源

2015-02-07 16:32:41

非常感谢。这样可行！ – 2015-02-07 17:24:24

亲爱的本，我的目标是定义一个非齐次泊松过程，其中dmyf中的lambda是协变量的函数。再次，我已经尝试过了，并且正确估计了协变量的系数。但是我想知道，在这种情况下，我在smyf中指定的时刻是否有意义，如果我能够正确预测过程的结果并对其进行一些诊断。你有什么建议，我可以找到更多信息？我已经尝试过NHPoisson产生很棒的情节，但我不确定我可以根据自己的需要量身定做。所以我想知道这是否可以用bbmle来完成。 – 2015-02-07 17:56:40

好吧，如果你想计算残差，那么你需要定义*一些*预测值，这样你就可以比较预测值和观测值... – 2015-02-07 18:03:37

不是一个真正的答案，但需要更多这方面的帮助：

我用这个来尝试做一个“自定义”β-二项功能模仿一个在第一位的小插曲。

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-15 19:20:35 emudrak

你可以使这个新的问题...我相信你的问题是你必须提供一个'log = TRUE'选项，因为这就是'mle2'要求的。我将在我的回答中详细说明这一点。 – 2016-02-15 20:46:48

谢谢！这解决了它。 – emudrak 2016-02-15 21:11:59

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

回答

相关问题