Как мне разместить набор данных для распределения Парето в R?

Хорошо, если у вас есть образец из распределения Парето с параметрами и (где - параметр нижней границы, а - параметр формы), логарифмическая вероятность этого образец является: $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)$

это монотонно возрастает в , поэтому максимальный максимальный показатель является наибольшим значением, которое согласуется с наблюдаемыми данными. Поскольку параметр определяет нижнюю границу носителя для распределения Парето, оптимум $m$ $m$

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

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

который не зависит от . Далее, используя обычные трюки исчисления, MLE для должно удовлетворять $\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$

некоторая простая алгебра говорит нам ОМП является $\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})}$

Во многих важных смыслах (например, оптимальная асимптотическая эффективность в том смысле, что достигается нижняя граница Крамера-Рао), это лучший способ согласовать данные с распределением Парето. Приведенный ниже код R вычисляет MLE для данного набора данных 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

Редактировать: Основываясь на комментариях @cardinal и I ниже, мы также можем заметить, что является обратной величиной среднего значения выборки , которые происходят с имеют экспоненциальное распределение. Следовательно, если у нас есть доступ к программному обеспечению, которое может соответствовать экспоненциальному распределению (что более вероятно, так как оно возникает во многих статистических задачах), то подгонка распределения Парето может быть достигнута путем преобразования набора данных таким образом и подгонки его экспоненциальному распределению в преобразованном масштабе. $\hat{\alpha}$ $\log(X_i /\hat{m})$

макрос
источник

(+1) Мы можем написать что-то более внушительно, заметив, что распределяется экспоненциально со скоростью . Исходя из этого и неизменности MLE при преобразовании, мы сразу же заключаем, что , где мы заменяем на в последнем выражении. Это также намекает на то, как мы могли бы использовать стандартное программное обеспечение для соответствия Парето, даже если нет явной опции.

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$

кардинал

@cardinal - Таким образом, является обратной величиной среднего значения выборки для , которые имеют экспоненциальное распределение. Как это поможет нам?

\hat{α}

$\hat{\alpha}$

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

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

Макро

Привет Макро. Я пытался подчеркнуть, что проблема оценки параметров Парето может быть (по существу) сведена к оценке скорости экспоненты: с помощью приведенного выше преобразования мы можем преобразовать наши данные и задачу в (возможно) более знакомый и сразу же извлеките ответ (при условии, что мы или наше программное обеспечение уже знаем, что делать с образцом экспонент).

кардинал

Как я могу измерить ошибку такого рода подгонки?

Эмануэле

@emanuele, приблизительная дисперсия MLE является инверсией информационной матрицы Фишера, которая потребует от вас вычисления хотя бы одной производной логарифмического правдоподобия. Или вы могли бы использовать своего рода загрузочную передискретизацию для оценки стандартной ошибки.

Макро

Вы можете использовать fitdistфункцию, представленную в fitdistrplusпакете:

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
источник

Должно ли это быть library(fitdistrplus)?

Шон

@ Да, да, редактирую ответ соответственно

Кевин Л Киз

Обратите внимание, что library(actuar)для этого требуется вызов .

Jsta

Что представляет собой fp $ оценка ["форма"] в этом случае? Возможно, это примерная альфа? Или бета?

Альберт Хендрикс

Как мне разместить набор данных для распределения Парето в R?

Ответы: