Jak dopasować zestaw danych do rozkładu Pareto w R?

Mieć, powiedzmy, następujące dane:

8232302  684531  116857   89724   82267   75988   63871   
  23718    1696     436     439     248     235

Potrzebujesz prostego sposobu dopasowania tego (i kilku innych zestawów danych) do dystrybucji Pareto. Idealnie byłoby wyprowadzić pasujące wartości teoretyczne, mniej idealnie parametry.

r pareto-distribution

— Felix
źródło

cran.r-project.org/web/packages/fitdistrplus/fitdistrplus.pdf

— Stéphane Laurent

Co należy rozumieć przez „dopasowanie wartości teoretycznych”? Oczekiwania dotyczące statystyk zamówień, biorąc pod uwagę oszacowania parametrów? Albo coś innego?

— Glen_b

Cóż, jeśli masz próbkę z rozkładu pareto z parametrami i (gdzie jest dolnym ograniczonym parametrem, a jest parametrem kształtu), prawdopodobieństwo logarytmu tego próbka to: $X_1, ..., X_n$ $m>0$ $\alpha>0$ $m$ $\alpha$

n \log (α) + n α \log (m) - (α + 1) \sum_{i = 1}^{n} \log (X_{i})

$n \log(\alpha) + n \alpha \log(m) - (\alpha+1) \sum_{i=1}^{n} \log(X_i)$

jest to monotonicznie rosnąca , więc maksymalizator jest największą wartością zgodną z obserwowanymi danymi. Ponieważ parametr określa dolną granicę podparcia dla rozkładu Pareto, optymalne jest $m$ $m$

\hat{m} = min_{i} X_{i}

$\hat{m} = \min_{i} X_i$

co nie zależy od . Następnie, używając zwykłych sztuczek rachunku różniczkowego, MLE dla musi spełniać $\alpha$ $\alpha$

\frac{n}{α} + n \log (\hat{m}) - \sum_{i = 1}^{n} \log (X_{i}) = 0

$\frac{n}{\alpha} + n \log( \hat{m} ) - \sum_{i=1}^{n} \log(X_i) = 0$

jakaś prosta algebra mówi nam, że MLE z jest $\alpha$

\hat{α} = \frac{n}{\sum_{i = 1}^{n} \log (X_{i} / \hat{m})}

$\hat{\alpha} = \frac{n}{\sum_{i=1}^{n} \log(X_i/\hat{m})}$

W wielu ważnych aspektach (np. Optymalna wydajność asymptotyczna, ponieważ osiąga dolną granicę Cramer-Rao), jest to najlepszy sposób dopasowania danych do rozkładu Pareto. Kod R poniżej oblicza MLE dla danego zestawu danych X.

pareto.MLE <- function(X)
{
   n <- length(X)
   m <- min(X)
   a <- n/sum(log(X)-log(m))
   return( c(m,a) ) 
}

# example. 
library(VGAM)
set.seed(1)
z = rpareto(1000, 1, 5) 
pareto.MLE(z)
[1] 1.000014 5.065213

Edycja: Na podstawie komentarza @cardinal i I poniżej możemy również zauważyć, że jest odwrotnością średniej próbki z , które zdarzają się mieć rozkład wykładniczy. Dlatego jeśli mamy dostęp do oprogramowania, które może pasować do rozkładu wykładniczego (co jest bardziej prawdopodobne, ponieważ wydaje się, że pojawia się w wielu problemach statystycznych), wówczas dopasowanie rozkładu Pareto można osiągnąć poprzez transformację zestawu danych w ten sposób i dopasowanie go do rozkładu wykładniczego w przekształconej skali. $\hat{\alpha}$ $\log(X_i /\hat{m})$

— Makro
źródło

(+1) Możemy pisać rzeczy bardziej sugestywnie, zauważając, że jest rozkładany wykładniczo ze współczynnikiem . z tego i niezmienności transformowanych MLE wywodzimy od razu, że , gdzie zamieniamy na w tym ostatnim wyrażeniu. Wskazuje to również, w jaki sposób możemy użyć standardowego oprogramowania, aby dopasować Pareto, nawet jeśli nie jest dostępna żadna wyraźna opcja.

Y_{i} = \log (X_{i} / m)

$Y_i = \log(X_i/m)$

α

$\alpha$

\hat{α} = 1 / \bar{Y}

$\hat\alpha = 1/\bar Y$

m

$m$

\hat{m}

$\hat m$

— kardynał

@cardinal - Tak więc jest odwrotnością średniej próbki , które mają rozkład wykładniczy. Jak nam to pomaga?

\hat{α}

$\hat{\alpha}$

\log (X_{i} / \hat{m})

$\log(X_i/\hat{m})$

— Makro

Cześć, Makro. Chodziło mi o to, że problem oszacowania parametrów Pareto można (zasadniczo) zredukować do problemu oszacowania współczynnika wykładniczego: za pomocą powyższej transformacji możemy przekonwertować nasze dane i problem na (być może) bardziej znajomy i natychmiast wyodrębnij odpowiedź (zakładając, że my lub nasze oprogramowanie wiemy już, co zrobić z próbką wykładniczą).

— kardynał

Jak mogę zmierzyć błąd tego rodzaju dopasowania?

— emanuele

@emanuele, przybliżona wariancja MLE jest odwrotnością macierzy informacji Fishera, która będzie wymagać obliczenia co najmniej jednej pochodnej prawdopodobieństwa log. Lub możesz użyć pewnego rodzaju próbkowania ładowania początkowego, aby oszacować standardowy błąd.

— Makro

Możesz użyć fitdistfunkcji dostarczonej w fitdistrpluspakiecie:

library(MASS)
library(fitdistrplus)
library(actuar)

# suppose data is in dataPar list
fp <- fitdist(dataPar, "pareto", start=list(shape = 1, scale = 500))
#the mle parameters will be stored in fp$estimate

— akashrajkn
źródło

Tak powinno być library(fitdistrplus)?

— Sean

@Sean tak, odpowiednio edytujesz odpowiedź

— Kevin L Keys

Pamiętaj, że połączenie library(actuar)jest wymagane, aby to zadziałało.

— jsta

Co w tym przypadku reprezentuje fp $ oszacowanie [„kształt”]? Czy to może oszacowana alfa? A może beta?

— Albert Hendriks,