Пакет 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
Да, это неуклюже, но это работает на удивление хорошо!... гм, большую часть времени. Недавно пользователь сообщил, что следующие разрывы:
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, и там есть масса вещей, которые не кажутся очень хорошо организованными.
В рамках этого поиска мне пришло в голову преобразовать область интегрирования в конечный интервал через обратную касательную и вуаля! Проверьте это:
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
Вопросы:
- Вы можете сделать лучше, чем это?
- Есть ли что-то в процедурах цифровой интеграции, что люди, которые разбираются в таких вещах, могут пролить свет на то, что происходит здесь? У меня есть тайное подозрение, что при больших
tкосинус быстро флуктуирует, что вызывает проблемы...? - Существуют ли существующие подпрограммы/пакеты R, которые лучше подходят для такого типа проблем, и может ли кто-нибудь указать мне на хорошо расположенную позицию (на горе), чтобы начать подъем?
Комментарии:
- Да, использование
tв качестве аргумента функции - это плохая практика. - Я вычислил точный ответ для
shape > 1, используя опубликованный результат с помощью Maple, иbrute-force-integrate-by-the-definition-with-Rногами клена. То есть, я получаю один и тот же ответ (вплоть до числовой точности) за небольшую долю секунды и даже меньшую часть цены.
Изменить:
Я собирался записать точные интегралы, которые я ищу, но кажется, что этот конкретный сайт не поддерживает MathJAX, поэтому я дам ссылки. Я хочу, чтобы численно оценить характеристическую функцию распределения Weibull для разумных входов t (что угодно это значит). Значение представляет собой комплексное число, но мы можем разбить его на его действительную и мнимую части и на то, что я называл Rp и Ip выше.
Один заключительный комментарий: Wikipedia имеет формулу, указанную (бесконечную серию) для Weibull c.f. и эта формула соответствует той, которая доказана в документе, на который я ссылался выше, однако, эта серия доказана для shape > 1. Случай 0 < shape < 1 остается открытой проблемой; Подробнее см. в документе.
