Dodatni rozkład stabilny w R.

9

Pozytywne rozkłady stabilne są opisane przez cztery parametry: parametr skośności , parametr skali , parametr lokalizacji i tak dalej - wywołany parametr indeksu . Gdy wynosi zero, rozkład jest symetryczny wokół , gdy jest dodatni (względnie ujemny), rozkład jest przekrzywiony w prawo (odpowiednio w lewo) Stabilne rozkłady pozwalają na grube ogony, gdy maleje.β[1,1]σ>0μ(,)α(0,2]βμα

Gdy jest ściśle mniejsze niż jeden, a obsługa dystrybucji ogranicza się do .αβ=1(μ,)

Funkcja gęstości ma tylko wyrażenie w formie zamkniętej dla niektórych szczególnych kombinacji wartości parametrów. Gdy , , , a to jest (patrz wzór (4.4) tutaj ):μ=0α<1β=1σ=α

f(y)=1πyk=1Γ(kα+1)k!(yα)ksin(αkπ)

Ma nieskończoną średnią i wariancję.

Pytanie

Chciałbym użyć tej gęstości w R. Używam

> alpha <- ...
> dstable(y, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)

gdzie funkcja dstable jest dostarczana z pakietem fBasics.

Czy możesz potwierdzić, że jest to właściwy sposób obliczenia tej gęstości w R?

Z góry dziękuję!

EDYTOWAĆ

Jednym z powodów, dla których jestem podejrzliwy, jest to, że na wyjściu wartość delta różni się od wartości na wejściu. Przykład:

> library(fBasics)
> alpha <- 0.4
> dstable(4, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)
[1] 0.02700602
attr(,"control")
   dist alpha beta gamma    delta pm
stable   0.4    1   0.4 0.290617  1
ocram
źródło

Odpowiedzi:

6

Krótka odpowiedź brzmi: jest w porządku, ale nie jest poprawna. Aby uzyskać dodatni rozkład stabilny podany przez twoją formułę w R, musisz ustawić δγ

γ=|1itan(πα/2)|1/α.

Najwcześniejszy przykład, jaki mogłem znaleźć podaną przez ciebie formułę, to (Feller, 1971), ale znalazłem tę książkę tylko w formie fizycznej. Jednak (Hougaard, 1986) podaje tę samą formułę wraz z transformacją Laplace'a Z podręcznika ( jest używany w ), parametryzacja pochodzi z (Samorodnitsky i Taqqu, 1994), innego zasobu, którego odtworzenie online mi umknęło. Jednak (Weron, 2001) daje charakterystyczną cechę Samorodnitsky and Taqqu za parametryzację za

L(s)=E[exp(sX)]=exp(sα).
stablediststabledistfBasicspm=1α1
φ(t)=E[exp(itX)]=exp[iδtγα|t|α(1iβsign(t)tanπα2)].
Zmieniłem nazwy niektórych parametrów z papieru Werona na monety przy użyciu notacji, której używamy. Używa dla i dla . W każdym razie, podłączając i , otrzymujemy μδσγβ=1δ=0
φ(t)=exp[γα|t|α(1isign(t)tanπα2)].

Zauważ, że dla i że . Formalnie , więc ustawiając w otrzymujemy Warto zwrócić uwagę na to, że które odpowiada , również wynosi , więc jeśli spróbujesz lub(1itan(πα/2))/|1itan(πα/2)|=exp(iπα/2)α(0,1)iα=exp(iπα/2)L(s)=φ(is)γ=|1itan(πα/2)|1/αφ(t)

φ(is)=exp(sα)=L(s).
γα=1/21/2γ=αγ=1α, co w rzeczywistości nie jest złym przybliżeniem, kończysz dokładnie na .α=1/2

Oto przykład w R, aby sprawdzić poprawność:

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)
}

Wyjście z wykresu

  1. Feller, W. (1971). Wprowadzenie do teorii prawdopodobieństwa i jej zastosowania , 2 , 2 wyd. Nowy Jork: Wiley.
  2. Hougaard, P. (1986). Modele przeżycia dla heterogenicznych populacji pochodzących ze stabilnych rozkładów , Biometrika 73 , 387-396.
  3. Samorodnitsky, G., Taqqu, MS (1994). Stable Non-Gaussian Random Process , Chapman & Hall, New York, 1994.
  4. Weron, R. (2001). Sprawdzone rozkłady stabilne podatkowo: indeks ogona> 2 nie wyklucza systemu stabilnego podatkowo , International Journal of Modern Physics C, 2001, 12 (2), 209-223.
P. Schnell
źródło
1
Cała przyjemność po mojej stronie. Temat pozytywnych stabilnych parametryzacji spowodował u mnie wiele bólu głowy na początku tego roku (to naprawdę bałagan), więc zamieszczam to, co wymyśliłem. Ta szczególna postać jest przydatna w analizie przeżycia, ponieważ forma Laplaciana pozwala na prosty związek między parametrami regresji warunkowej i brzeżnej w modelach proporcjonalnych zagrożeń, gdy istnieje kruchy termin po dodatnim rozkładzie stabilnym (patrz artykuł Hougaarda).
P Schnell,
6

Myślę, że dzieje się tak, że w danych wyjściowych deltamoże być zgłaszana wartość wewnętrznej lokalizacji, podczas gdy w danych wejściowych deltajest opisana zmiana. [Wydaje się, że jest podobny problem z gammakiedy pm=2.] Więc jeśli spróbujesz zwiększyć przesunięcie do 2

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1)
[1] 0.06569375
attr(,"control")
   dist alpha beta gamma    delta pm
 stable   0.4    1   0.4 2.290617  1

następnie dodajesz 2 do wartości lokalizacji.

Z beta=1i pm=1masz dodatnią zmienną losową z dolną granicą rozkładu na poziomie 0.

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=1))
[1] 0.002666507

Przesuń o 2, a dolna granica wzrośnie o tę samą wartość

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1))
[1] 2.003286

Ale jeśli chcesz, aby dane deltawejściowe były wewnętrzną wartością położenia, a nie przesunięciem lub dolną granicą, musisz użyć innej specyfikacji parametrów. Na przykład, jeśli spróbujesz wykonać następujące czynności (za pomocą pm=3i próbujesz delta=0i delta=0.290617znalazłeś wcześniej), wydaje się, że dostajesz to samo delta. Z pm=3i delta=0.290617otrzymujesz taką samą gęstość 0,02700602, którą znalazłeś wcześniej, a dolną granicę na 0. Z pm=3i delta=0otrzymujesz ujemną dolną granicę (w rzeczywistości -0,290617).

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3)
[1] 0.02464434
attr(,"control")
   dist alpha beta gamma delta pm
 stable   0.4    1   0.4     0  3
> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3)
[1] 0.02700602
attr(,"control")
   dist alpha beta gamma    delta pm
 stable   0.4    1   0.4 0.290617  3
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3))
[1] -0.2876658
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3))
[1] 0.004303485

Łatwiej może być po prostu zignorować deltadane wyjściowe i tak długo, jak długo będziesz beta=1używać pm=1środków deltawejściowych, dolna granica rozkładu, która wydaje się mieć wartość 0.

Henz
źródło
5

Warto również zauważyć: Martin Maechler właśnie dokonał zmiany kodu stabilnej dystrybucji i dodał kilka ulepszeń.

Jego nowy pakiet stabledist będzie również używany przez fBasics, więc możesz też chcieć to sprawdzić.

Dirk Eddelbuettel
źródło