Distribusi stabil positif di R

9

Distribusi stabil positif dijelaskan oleh empat parameter: parameter skewness , parameter skala , parameter lokasi , dan sebagainya -called index parameter . Ketika adalah nol distribusi simetris di sekitar , ketika itu positif (resp. negatif) distribusi miring ke kanan (resp. ke kiri) Distribusi yang stabil memungkinkan ekor gemuk ketika berkurang.β[1,1]σ>0μ(,)α(0,2]βμα

Ketika benar-benar kurang dari satu dan , dukungan distribusi terbatas pada .αβ=1(μ,)

Fungsi kerapatan hanya memiliki ekspresi bentuk-tertutup untuk beberapa kombinasi nilai tertentu untuk parameter. Ketika , , , dan σ = α adalah (lihat rumus (4.4) di sini ):μ=0α<1β=1σ=α

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

Ini memiliki mean dan varian yang tak terbatas.

Pertanyaan

Saya ingin menggunakan kerapatan itu dalam R. Saya menggunakan

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

di mana fungsi dstable datang dengan paket fBasics.

Bisakah Anda mengonfirmasi ini adalah cara yang tepat untuk menghitung kepadatan dalam R?

Terima kasih sebelumnya!

EDIT

Salah satu alasan mengapa saya curiga adalah bahwa, dalam output, nilai delta berbeda dari yang ada di input. Contoh:

> 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
okram
sumber

Jawaban:

6

Jawaban singkatnya adalah Anda baik-baik saja, tetapi γ Anda salah. Untuk mendapatkan distribusi stabil positif yang diberikan oleh rumus Anda di R, Anda harus menetapkan γ = | 1 - i tan ( π α / 2 ) | - 1 / α .δγ

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

Contoh paling awal yang dapat saya temukan dari formula yang Anda berikan ada di (Feller, 1971), tetapi saya hanya menemukan buku itu dalam bentuk fisik. Namun (Hougaard, 1986) memberikan rumus yang sama, bersama dengan Transformasi Laplace Dari manual ( digunakan dalam ), the

L(s)=E[exp(sX)]=exp(sα).
stablediststabledistfBasicspm=1parameterisasi adalah dari (Samorodnitsky dan Taqqu, 1994), sumber daya lain yang reproduksi daringnya telah menghindarkan saya. Namun (Weron, 2001) memberikan fungsi karakteristik dalam Samorodnitsky dan parameterisasi Taqqu untuk menjadi φ ( t ) = E [ exp ( i t X ) ] = exp [ i δ t - γ α | t | α ( 1 - i β s i g n ( t )α1 Saya telah mengganti nama beberapa parameter dari kertas Weron menjadi coinide dengan notasi yang kami gunakan. Dia menggunakanμuntukδdanσuntukγ. Bagaimanapun, memasukkanβ=1danδ=0, kita mendapatkan φ(t)=exp[-γα| t| α(1-isign(t)tanπα
φ(t)=E[exp(itX)]=exp[iδtγα|t|α(1iβsign(t)tanπα2)].
μδσγβ=1δ=0
φ(t)=exp[γα|t|α(1isign(t)tanπα2)].

Perhatikan bahwa untuk α ( 0 , 1 ) dan i α = exp ( i π α / 2 ) . Secara formal, L(1itan(πα/2))/|1itan(πα/2)|=exp(iπα/2)α(0,1)iα=exp(iπα/2) , jadi dengan mengatur γ = | 1 - i tan ( π α / 2 ) | - 1 / α dalam φ ( t ) kita dapatkan φ ( i s ) = exp ( - s α ) = L ( s ) . Satu hal yang menarik untuk dicatat adalah bahwa γL(s)=φ(is)γ=|1itan(πα/2)|1/αφ(t)

φ(is)=exp(sα)=L(s).
γyang bersesuaian dengan juga 1 / 2 , jadi jika Anda mencoba γ = α atau γ = 1 - α , yang sebenarnya bukan pendekatan yang buruk, Anda berakhir tepat benar untuk α = 1 / 2 .α=1/21/2γ=αγ=1αα=1/2

Berikut ini contoh dalam R untuk memeriksa kebenaran:

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

Keluaran plot

  1. Feller, W. (1971). Pengantar Teori Probabilitas dan Penerapannya , 2 , 2nd ed. New York: Wiley.
  2. Hougaard, P. (1986). Model Kelangsungan Hidup untuk Populasi Heterogen Berasal dari Distribusi Stabil , Biometrika 73 , 387-396.
  3. Samorodnitsky, G., Taqqu, MS (1994). Proses Acak Non-Gaussian Stabil , Chapman & Hall, New York, 1994.
  4. Weron, R. (2001). Distribusi Levy-stable ditinjau kembali: indeks ekor> 2 tidak mengecualikan rezim Levy-stable , International Journal of Modern Physics C, 2001, 12 (2), 209-223.
P Schnell
sumber
1
Dengan senang hati. Topik parameterisasi stabil positif menyebabkan banyak sakit kepala untuk saya awal tahun ini (benar-benar berantakan), jadi saya memposting apa yang saya hasilkan. Bentuk khusus ini berguna dalam analisis kelangsungan hidup karena bentuk Laplacian memungkinkan hubungan sederhana antara parameter regresi kondisional dan marginal dalam model bahaya proporsional ketika ada istilah kelemahan berikut distribusi stabil yang positif (lihat makalah Hougaard).
P Schnell
6

Apa yang saya pikirkan sedang terjadi adalah bahwa dalam output deltamungkin melaporkan nilai lokasi internal, sedangkan dalam input deltamenggambarkan pergeseran. [Sepertinya ada masalah yang sama dengan gammakapan pm=2.] Jadi jika Anda mencoba meningkatkan shift menjadi 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

maka Anda menambahkan 2 ke nilai lokasi.

Dengan beta=1dan pm=1Anda memiliki variabel acak positif dengan batas bawah distribusi pada 0.

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

Bergeser 2 dan batas bawah naik dengan jumlah yang sama

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

Tetapi jika Anda ingin deltainput menjadi nilai lokasi internal daripada pergeseran atau batas bawah, maka Anda perlu menggunakan spesifikasi yang berbeda untuk parameter. Misalnya jika Anda mencoba yang berikut ini (dengan pm=3dan mencoba delta=0dan yang delta=0.290617Anda temukan sebelumnya), Anda sepertinya mendapatkan yang deltamasuk dan keluar yang sama. Dengan pm=3dan delta=0.290617Anda mendapatkan kerapatan yang sama dengan 0,02700602 yang Anda temukan sebelumnya dan batas bawah pada 0. Dengan pm=3dan delta=0Anda mendapatkan batas bawah negatif (pada kenyataannya -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

Anda mungkin merasa lebih mudah untuk mengabaikan deltadalam output, dan selama Anda tetap beta=1menggunakan pm=1berarti deltadalam input adalah distribusi batas bawah, yang tampaknya Anda ingin menjadi 0.

Henry
sumber
4

Juga dari catatan: Martin Maechler hanya refactored kode untuk distabilkan yang didistribusikan dan menambahkan beberapa perbaikan.

Stabilis paket barunya akan digunakan oleh fBasics juga, jadi Anda mungkin ingin melihatnya juga.

Dirk Eddelbuettel
sumber