R中涉及积分和贝塞尔函数的Newton Raphson代码答案

【问题标题】：Newton Raphson code in R involving integration and Bessel functionR中涉及积分和贝塞尔函数的Newton Raphson代码
【发布时间】：2014-10-10 13:14:33
【问题描述】：

我想估计涉及贝塞尔函数和积分的函数参数。但是，当我尝试运行它时，我收到一条消息“f(x, ...) 中的错误：找不到函数“BesselI””。我不知道如何修复它，如果有任何相关建议，我将不胜感激。

library(Bessel)
library(maxLik)
library(miscTools)


K<-300

f  <- function(theta,lambda,u) {exp(-u*theta)*Vectorize(BesselI(2*sqrt(t*u*theta*lambda),1))/u^0.5}   
F  <- function(theta,lambda){integrate(f,0,K,theta=theta,lambda=lambda)$value}
tt <- function(theta,lambda){(sqrt(lambda)*exp(-t*lambda)/(2*sqrt(t*theta)))*
                                   (theta*(2*t*lambda-1)*F(theta,lambda))}
loglik <- function(param) {
   theta <- param[1]
   lambda <- param[2]
   ll <-sum(log(tt(theta,lambda)))
}
t <- c(24,220,340,620,550,559,689,543)
res <- maxNR(loglik, start=c(0.001,0.0005),print.level=1,tol = 1e-08) 
summary(res)

牛顿-拉夫森最大化迭代次数：0 返回码：100 初始值超出范围。

我收到“有 50 个或更多警告（使用 warnings() 查看前 50 个）”，当我使用 warnings() 时，以下是警告。

In t * u : longer object length is not a multiple of shorter object length.

sessionInfo()
R version 2.14.2 (2012-02-29)
Platform: i386-pc-mingw32/i386 (32-bit)
locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252  
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                    
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] maxLik_1.1-2     miscTools_0.6-16   Bessel_0.5-4     Rmpfr_0.5-1  
[5] gmp_0.5-4       
loaded via a namespace (and not attached):
[1] sandwich_2.2-10

【问题讨论】：

Bessel 包真的安装了吗？
Yes Bessel 包已安装。
你能发布sessionInfo()的输出吗？为了清楚起见，我的意思是请将其添加到问题的末尾 - 而不是在 cmets 中。
BesselI 函数未矢量化...
即使在矢量化之后，我也有同样的错误，即“In t * u：较长的对象长度不是较短对象长度的倍数”

标签： r integration newtons-method bessel-functions

【解决方案1】：

警告：不完整的答案/探索。

让我们稍微简化一下，看看警告来自哪里，这可能会提供进一步的洞察力。至少我们可以消除这种可能性。

K <- 300 ## global variable, maybe a bad idea
t <- c(24,220,340,620,550,559,689,543)  ## global/same name as t(), ditto
library(Bessel)
f  <- function(theta,lambda,u) {
    exp(-u*theta)*BesselI(2*sqrt(t*u*theta*lambda),1)/u^0.5}
## Vectorize() isn't doing any good here anyway, take it out ...
F  <- function(theta,lambda){
       integrate(f,0,K,theta=theta,lambda=lambda)$value}

F() 有效，但仍然给出矢量化警告：

F(theta=1e-3,lambda=5e-4)
## [1] 3.406913
## There were 50 or more warnings (use warnings() to see the first 50)
f(theta=1e-3,lambda=5e-4,u=0:10)
##  [1]         NaN 0.010478182 0.013014566
##  [4] 0.017562239 0.016526010 0.016646497
##  [7] 0.018468755 0.016377872 0.003436664
## [10] 0.010399265 0.012919646
## Warning message:
## In t * u : longer object length is not a multiple of shorter object length

我们仍然收到警告。我们还可以看到，在 0 处评估被积函数很可能是一个问题（我们在分母中有一个 u 的平方根......）

看起来我们可能需要在外积（t 和 u 的所有组合）上评估 f()，然后可能对 @ 的值求和 f 987654329@?这确实需要解决，因为t 具有固定长度（可能是某种数据样本），而积分变量u 将具有任意多个值...

您能否提供对数似然函数的原始推导的链接？

【讨论】：

我没有这个派生的链接。 “tt”函数是密度函数，因为它涉及贝塞尔函数和积分，这就是为什么我只使用 sum(log(tt)) 因为可能性是 tt 的乘积，使用 log 乘积后将是 sum。
你需要重新安排一些事情（我现在没有时间）——你必须为 t 的单个值整合 u 的值，然后总结它们。

【解决方案2】：

这就是我所拥有的。不幸的是，我无法弄清楚如何在 t 向量上运行它，但我想你可以只做一个 for 循环。

library(base)
install.packages("maxLik")
library(maxLik)

#setting some initial values to test as I go
K <- 300
u <- 1
theta <- 1
T <- c(1,2,3,4)
lambda <- 1

#I got it to work with besselI instead of BesselI
#I also dropped Vectorize and instead do a for loop at the end
f  <- function ( theta, lambda, u ) {
    exp(-u * theta) * besselI(2 * sqrt(t * u * theta * lambda), 1) / u^0.5
}
#testing to see if everything is working
f(1,1,1)

#making a function with only one input so it can be integrated 
F <- function ( u ) {
    f(theta, lambda, u)
}
F(1)

#testing for integration
t <- T[1]
integrate(F, 0, K)

#entered the integration function inside here at the bottom
tt <- function ( theta, lambda ) {
    (sqrt(lambda) * exp(-t * lambda) / (2 * sqrt(t * theta))) *
        (theta * (2 * t * lambda - 1) * integrate(F, 0, K)$value)
}
tt(1,1)

#took the absolute value of theta and lambda just in case they turn out negative
loglik <- function(param) {
    theta <- param[1]
    lambda <- param[2]
    ll <-sum(log(tt(abs(theta),abs(lambda))))
    return(ll)
}

#example for using the for loop
for ( i in 1 : length(T) ) {
    t <- T[i]
    print(loglik(c(1,1)))
}

maxNR(loglik, start = c(1, 5), print.level = 1, tol = 1e-08)

我不确定您想要什么参数 theta 和 lamba，但我的设置有效并且使用了 t 的值。就奇异矩阵而言，t 和 theta/lambda 的大小可能存在一些依赖性。不确定 maxNR 是否使用广义逆，但我确定确实如此。

【讨论】：