Код Ньютона-Рафсона в R, включающий интеграцию и функцию Бесселя

Предупреждение: неполный ответ / исследование.

Давайте немного сократим это, чтобы попытаться понять, откуда приходят предупреждения, что может дать дополнительное понимание. По крайней мере, мы можем исключить такую ​​возможность.

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 и lamba, но мои настройки работали и использовались значения t. Могут быть некоторые зависимости от величины t и тета / лямбда в отношении сингулярных матриц. Не уверен, что maxNR использует обобщенные инверсии, но я уверен, что это так.