prob Пакет численно оценивает характеристические функции для базовых распределений R. Почти для всех распределений существуют существующие формулы. Однако для некоторых случаев решение в замкнутой форме неизвестно. Показательный пример: распределение Вейбулла (но см. Ниже).

Для характеристической функции Вейбулла я по существу вычисляю два интеграла и складываю их вместе:

fr <- function(x) cos(t * x) * dweibull(x, shape, scale)
fi <- function(x) sin(t * x) * dweibull(x, shape, scale)
Rp <- integrate(fr, lower = 0, upper = Inf)$value
Ip <- integrate(fi, lower = 0, upper = Inf)$value
Rp + (0+1i) * Ip

Да, это неуклюже, но работает на удивление хорошо! ...ahem, большую часть времени. Пользователь недавно сообщил, что следующие перерывы:

cfweibull(56, shape = 0.5, scale = 1)

Error in integrate(fr, lower = 0, upper = Inf) : 
  the integral is probably divergent

Теперь мы знаем, что интеграл не расходится, поэтому это должна быть численная задача. Немного поиграв, я мог заставить работать следующее:

fr <- function(x) cos(56 * x) * dweibull(x, 0.5, 1)

integrate(fr, lower = 0.00001, upper = Inf, subdivisions=1e7)$value
[1] 0.08024055

Это нормально, но это не совсем правильно, плюс требуется немало хлопот, которые плохо масштабируются. Я исследовал это для лучшего решения. Я нашел недавно опубликованную «закрытую форму» для характеристической функции сscale > 1 (посмотреть здесь), но это предполагаетОбобщенная гипергеометрическая функция Райта который не реализован в R (пока). Я заглянул в архивintegrate альтернативы, и есть масса вещей, которые, кажется, не очень хорошо организованы.

В рамках этого поиска мне пришло в голову перевести область интегрирования в конечный интервал через обратную касательную, иvoila! Проверьте это:

cfweibull3 <- function (t, shape, scale = 1){
  if (shape <= 0 || scale <= 0) 
    stop("shape and scale must be positive")
  fr <- function(x) cos(t * tan(x)) * dweibull(tan(x), shape, scale)/(cos(x))^2
  fi <- function(x) sin(t * tan(x)) * dweibull(tan(x), shape, scale)/(cos(x))^2
  Rp <- integrate(fr, lower = 0, upper = pi/2, stop.on.error = FALSE)$value
  Ip <- integrate(fi, lower = 0, upper = pi/2, stop.on.error = FALSE)$value
  Rp + (0+1i) * Ip
}

> cfweibull3(56, shape=0.5, scale = 1)
[1] 0.08297194+0.07528834i

Questions: Can you do better than this? Is there something about numerical integration routines that people who are expert about such things could shed some light on what's happening here? I have a sneaking suspicion that for large t the cosine fluctuates rapidly which causes problems...? Are there existing R routines/packages which are better suited for this type of problem, and could somebody point me to a well-placed position (on the mountain) to start the climb? Comments: Yes, it is bad practice to use t as a function argument. I calculated the exact answer for shape > 1 using the published result with Maple, and the brute-force-integrate-by-the-definition-with-R kicked Maple's ass. That is, I get the same answer (up to numerical precision) in a small fraction of a second and an even smaller fraction of the price. Edit:

Я собирался записать точные интегралы, которые я искал, но кажется, что этот конкретный сайт не поддерживает MathJAX, поэтому я вместо этого дам ссылки. Я рассчитываю на численную оценкухарактеристическая функция изРаспределение Вейбулла для разумных затратt (что бы это ни значило). Значение представляет собой комплексное число, но мы можем разделить его на его действительные и мнимые части, и это то, что я назвалRp а такжеIp выше.

Последний комментарий: в Википедии есть формула (бесконечный ряд) для Вейбулла c.f. и эта формула соответствует доказанной в статье, на которую я ссылался выше,howeverЭта серия доказана только дляshape > 1, Дело0 < shape < 1 остается открытой проблемой; см. бумагу для деталей.