Проблемы с использованием интеграции при записи функции правдоподобия

Question

Я пытаюсь написать функцию правдоподобия для распределения, которое определяется интеграцией. Я использую функцию integrate(), но когда я пытаюсь использовать это в остальной части функции, я получаю ошибку:

«Ошибка в B (alpha + i, beta + 6 - i)/B (alpha, beta): нечисловой аргумент двоичному оператору «

Значение интеграции, например,« 9.501501 с абсолютной ошибкой < 0.00078 ». Я пытался использовать trunc(), но это тоже не помогает. Я относительно новичок в R, так есть ли известное решение? Любая помощь будет оценена!

B <- function(a,b){ 
    integrand <- function(t){(t^(a-1))*((1-t)^(b-1))} 
    integrate(integrand,lower=0,upper=1) 
} 
betalik <- function(alpha,beta){ 
    likelihood <- 0 Z <- c(37,22,25,29,34,49) 
    for(i in 1:6) 
     likelihood <- likelihood + 
     Z[i]*log((B(alpha+i,beta+6-i))/B(alpha,beta)) 
    return(likelihood) 

}

Дориан,

Answer 1

Вслед за комментарии и ответы выше ...

Вот ваши первоначальные функции, отформатированных более красиво, интегрирующих @ комментарий Fhnuzoag в об извлечении $value компонента из integrate() результата:

B <- function(a,b){ 
    integrand <- function(t){(t^(a-1))*((1-t)^(b-1))} 
    integrate(integrand,lower=0,upper=1)$value 
} 
betalik <- function(alpha,beta){ 
    likelihood <- 0 
    Z <- c(37,22,25,29,34,49) 
    for(i in 1:6) likelihood <- likelihood + 
     Z[i]*log((B(alpha+i,beta+6-i))/B(alpha,beta)) 
    return(likelihood) 
} 

Здесь мы проверяем комментарий @ dason о том, что ваша функция B эквивалентна функции R(но функция R, безусловно, быстрее и, возможно, больше точный):

all.equal(B(1.1,2.7),beta(1.1,2.7)) ## TRUE 

Я предпочитаю, чтобы указать «данные» отдельно:

Z <- c(37,22,25,29,34,49) 

Новая версия функции правдоподобия, которая использует встроенный lbeta (лог-бета) функции, и векторизация:

blik2 <- function(alpha,beta,Z) { 
    index <- 1:6 
    sum(Z*(lbeta(alpha+index,beta+6-index)-lbeta(alpha,beta))) 
} 

all.equal(blik2(1.1,2.7,Z),betalik(1.1,2.7)) 
## small difference, blik2 is *probably* more 
## accurate ... "Mean relative difference: 6.406495e-08" 

(это может быть хорошо, чтобы обобщить это немного дальше и заменить значения 6 в коде с length(Z) ...)

Answer 2

Используя пример интеграции с ?integrate ....

integrate(dnorm, -1.96, 1.96) 
# 0.9500042 with absolute error < 1e-11 

Что на самом деле здесь происходит? Ну, интегральная функция создает объект «интегрировать» класс S3, который в основном представляет собой список с несколькими полями, а затем применяет к нему print(), выполняя в свою очередь print.integrate(). Выходной сигнал функции не фактически строка «0.9500042 с абсолютной ошибкой < 1e-11», она просто отображается таким образом.

Чтобы получить на то, что объект R, созданный integrate() на самом деле, сделать

obj = integrate(dnorm, -1.96, 1.96) 
str(obj) 
# List of 5 
# $ value  : num 0.95 
# $ abs.error : num 1.05e-11 
# $ subdivisions: int 1 
# $ message  : chr "OK" 
# $ call  : language integrate(f = dnorm, lower = -1.96, upper = 1.96) 
# - attr(*, "class")= chr "integrate" 

Так что, если вы просто хотите, значение интеграла, то вам придется извлечь value поле списка, созданного по этой функции, для ваших вычислений. Например.

10*integrate(dnorm, -1.96, 1.96)$value 
# [1] 9.500042