我已经想估计涉及贝塞尔函数和集成功能的参数。然而,当我试图运行它,我得到的信息是“错误的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级或更多的警告(使用警告()看到的第一个50)”,当我使用警告(),下面是警告。
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
警告:不完整的答案/探索。
让我们剥离下来一点点尝试,看看那里的警告是从,这可能会进一步洞察未来。至少我们可以消除的可能性。
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
的平方根在分母中...)
看起来我们可能需要在产品外,以评估f()
(t
和u
的所有组合),然后也许总结f
在t
的价值?这确实需要解决,因为t
具有固定长度(据推测这是某种形式的数据样本),而积分变量u
都将有任意多的值...
你能提供的链接数似然函数的原始推导?
下面是我得到了什么。不幸的是我无法弄清楚如何跨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和兰巴但我的设置工作,为t使用的值。有可能是对T和θ/λ的幅度有的依赖尽可能奇异矩阵去。不知道是否maxNR使用广义逆但我敢肯定它。