R: проблемы сходимости с числовым интегрированием

Question

Не уверен, что эта проблема с числовыми методами действительно должна быть здесь или в перекрестной проверке, но, поскольку у меня есть хороший воспроизводимый пример, я бы начал здесь.

Я собираюсь оценивать и подбирать группу распределений как для некоторых больших наборов данных, так и для наборов данных, генерируемых случайным образом из аналогичных распределений. В рамках этого процесса я буду генерировать оценки условного среднего для различных диапазонов значений, в том числе усеченных и неусеченных значений правого хвоста.

Функция cr_moment ниже, учитывая функцию PDF для dfun и параметры для этой функции в params вычисляет безусловное среднее значение этого распределения. Учитывая верхнюю, нижнюю или обе границы, он вычисляет условное среднее для диапазона, заданного этими границами, используя одно- или двукратно усеченное распределение для этих границ. Функция под ним, cr_gb2, специализирует cr_moment на обобщенном бета-распределении второго рода. И, наконец, значения параметров, приведенные ниже, которые аппроксимируют текущее долларовое распределение доходов домохозяйств в США по данным переписи населения США /BLS за 2000 год. McDonald & Ransom 2008. (Кроме того, спасибо Микко Мартилле в этом списке за помощь в кодировании эта функция).

Эта функция дает мне ошибку сближения, скопированную ниже, для различных нижних границ и верхней границы, равной 4.55e8 или выше, но не в 4.54e8. K-й момент GB2 существует для k

Error in integrate(f = prob_interval, lower = lb, upper = ub, subdivisions = 100L):
   the integral is probably divergent

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

cr_moment <- function(lb = -Inf, ub = Inf, dfun, params, v=1, ...){
  x_pdf         <- function(X){
    X^v * do.call(what=dfun, args=c(list(x=X), params))
    }
  prob_interval <- function(X){
    do.call(what=dfun, args=c(list(x=X), params))
    }
  integral_val  <- integrate(f = x_pdf, lower = lb, upper = ub)
  integral_prob <- integrate(f = prob_interval, lower = lb, upper = ub)
  crm <-  interval_val[[1]] / interval_prob[[1]]
  out <- list(value = integral_val[[1]], probability = integral_prob[[1]],  
              cond_moment = crm)
  out
}

library(GB2)
cr_gb2 <- function(lb = -Inf, ub = Inf, v = 1, params){
  cr_moment(lb, ub, dfun = dgb2, params = get("params"))
}

GB2_params <- list(shape1 = 2.2474, scale = 58441.5, shape2 = 0.6186, shape3 = 1.118)

cr_gb2(lb=1, ub= 4.55e8, params = GB2_params)