Как определить, какое распределение лучше всего подходит для моих данных?

133

У меня есть набор данных, и я хочу выяснить, какое распределение лучше всего подходит для моих данных.

Я использовал fitdistr()функцию для оценки необходимых параметров для описания предполагаемого распределения (т. Е. Вейбулла, Коши, Нормаль). Используя эти параметры, я могу провести тест Колмогорова-Смирнова, чтобы оценить, соответствуют ли мои выборочные данные тому же распределению, что и предполагаемое распределение.

Если значение р> 0,05, я могу предположить, что данные выборки взяты из того же распределения. Но значение р не дает никакой информации о божественности посадки, не так ли?

Таким образом, если значение p в моих выборочных данных составляет> 0,05 для нормального распределения, а также распределения Вейбулла, как я могу узнать, какое распределение лучше соответствует моим данным?

Это в основном то, что я сделал:

> mydata
 [1] 37.50 46.79 48.30 46.04 43.40 39.25 38.49 49.51 40.38 36.98 40.00
[12] 38.49 37.74 47.92 44.53 44.91 44.91 40.00 41.51 47.92 36.98 43.40
[23] 42.26 41.89 38.87 43.02 39.25 40.38 42.64 36.98 44.15 44.91 43.40
[34] 49.81 38.87 40.00 52.45 53.13 47.92 52.45 44.91 29.54 27.13 35.60
[45] 45.34 43.37 54.15 42.77 42.88 44.26 27.14 39.31 24.80 16.62 30.30
[56] 36.39 28.60 28.53 35.84 31.10 34.55 52.65 48.81 43.42 52.49 38.00
[67] 38.65 34.54 37.70 38.11 43.05 29.95 32.48 24.63 35.33 41.34

# estimate shape and scale to perform KS-test for weibull distribution
> fitdistr(mydata, "weibull")
     shape        scale   
   6.4632971   43.2474500 
 ( 0.5800149) ( 0.8073102)

# KS-test for weibull distribution
> ks.test(mydata, "pweibull", scale=43.2474500, shape=6.4632971)

        One-sample Kolmogorov-Smirnov test

data:  mydata
D = 0.0686, p-value = 0.8669
alternative hypothesis: two-sided

# KS-test for normal distribution
> ks.test(mydata, "pnorm", mean=mean(mydata), sd=sd(mydata))

        One-sample Kolmogorov-Smirnov test

data:  mydata
D = 0.0912, p-value = 0.5522
alternative hypothesis: two-sided

Значения p составляют 0,8669 для распределения Вейбулла и 0,5522 для нормального распределения. Таким образом, я могу предположить, что мои данные следуют как по Вейбуллу, так и по нормальному распределению. Но какая функция распределения лучше описывает мои данные?


Обращаясь к elevendollar, я нашел следующий код, но не знаю, как интерпретировать результаты:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
tobibo
источник
5
Почему вы хотите выяснить, какой дистрибутив лучше всего подходит для ваших данных?
Роланд
6
Потому что я хочу генерировать псевдослучайные числа, следующие данному распределению.
тобибо
6
Вы не можете использовать KS, чтобы проверить, соответствует ли распределение с параметрами, найденными из набора данных, набором данных. См. № 2 на этой странице, например, плюс альтернативы (и другие способы, которыми тест KS может вводить в заблуждение).
tpg2114
Другое обсуждение здесь с примерами кода о том, как применить тест KS, когда параметры оцениваются из образца.
Аксакал
1
I used the fitdistr() function ..... что такое fitdistrфункция? Что-то из Excel? Или вы что-то написали сами на С?
волки

Ответы:

163

Во-первых, вот несколько быстрых комментариев:

  • В -значение из в Kolmovorov-Смирнов-Test (KS-Test) с расчетными параметрами будет совершенно неправильно. Так что, к сожалению, вы не можете просто подобрать распределение, а затем использовать оценочные параметры в тесте Колмогорова-Смирнова для проверки вашей выборки.p
  • Ваш образец никогда не будет точно следовать определенному распределению. Таким образом, даже если ваши из KS-Test будут действительными и , это будет означать, что вы не можете исключить, что ваши данные соответствуют этому конкретному распределению. Другая формулировка будет заключаться в том, что ваш образец совместим с определенным распределением. Но ответ на вопрос "Мои данные точно соответствуют распределению xy?" всегда нет.p>0.05
  • Цель здесь не может состоять в том, чтобы с уверенностью определить, за каким распределением следует ваша выборка. Цель - то, что @whuber (в комментариях) называет скупыми приблизительными описаниями данных. Наличие конкретного параметрического распределения может быть полезным в качестве модели данных.

Но давайте сделаем некоторые исследования. Я буду использовать отличный fitdistrplusпакет, который предлагает несколько хороших функций для настройки распределения. Мы будем использовать функцию, descdistчтобы получить некоторые идеи о возможных распределениях кандидатов.

library(fitdistrplus)
library(logspline)

x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00,
38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40,
42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40,
49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60,
45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30,
36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00,
38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34)

Теперь давайте используем descdist:

descdist(x, discrete = FALSE)

Descdist

Куртоз и квадратная асимметрия вашего образца представлены в виде синей точки с именем «Наблюдение». Кажется, что возможные распределения включают распределение Вейбулла, Логнормального и, возможно, гамма-распределения.

Давайте подойдем к распределению Вейбулла и нормальному распределению:

fit.weibull <- fitdist(x, "weibull")
fit.norm <- fitdist(x, "norm")

Теперь осмотрите на предмет соответствия нормальному:

plot(fit.norm)

Нормальная посадка

А для Вейбулла подойдет:

plot(fit.weibull)

Weibull Fit

Оба выглядят хорошо, но, судя по QQ-Plot, Weibull, возможно, выглядит немного лучше, особенно на хвостах. Соответственно, AIC соответствия Вейбулла ниже по сравнению с нормальным соответствием:

fit.weibull$aic
[1] 519.8537

fit.norm$aic
[1] 523.3079

Тестовое моделирование Колмогорова-Смирнова

Я буду использовать процедуру @ Aksakal, описанную здесь, чтобы смоделировать статистику KS под нулем.

n.sims <- 5e4

stats <- replicate(n.sims, {      
  r <- rweibull(n = length(x)
                , shape= fit.weibull$estimate["shape"]
                , scale = fit.weibull$estimate["scale"]
  )
  estfit.weibull <- fitdist(r, "weibull") # added to account for the estimated parameters
  as.numeric(ks.test(r
                     , "pweibull"
                     , shape= estfit.weibull$estimate["shape"]
                     , scale = estfit.weibull$estimate["scale"])$statistic
  )      
})

ECDF моделируемой статистики KS выглядит следующим образом:

plot(ecdf(stats), las = 1, main = "KS-test statistic simulation (CDF)", col = "darkorange", lwd = 1.7)
grid()

Имитация KS-статистики

Наконец, наше значение, использующее смоделированное нулевое распределение KS-статистики:p

fit <- logspline(stats)

1 - plogspline(ks.test(x
                       , "pweibull"
                       , shape= fit.weibull$estimate["shape"]
                       , scale = fit.weibull$estimate["scale"])$statistic
               , fit
)

[1] 0.4889511

Это подтверждает наш графический вывод о том, что выборка совместима с распределением Вейбулла.

Как объяснено здесь , мы можем использовать начальную загрузку, чтобы добавить точечные доверительные интервалы к оцененному Weibull PDF или CDF:

xs <- seq(10, 65, len=500)

true.weibull <- rweibull(1e6, shape= fit.weibull$estimate["shape"]
                         , scale = fit.weibull$estimate["scale"])

boot.pdf <- sapply(1:1000, function(i) {
  xi <- sample(x, size=length(x), replace=TRUE)
  MLE.est <- suppressWarnings(fitdist(xi, distr="weibull"))  
  dweibull(xs, shape=MLE.est$estimate["shape"],  scale = MLE.est$estimate["scale"])
}
)

boot.cdf <- sapply(1:1000, function(i) {
  xi <- sample(x, size=length(x), replace=TRUE)
  MLE.est <- suppressWarnings(fitdist(xi, distr="weibull"))  
  pweibull(xs, shape= MLE.est$estimate["shape"],  scale = MLE.est$estimate["scale"])
}
)   

#-----------------------------------------------------------------------------
# Plot PDF
#-----------------------------------------------------------------------------

par(bg="white", las=1, cex=1.2)
plot(xs, boot.pdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.pdf),
     xlab="x", ylab="Probability density")
for(i in 2:ncol(boot.pdf)) lines(xs, boot.pdf[, i], col=rgb(.6, .6, .6, .1))

# Add pointwise confidence bands

quants <- apply(boot.pdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.pdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.pdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)

CI_Density

#-----------------------------------------------------------------------------
# Plot CDF
#-----------------------------------------------------------------------------

par(bg="white", las=1, cex=1.2)
plot(xs, boot.cdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.cdf),
     xlab="x", ylab="F(x)")
for(i in 2:ncol(boot.cdf)) lines(xs, boot.cdf[, i], col=rgb(.6, .6, .6, .1))

# Add pointwise confidence bands

quants <- apply(boot.cdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.cdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.cdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)
#lines(xs, min.point, col="purple")
#lines(xs, max.point, col="purple")

CI_CDF


Автоматическая распределительная арматура с GAMLSS

gamlssПакет Rпредлагает возможность попробовать множество различных распределений и выбрать «лучший» в соответствии с GAIC (обобщенная информация Akaike критерий). Основная функция есть fitDist. Важной опцией в этой функции является тип проверяемых дистрибутивов. Например, setting type = "realline"будет пробовать все реализованные распределения, определенные на всей реальной строке, тогда как type = "realsplus"будут пытаться только распределения, определенные на реальной положительной строке. Другим важным вариантом является параметр , который является штрафом для GAIC. В приведенном ниже примере я установил параметр который означает, что «лучшее» распределение выбирается в соответствии с классическим AIC. Вы можете установить на что угодно, напримерkk=2klog(n) для BIC.

library(gamlss)
library(gamlss.dist)
library(gamlss.add)

x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00,
       38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40,
       42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40,
       49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60,
       45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30,
       36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00,
       38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34)

fit <- fitDist(x, k = 2, type = "realplus", trace = FALSE, try.gamlss = TRUE)

summary(fit)

*******************************************************************
Family:  c("WEI2", "Weibull type 2") 

Call:  gamlssML(formula = y, family = DIST[i], data = sys.parent()) 

Fitting method: "nlminb" 


Coefficient(s):
             Estimate  Std. Error  t value   Pr(>|t|)    
eta.mu    -24.3468041   2.2141197 -10.9962 < 2.22e-16 ***
eta.sigma   1.8661380   0.0892799  20.9021 < 2.22e-16 ***

Согласно AIC, распределение Вейбулла (точнее WEI2, его специальная параметризация) наилучшим образом соответствует данным. Точная параметризация распределения WEI2подробно описана в этом документе на стр. 279. Давайте проверим соответствие, посмотрев на остатки на графике червя (в основном, детализированный график QQ):

WormPlot

Мы ожидаем, что остатки будут близки к средней горизонтальной линии, и 95% из них будут лежать между верхней и нижней пунктирными кривыми, которые действуют как 95% точечные доверительные интервалы. В этом случае график червя выглядит хорошо для меня, указывая на то, что распределение Вейбулла адекватно подходит.

COOLSerdash
источник
1
+1 Хороший анализ. Хотя один вопрос. Позволяет ли положительное заключение о совместимости с конкретным основным распределением (в данном случае Вейбуллом) исключить возможность наличия распределения смеси? Или нам нужно выполнить надлежащий анализ смеси и проверить GoF, чтобы исключить этот вариант?
Александр Блех
18
@AleksandrBlekh Невозможно обладать достаточной мощностью, чтобы исключить смесь: когда смесь имеет два почти идентичных распределения, ее невозможно обнаружить, а когда все компоненты, кроме одного, имеют очень малые пропорции, они также не могут быть обнаружены. Как правило (в отсутствие теории, которая могла бы предложить распределительную форму), каждый подходит для параметрических распределений, чтобы получить скупое приблизительное описание данных. Смеси не являются ни одним из них: они требуют слишком много параметров и являются слишком гибкими для этой цели.
whuber
4
@whuber: +1 Ценю ваше отличное объяснение!
Александр Блех
1
@Lourenco Я посмотрел на график Каллена и Фей. Синяя точка обозначает наш образец. Вы видите, что точка близка к линиям Вейбулла, Логнормала и Гаммы (то есть между Вейбуллом и Гаммой). После подбора каждого из этих распределений я сравнил статистику соответствия качества, используя функцию gofstatи AIC. Не существует единого мнения о том, как лучше определить «лучший» дистрибутив. Мне нравятся графические методы и АПК.
COOLSerdash
1
@Lourenco Ты имеешь в виду логнорм? Логистическое распределение (знак «+») немного отличается от наблюдаемых данных. Логнормаль также будет кандидатом, на которого я обычно смотрю. В этом уроке я решил не показывать его, чтобы пост был коротким. Логнормальное показывает худшее соответствие по сравнению как с распределением Вейбулла, так и с нормальным распределением. AIC составляет 537,59, и графики также выглядят не очень хорошо.
COOLSerdash
15

Графики - это, в основном, хороший способ лучше понять, как выглядят ваши данные. В вашем случае я бы порекомендовал построить эмпирическую кумулятивную функцию распределения (ecdf) относительно теоретических cdfs с параметрами, которые вы получили из fitdistr ().

Я сделал это один раз для своих данных, а также включил доверительные интервалы. Вот картинка, которую я получил, используя ggplot2 ().

введите описание изображения здесь

Черная линия - это эмпирическая кумулятивная функция распределения, а цветные линии - это файлы cdf из разных распределений с использованием параметров, которые я получил, используя метод максимального правдоподобия. Легко видеть, что экспоненциальное и нормальное распределение не очень хорошо подходят для данных, потому что линии имеют другую форму, чем ecdf, и линии довольно далеко от ecdf. К сожалению, другие дистрибутивы довольно близки. Но я бы сказал, что линия logNormal ближе всего к черной линии. Используя меру расстояния (например, MSE), можно проверить предположение.

Если у вас есть только два конкурирующих распределения (например , собирание те , которые , кажется, подходят лучше всего в сюжете) , вы можете использовать вероятностно-Ratio-тест , чтобы проверить , какие дистрибутивы лучше подходит.

elevendollar
источник
20
Добро пожаловать в CrossValidated! Ваш ответ мог бы быть более полезным, если бы вы могли отредактировать его, включив в него (а) код, который вы использовали для создания графики, и (б) способ ее чтения.
Стефан Коласса
2
Что там замышляется? Это что-то вроде сюжета экспоненциальности?
Glen_b
1
Но как вы решаете, какой дистрибутив лучше всего подходит для ваших данных? Только по графику я не могу сказать вам, подходит ли logNormal или weibull к вашим данным.
тобибо
4
Если вы хотите создать генератор псевдослучайных чисел, почему бы не использовать эмпирический cdf? Вы хотите нарисовать цифры, которые выходят за рамки вашего наблюдаемого распределения?
одиннадцать долларов
6
Принимая ваш график по номиналу, может показаться, что ни одно из ваших возможных распределений вообще не вписывается в данные. Кроме того, ваш ecdf имеет горизонтальную асимптоту менее 0,03, что не имеет смысла, поэтому я не уверен, что это действительно ecdf.
Хонг Ой