Bagaimana cara mengetahui apakah data saya cocok dengan distribusi Pareto?

10

Saya punya sampel yang merupakan vektor dengan 220 angka. Berikut ini tautan ke histogram data saya. . Dan saya ingin memeriksa apakah data saya cocok dengan distribusi Pareto, tetapi saya tidak ingin melihat plot QQ dengan distribusi itu, tetapi saya membutuhkan jawaban yang tepat dengan nilai-p dalam R, seperti uji Anderson-Darling untuk normalitas ( ad.test) . Bagaimana saya bisa melakukan itu? Harap sespesifik mungkin.

bersabar
sumber
1
Hasil uji statistik tidak akan memberi tahu Anda bahwa data Anda memiliki distribusi Pareto . Bahkan Anda bisa sangat yakin bahwa jika itu adalah data nyata, mereka tidak memiliki distribusi Pareto. Semua tes akan menunjukkan kepada Anda apakah Anda memiliki cukup data untuk mengambil deviasi jumlah dari Pareto yang Anda miliki. Artinya, jika menolak semua yang dikatakannya adalah 'yeah ukuran sampel cukup besar untuk memberi tahu Anda apa yang sudah Anda ketahui'. Mengapa Anda melakukan latihan semacam itu, yang tidak dapat menjawab pertanyaan aktual yang Anda miliki?
Glen_b -Reinstate Monica
Apakah pertanyaan Anda benar-benar tidak lebih dari 'baris kode apa yang saya tulis untuk membuat program R melakukan prosedur X'? Maka itu di luar topik di sini. Ini mungkin memenuhi syarat sebagai pertanyaan pemrograman. Jika ada aspek statistik untuk pertanyaan Anda (seperti 'apakah ini masuk akal untuk dilakukan?') Maka Anda harus mengklarifikasi dan menekankan aspek-aspek itu
Glen_b -Reinstate Monica
1
Sekarang untuk tes Anderson-Darling (atau, dalam hal ini, Kolmogorov-Smirnov yang @Zen sarankan di atas). Itu adalah tes untuk distribusi yang ditentukan sepenuhnya . Artinya, agar pengujian memiliki sifat yang diinginkan, Anda harus menentukan apriori ( BUKAN perkiraan ) semua parameter. Jadi Anda tidak dapat menggunakan keduanya untuk latihan ini karena Anda tidak memiliki parameter yang ditentukan sebelumnya. (Agaknya Anda melakukan ini atas saran orang lain. Sangat sulit untuk menjelaskan kesalahpahaman kepada seseorang melalui perantara.)
Glen_b -Reinstate Monica
Untuk apa pengujian ini dilakukan? mis. tindakan apa yang akan berubah tergantung pada apakah Anda menolak atau gagal menolak?
Glen_b -Reinstate Monica
Anda harus selalu melihat plot QQ, terlepas dari motif Anda. Dan Anda seharusnya tidak memfitnah nilai-P "tepat". Tes yang berbeda akan memberi Anda nilai P "tepat" yang berbeda.
Nick Cox

Jawaban:

13

(PS) Pertama-tama saya pikir Glen_b tepat di komentar di atas nya pada kegunaan tes seperti: data riil yang pasti tidak persis Pareto didistribusikan, dan untuk aplikasi praktis yang paling pertanyaan akan "seberapa baik adalah pendekatan Pareto?" - dan plot QQ adalah cara yang baik untuk menunjukkan kualitas perkiraan seperti itu.

Cara apa pun yang dapat Anda lakukan dengan statistik Kolmogorov-Smirnov, setelah memperkirakan parameter dengan kemungkinan maksimum. Estimasi parameter ini mencegah untuk menggunakan -value , sehingga Anda dapat melakukan bootstrap parametrik untuk memperkirakannya. Seperti yang Glen_b katakan dalam komentar, ini dapat dihubungkan ke tes Lilliefors .pks.test

Berikut adalah beberapa baris kode R.

Pertama-tama tentukan fungsi dasar untuk menangani distribusi Pareto.

# distribution, cdf, quantile and random functions for Pareto distributions
dpareto <- function(x, xm, alpha) ifelse(x > xm , alpha*xm**alpha/(x**(alpha+1)), 0)
ppareto <- function(q, xm, alpha) ifelse(q > xm , 1 - (xm/q)**alpha, 0 )
qpareto <- function(p, xm, alpha) ifelse(p < 0 | p > 1, NaN, xm*(1-p)**(-1/alpha))
rpareto <- function(n, xm, alpha) qpareto(runif(n), xm, alpha)

Fungsi berikut menghitung MLE dari parameter (justifikasi dalam Wikipedia ).

pareto.mle <- function(x)
{
  xm <- min(x)
  alpha <- length(x)/(sum(log(x))-length(x)*log(xm))
  return( list(xm = xm, alpha = alpha))
}

Dan fungsi ini menghitung statistik KS, dan menggunakan bootstrap parametrik untuk memperkirakan nilai- .p

pareto.test <- function(x, B = 1e3)
{
  a <- pareto.mle(x)

  # KS statistic
  D <- ks.test(x, function(q) ppareto(q, a$xm, a$alpha))$statistic

  # estimating p value with parametric bootstrap
  B <- 1e5
  n <- length(x)
  emp.D <- numeric(B)
  for(b in 1:B)
  {
    xx <- rpareto(n, a$xm, a$alpha);
    aa <- pareto.mle(xx)
    emp.D[b] <- ks.test(xx, function(q) ppareto(q, aa$xm, aa$alpha))$statistic
  }

  return(list(xm = a$xm, alpha = a$alpha, D = D, p = sum(emp.D > D)/B))
}

Sekarang, misalnya, sampel yang berasal dari distribusi Pareto:

> # generating 100 values from Pareto distribution
> x <- rpareto(100, 0.5, 2)
> pareto.test(x)
$xm
[1] 0.5007593

$alpha
[1] 2.080203

$D
         D 
0.06020594 

$p
[1] 0.69787

... dan dari a :χ2(2)

> # generating 100 values from chi square distribution
> x <- rchisq(100, df=2)
> pareto.test(x)
$xm
[1] 0.01015107

$alpha
[1] 0.2116619

$D
        D 
0.4002694 

$p
[1] 0

Perhatikan bahwa saya tidak mengklaim bahwa tes ini tidak bias: ketika sampel kecil, beberapa bias mungkin ada. Bootstrap parametrik tidak memperhitungkan ketidakpastian pada estimasi parameter (pikirkan apa yang akan terjadi ketika menggunakan strategi ini untuk menguji secara naif jika rata-rata dari beberapa variabel normal dengan varian yang tidak diketahui adalah nol).

PS Wikipedia mengatakan beberapa kata tentang ini. Berikut adalah dua pertanyaan lain yang disarankan untuk strategi yang serupa: Uji goodness of fit untuk campuran , uji goodness of fit untuk distribusi gamma .

Elvis
sumber
3
Saat Anda menyesuaikan distribusi statistik uji untuk estimasi parameter dengan cara ini, ini bukan tes KS (meskipun berdasarkan statistik KS) - ini adalah jenis tes Lilliefors tertentu . Ini bukan lagi nonparametrik, tetapi orang dapat membangun satu melalui simulasi untuk setiap distribusi yang diberikan. Lilliefors melakukan ini khusus untuk yang normal dan eksponensial ... di tahun 1960-an.
Glen_b -Reinstate Monica
Terima kasih atas komentar ini @Glen_b saya tidak tahu itu.
Elvis
Tidak masalah; tidak mengubah apa-apa tentang konten dari apa yang Anda lakukan (yang baik-baik saja), hanya apa yang seharusnya disebut.
Glen_b -Reinstate Monica
@ Glen_b Saya membuat beberapa perubahan substansial dalam jawaban saya, terima kasih lagi!
Elvis