Короткий ответ: ваша в порядке, но ваша неверна. Чтобы получить положительное стабильное распределение, заданное вашей формулой в R, вам нужно установить
γ γ = | 1 - я загар ( π α / 2 ) | - 1 / α .δγ
γ=|1−itan(πα/2)|−1/α.
Самый ранний пример формулы, которую я смог найти, был в (Feller, 1971), но я нашел эту книгу только в физической форме. Однако (Hougaard, 1986) дает ту же формулу вместе с преобразованием Лапласа
Из руководства ( используется в ) параметризация взята из (Самородницкий и Такк, 1994), другого ресурса, чье онлайн-воспроизведение ускользнуло от меня. Однако (Weron, 2001) дает характеристическую функцию в Самородницком и параметризацию Taqqu для чтобы быть
L(s)=E[exp(−sX)]=exp(−sα).
stabledist
stabledist
fBasics
pm=1
α≠1φ(t)=E[exp(itX)]=exp[iδt−γα|t|α(1−iβsign(t)tanπα2)].
Я переименовал некоторые параметры из статьи Верона в другую сторону с помощью обозначения, которое мы используем. Он использует для и для . В любом случае, подключив и , мы получим
μδσγβ=1δ=0φ(t)=exp[−γα|t|α(1−isign(t)tanπα2)].
Обратите внимание, что для и что . Формально, , поэтому, установив в мы получаем
Интересно отметить, что которая соответствует , также равна , поэтому, если вы попробуете или(1−itan(πα/2))/|1−itan(πα/2)|=exp(−iπα/2)α∈(0,1)iα=exp(iπα/2)γ = | 1 - я загар ( π α / 2 ) | - 1 / α φ ( t ) φ ( i s ) = exp ( - s α ) = L ( s ) . & gamma & alpha ; = 1 / 2 1 / 2 & gamma = & alpha ; & gamma = 1 - & alpha ; & alpha ;L(s)=φ(is)γ=|1−itan(πα/2)|−1/αφ(t)
φ(is)=exp(−sα)=L(s).
γα=1/21/2γ=αγ=1−α, что на самом деле не плохое приближение, вы в конечном итоге совершенно правильно для .
α=1/2
Вот пример в R для проверки правильности:
library(stabledist)
# Series representation of the density
PSf <- function(x, alpha, K) {
k <- 1:K
return(
-1 / (pi * x) * sum(
gamma(k * alpha + 1) / factorial(k) *
(-x ^ (-alpha)) ^ k * sin(alpha * k * pi)
)
)
}
# Derived expression for gamma
g <- function(a) {
iu <- complex(real=0, imaginary=1)
return(abs(1 - iu * tan(pi * a / 2)) ^ (-1 / a))
}
x=(1:100)/100
plot(0, xlim=c(0, 1), ylim=c(0, 2), pch='',
xlab='x', ylab='f(x)', main="Density Comparison")
legend('topright', legend=c('Series', 'gamma=g(alpha)'),
lty=c(1, 2), col=c('gray', 'black'),
lwd=c(5, 2))
text(x=c(0.1, 0.25, 0.7), y=c(1.4, 1.1, 0.7),
labels=c(expression(paste(alpha, " = 0.4")),
expression(paste(alpha, " = 0.5")),
expression(paste(alpha, " = 0.6"))))
for(a in seq(0.4, 0.6, by=0.1)) {
y <- vapply(x, PSf, FUN.VALUE=1, alpha=a, K=100)
lines(x, y, col="gray", lwd=5, lty=1)
lines(x, dstable(x, alpha=a, beta=1, gamma=g(a), delta=0, pm=1),
col="black", lwd=2, lty=2)
}
- Феллер В. (1971). Введение в теорию вероятностей и ее приложения , 2 , 2-е изд. Нью-Йорк: Уайли.
- Хугаард, П. (1986). Модели выживания гетерогенных популяций, полученных из стабильных распределений , Biometrika 73 , 387-396.
- Самородницкий Г., Такк М.С. (1994). Стабильные негауссовские случайные процессы , Чепмен и Холл, Нью-Йорк, 1994.
- Верон, Р. (2001). Пересмотрено стабильное распределение Леви: хвостовой индекс> 2 не исключает стабильный режим Леви , Международный журнал современной физики, 2001, 12 (2), 209-223.
Я думаю, что происходит, что в выходных данных
delta
может сообщаться значение внутреннего местоположения, в то время как во входных данныхdelta
описывается сдвиг. [Похоже, есть аналогичная проблема с тем,gamma
когдаpm=2
.] Так что, если вы попытаетесь увеличить сдвиг до 2затем вы добавляете 2 к значению местоположения.
С
beta=1
и уpm=1
вас есть положительная случайная величина с нижней границей распределения в 0.Сдвиг на 2 и нижняя граница возрастает на ту же величину
Но если вы хотите, чтобы
delta
входное значение было внутренним значением местоположения, а не сдвигом или нижней границей, то вам нужно использовать другую спецификацию для параметров. Например, если вы попробуете следующее (сpm=3
и с попыткой,delta=0
и с тем, чтоdelta=0.290617
вы нашли ранее), вы, похоже, получите то же самоеdelta
в и из. Сpm=3
иdelta=0.290617
вы получите ту же плотность 0.02700602, которую вы нашли ранее, и нижнюю границу в 0. С,pm=3
иdelta=0
вы получите отрицательную нижнюю границу (на самом деле -0.290617).Возможно, вам будет проще просто проигнорировать
delta
в выводе, и пока вы продолжаетеbeta=1
использоватьpm=1
средстваdelta
во входных данных, это нижняя граница распределения, которая, по-видимому, должна быть равна 0.источник
Также следует отметить, что Мартин Мачлер просто переработал код для стабильного дистрибутива и добавил некоторые улучшения.
Его новый пакет stabledist будет также использоваться fBasics, так что вы можете захотеть взглянуть и на это.
источник