Показать сообщение отдельно
Старый 11.11.2012, 01:32   #46
Hogfather
Platinum Member
 
Аватар для Hogfather
 
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,298
По умолчанию

"Ладно. пора кончать этот бардак. Давайте её закопаем"


Итак, коллеги. Товарищ Дмитрий В. получил интересные результаты, стал мучать их в Excel и получил картинки, которые никуда не годятся. Дла начала, у нас распределение явно дискретное, а мы рисуем график как для непрерывного. Зачем точки соединять то?
Плюнем на Excel слюною, пусть в нем, товарищи, успешные менеджеры отчеты делают, нам путь в нормальный статистический пакет, поэтому только хардкор, только R.
Устанавливаем R, создаем вектор данных.
Код:
> LT<-c(rep(1,9),rep(2,267),rep(3,2843),rep(4,5450),rep(5,6564),rep(6,7044),rep(7,7518),rep(8,7071),rep(9,5620),rep(10,4016),rep(11,2545),rep(12,1494),rep(13,854),rep(14,416),rep(15,214),rep(16,122),rep(17,53),rep(18,16),rep(19,7),rep(20,2),21,22)
Функция rep повторяет первый аргумент число раз, равное второму аргументу, поэтому для rep(1,9) имеем в результате вектор [1,1,1,1,1,1,1,1,1].
Данные берем с графика, любезно предоставленного нам.
Смотрим на результат и радуемся
Код:
> summary(LT)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     5.0     7.0     7.2     9.0    22.0 
> length(LT)
[1] 52127
Дальше мы просто обязаны поффтыкать на картинки, иначе мы не ученые, а кот начхал.
Сказано-сделано, строим 4 графика в одном.
Код:
> old.par <- par(mfrow=c(2,2))
> hist(LT,main="Распределение букв",ylab="Число наблюдений",xlab="Число букв")
> hist(LT,freq=F,ylab="Вероятность",xlab="Число букв",main="Распределение букв")
> plot(ecdf(LT),verticals=T,main="График функции распределения")
> boxplot(LT,main="Диаграмма Ящик-с-Усами",xlab="Число букв",horizontal=T)
> par(old.par)


Что мы, собственно говоря видим. А видим, что распределение у нас вполне милое, да слегка несимметричное, но с кем не бывает.
Пытаемся натянуть сову на глобус. Для этого используем подгонку распределения методом максимального правдоподобия (maximum-likelihood estimation, MLE). Метод это весьма кошерен, но связан со сложными вычислениями. К счастью для нас, в R уже всё таки имеется. Достаточно подключить библиотеку MASS.

Резвимся по полной
Код:
> library(MASS)
Предупреждение
пакет ‘MASS’ был собран под R версии 2.14.2 

> fitdistr(LT, "gamma")
      shape         rate    
  7.257622928   1.008025740 
 (0.043960087) (0.006321817)
Предупреждения
1: In dgamma(x, shape, scale, log) : созданы NaN
2: In dgamma(x, shape, scale, log) : созданы NaN
3: In dgamma(x, shape, scale, log) : созданы NaN
4: In dgamma(x, shape, scale, log) : созданы NaN
5: In dgamma(x, shape, scale, log) : созданы NaN
6: In dgamma(x, shape, scale, log) : созданы NaN
7: In dgamma(x, shape, scale, log) : созданы NaN

> fitdistr(LT, "normal")
      mean           sd     
  7.199838855   2.628803586 
 (0.011514015) (0.008141638)

> fitdistr(LT,"lognormal")
      meanlog           sdlog 
  1.903586097   0.385993556 
 (0.001690630) (0.001195456)

> fitdistr(LT, "Poisson")
     lambda  
  7.19983886 
 (0.01175249)
Итак, мы что-то наподгоняли. Попробовали гамму, нормальное, логнормальное и Пуассона. В скобках, для удобства, дана ошибка параметров.

Неплохим графическим методом оценки качества подгонки распределения является график квантилей (quantile). Квантиль — это такое число, что заданная случайная величина не превышает его лишь с указанной вероятностью. Можно рассматривать квантиль как функцию вероятности Q(p), обратную функции распределения вероятностей. Если мы подогнали правильно, то точки на графике должны лежать рядом с прямой y = x. Строим четыре графика для наших распределений.

Код:
> old.par <- par(mfrow=c(2,2))
> qqplot(LT, rgamma(n = 52127, 7.257622928, 1.008025740), main = "Подгонка гамма-распределения, QQ-plot")
> abline(0, 1)
> qqplot(LT, rpois(n = 52127, 7.19983886), main = "Подгонка распределения Пуассона, QQ-plot")
> abline(0, 1)
> qqplot(LT, rnorm(n = 52127, 7.199838855,2.628803586), main = "Подгонка нормального распределения, QQ-plot")
> abline(0, 1)
> qqplot(LT, rlnorm(n = 52127, 1.903586097, 0.385993556), main = "Подгонка Логнормального распределения, QQ-plot")
> abline(0, 1)
> par(old.par)


Кому как, а мне больше нравится старик Пуассон. Попробуем нарисовать график аппроксимирующих распределений.

Код:
> plot(ecdf(LT),verticals=T,main="Аппроксимация функции распределения")
> lines(0:2200/100,pgamma(0:2200/100,7.257622928, 1.008025740),col="red")
> lines(0:2200/100,ppois(0:2200/100,lambda=7.19983886),col="blue",lwd=2)
> legend(15,0.2,c("Гамма распределение","Распределение Пуассона"),col=c("red","blue"),lwd=2)
Результат


Ну, пока хватит. Коню понятно, что здесь никакая не гамма, а обычный Пуассон, причем Лямбда равна среднему числу букв в слове.
Ну, теперь сам бог велел провести тест Колмогорова-Смирнова
Код:
> ks.test(LT,rpois(0:2200/100,lambda=7.19983886))

        Two-sample Kolmogorov-Smirnov test

data:  LT and rpois(0:2200/100, lambda = 7.19983886) 
D = 0.0261, p-value = 0.1137
alternative hypothesis: two-sided 

Предупреждение
In ks.test(LT, rpois(0:2200/100, lambda = 7.19983886)) :
  p-values будут примерными в присутствии повторяющихся значений
Ай, да Hogfather, хочется сказать, ай, молодец!

Бурные продолжительные аплодисменты.

А гамма ваша, кака редкая...

Код:
> ks.test(LT,rgamma(0:2200/100, 7.257622928, 1.008025740))

        Two-sample Kolmogorov-Smirnov test

data:  LT and rgamma(0:2200/100, 7.257622928, 1.00802574) 
D = 0.1011, p-value < 2.2e-16
alternative hypothesis: two-sided 

Предупреждение
In ks.test(LT, rgamma(0:2200/100, 7.257622928, 1.00802574)) :
  p-values будут примерными в присутствии повторяющихся значений



Согласен на соавторство


P.S. Ну, мои маленькие девиантные друзья, если кто хочет поподробнее почитать про подгонку распределений в R, рекомендую на сон грядущий статью "Fitting distributions with R"


P.P.S. А список наиболее распространенных распределений можно посмотреть вот тут, в вашей любимой Википедии

Последний раз редактировалось Hogfather; 11.11.2012 в 02:21.
---------
DNF is not an option
Hogfather вне форума   Ответить с цитированием
Реклама