Интеграция в R не работает, когда верхний предел равен Inf

Мне нужно численно аппроксимировать дисперсию логарифма суммы двух логарифмически нормальных случайных величин. Я хотел бы сделать это в R, вот пример:

################
## Setup data ##
################

sdX1 = 0.33
sdX2 = 0.70
muX1 = log(32765) - 0.5 * sdX1^2
muX2 = log(52650) - 0.5 * sdX2^2

####################################
## PDF for sum of 2 lognormal RVs ##
####################################

d2lnorm = function(z, muX1, muX2, sdX1, sdX2){

    #PDFs
    L1 = distr::Lnorm(meanlog = muX1, sdlog = sdX1)
    L2 = distr::Lnorm(meanlog = muX2, sdlog = sdX2)

    #Convlution integral
    L1plusL2 = distr::convpow(L1 + L2, 1)

    #Density function
    f.Z = distr::d(L1plusL2)

    #Evaluate
    return(f.Z(z))

}

############################################
## Expectation for sum of 2 lognormal RVs ##
############################################    

ex2lnorm = function(muX1, muX2, sdX1, sdX2){    

    #E(g(x)) = integral of g(x) * f(x) w.r.t x
    integrate(function(z) log(z) * d2lnorm(z, muX1 = muX1, muX2 = muX2, sdX1 = sdX1, sdX2 = sdX2), lower = 0, upper = +Inf)$value

}

##############
## Run code ##
##############

ex2lnorm(muX1, muX2, sdX1, sdX2)

Однако интеграл оценивается как 0. Если я изменю верхний предел integrate на большое число, это сработает. Тем не менее, я запускаю кучу симуляций и не могу возиться с верхним пределом, чтобы заставить его работать каждый раз. Есть ли другой способ заставить это работать последовательно?


r integrate
person user13317    schedule 25.07.2017    source источник
comment
Мне кажется, что значения критически зависят от значений верхнего. Попробуйте изменить верхнее значение на 10 ^ (3:10)   -  person IRTFM    schedule 26.07.2017

Ответы (1)


arrow_upward
1
arrow_downward

Что вы могли бы сделать, если бы нашли максимальное верхнее значение, которое не дает 0:

ex2lnorm = function(muX1, muX2, sdX1, sdX2){    

  f <- function(z) log(z) * d2lnorm(z, muX1 = muX1, muX2 = muX2, sdX1 = sdX1, sdX2 = sdX2)

  upper <- 2^10
  cond.lower <- FALSE
  repeat {
    new <- integrate(f, lower = 0, upper = upper)$value
    if (new > 0) cond.lower <- TRUE
    if (cond.lower && new == 0) break
    upper <- upper * 2
    prev <- new
  }

  prev
}
person F. Privé    schedule 25.07.2017
comment
Я надеялся, что будет волшебная пуля, но поиск кажется хорошим решением. - person user13317; 26.07.2017