以项的无限总和形式编写的 PDF 参数的 MLE

Question

我的问题涉及当概率分布以无限和的形式表示时，使用 R 推导参数的最大似然估计，例如下面由 Rao、Girija 等人提出的问题。

我想看看当模型应用于给定数据集时，我是否可以重现这些作者（他们使用 Matlab，而不是 R）获得的最大似然估计。下面给出了我的尝试，尽管这会引发几个警告，即“较长的对象长度不是较短的对象长度的倍数”。 我知道为什么会收到此警告，但我不知道如何补救。 我如何编辑我的代码来解决这个问题？

此外，有没有更好的方法来处理无限和？在这里，我只是为 n (1000) 使用任意大数。

library(bbmle)
svec <- list(c=1,lambda=1)
x <- scan(textConnection("0.1396263 0.1570796 0.2268928 0.2268928 0.2443461 0.3141593 0.3839724 0.4712389 0.5235988 0.5934119 0.6632251 0.6632251 0.6981317 0.7679449 0.7853982 0.8203047 0.8377580 0.8377580 0.8377580 0.8377580 0.8726646 0.9250245 0.9773844 0.9948377 1.0122910 1.0122910 1.0646508 1.0995574 1.1170107 1.1170107 1.1170107 1.1344640 1.1344640 1.1868239 1.2217305 1.2740904 1.3613568 1.3613568 1.3613568 1.4486233 1.4486233 1.5358897 1.5358897 1.5358897 1.5707963 1.6057029 1.6057029 1.6231562 1.6580628 1.6755161 1.7104227 1.7453293 1.7976891 1.8500490 1.9722221 2.0594885 2.4085544 2.6703538 2.6703538 2.7052603 3.5604717 3.7524579 3.8920842 3.9444441 4.1364303 4.1538836 4.2411501 4.2586034 4.3633231 4.3807764 4.4854962 4.6774824 4.9741884 5.5676003 5.9864793 6.1086524"))
dL <- function(x, c,lambda,n = 1000,log=TRUE) {
    k <- 0:n
    r <- log(sum(lambda*c*(x+2*k*pi)^(-c-1)*(exp(-(x+2*k*pi)^(-c))^(lambda))))
    if (log) return(r) else return(exp(r))
}
dat <- data.frame(x)
m1 <- mle2( x ~ dL(c,lambda),
            data=dat,
            start=svec,
            control=list(parscale=unlist(svec)),
            method="L-BFGS-B",
            lower=c(0,0)
      )

Answer 1

我建议从该算法开始，并制作一个密度函数，可以通过对其运行ge 定义 (c(0, 2*pi) 进行积分来测试其行为是否正确。您正在调用它是一个“概率函数”，但我将其与 CDF 而不是密度分布 (PDF) 相关联：

dL <- function(x, c=1,lambda=1,n = 1000, log=FALSE) {
    k <- 0:n
    r <- sum(lambda*c*(x+2*k*pi)^(-c-1)*(exp(-(x+2*k*pi)^(-c))^(lambda)))
    if (log) {log(r) }
}
vdL <- Vectorize(dL)
integrate(vdL, 0,2*pi)
#0.999841 with absolute error < 9.3e-06
LL <- function(x, c, lambda){ -sum( log( vdL(x, c, lambda))) }

（我认为您试图在 log-likelihood 函数中包含太多内容，因此我决定拆分这些步骤。）

当我运行那个版本时，我从最后的 mle2 步骤收到一条我不喜欢的警告消息，我认为可能是这个密度函数偶尔返回负数值，所以这是我的最终版本：

 dL <- function(x, c=1,lambda=1,n = 1000) {
     k <- 0:n
     r <- max( sum(lambda*c*(x+2*k*pi)^(-c-1)*(exp(-(x+2*k*pi)^(-c))^(lambda))), 0.00000001)
     
 }
 vdL <- Vectorize(dL)
 integrate(vdL, 0,2*pi)
#0.999841 with absolute error < 9.3e-06
 LL <- function(x, c, lambda){ -sum( log( vdL(x, c, lambda))) }
 (m0 <- mle2(LL,start=list(c=0.2,lambda=1),data=list(x=x)))
#------------------------
Call:
mle2(minuslogl = LL, start = list(c = 0.2, lambda = 1), data = list(x = x))

Coefficients:
        c    lambda 
0.9009665 1.1372237 

Log-likelihood: -116.96

（警告和 warning-free LL 编号相同。）

所以我想我认为您试图在 log-likelihood 函数的定义中包含太多内容，结果在某个地方被绊倒了。应该有两个求和，一个用于密度近似，第二个用于 log-likelihood 的求和。这些总和中的数字会有所不同，因此您会看到错误。解包步骤至少在不抛出错误的情况下取得了成功。我不确定那个密度代表什么，也无法验证正确性。

至于是否有更好的方法来逼近无限级数的问题，答案取决于对部分和的收敛速度的了解，以及是否可以设置公差运行ce 值比较连续的值并在较少的项数后停止计算。

当我查看密度时，我想知道它是否适用于某些散射过程：

curve(vdL(x, c=.9, lambda=1.137), 0.00001, 2*pi)

您可以通过查看连续项的比率来检查收敛速度。这是一个函数，它对任意 x:

的前 10 个项执行此操作

> ratios <- function(x, c=1, lambda=1) {lambda*c*(x+2*(1:11)*pi)^(-c-1)*(exp(-(x+2*(1:10)*pi)^(-c))^(lambda))/lambda*c*(x+2*(0:10)*pi)^(-c-1)*(exp(-(x+2*(0:10)*pi)^(-c))^(lambda)) }
> ratios(0.5)
 [1] 1.015263e-02 1.017560e-04 1.376150e-05 3.712618e-06 1.392658e-06 6.351874e-07 3.299032e-07 1.880054e-07
 [9] 1.148694e-07 7.409595e-08 4.369854e-08
Warning message:
In lambda * c * (x + 2 * (1:11) * pi)^(-c - 1) * (exp(-(x + 2 *  :
  longer object length is not a multiple of shorter object length
> ratios(0.05)
 [1] 1.755301e-08 1.235632e-04 1.541082e-05 4.024074e-06 1.482741e-06 6.686497e-07 3.445688e-07 1.952358e-07
 [9] 1.187626e-07 7.634088e-08 4.443193e-08
Warning message:
In lambda * c * (x + 2 * (1:11) * pi)^(-c - 1) * (exp(-(x + 2 *  :
  longer object length is not a multiple of shorter object length
> ratios(0.5)
 [1] 1.015263e-02 1.017560e-04 1.376150e-05 3.712618e-06 1.392658e-06 6.351874e-07 3.299032e-07 1.880054e-07
 [9] 1.148694e-07 7.409595e-08 4.369854e-08
Warning message:
In lambda * c * (x + 2 * (1:11) * pi)^(-c - 1) * (exp(-(x + 2 *  :
  longer object length is not a multiple of shorter object length

对我来说，这看起来收敛得非常快，所以我猜您可以只使用前 20 个术语并获得类似的结果。有 20 个术语，结果如下：

> integrate(vdL, 0,2*pi)
0.9924498 with absolute error < 9.3e-06
> (m0 <- mle2(LL,start=list(c=0.2,lambda=1),data=list(x=x)))

Call:
mle2(minuslogl = LL, start = list(c = 0.2, lambda = 1), data = list(x = x))

Coefficients:
        c    lambda 
0.9542066 1.1098169 

Log-likelihood: -117.83

由于您从不尝试孤立地解释 LL，而是查看差异，因此我猜测微小的差异不会对您的推论产生不利影响。

以项的无限总和形式编写的 PDF 参数的 MLE

MLE of the parameters of a PDF written as an infinite sum of terms

r

mle