Портал аспирантов - Параметрический анализ

"Нет, всё-таки откопаем..."

Итак, более корректная и интересная первая подгонка, поскольку во втором случае мы просто подгоняем функцию по 22 точкам. С интересом обнаружил, что в 2012 году для R вышел более мощный пакет подгонки fitdistrplus.

Попробуем в него поиграть. Опять берем Гамму.

Код:

>#Подключаем библиотеку

>library(fitdistrplus)

># Подгоняем гамма-распределение

> XX<-fitdist(LT, "gamma")

> summary(XX)

Fitting of the distribution ' gamma ' by maximum likelihood 

Parameters : 

      estimate  Std. Error

shape 7.258422 0.043965053

rate  1.008084 0.006322175

Loglikelihood:  -122729.7   AIC:  245463.4   BIC:  245481.2 

Correlation matrix:

          shape      rate

shape 1.0000000 0.9658169

rate  0.9658169 1.0000000

Ну и Пуассона

Код:

># Подгоняем распределение Пуассона

> XY<-fitdist(LT, "pois")

> summary(XY)

Fitting of the distribution ' pois ' by maximum likelihood 

Parameters : 

       estimate Std. Error

lambda 7.199839 0.01175249

Loglikelihood:  -123149.4   AIC:  246300.8   BIC:  246309.7

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

Пакет позволяет построить красивые картинки. Причем очень просто.

Код:

># Рисунок для гаммы

> plot(XX)

># Рисунок для Пуассона

> plot(XY)

http://aspirantura.spb.ru/forum/pict...&pictureid=975
Рисунок 1 -- Подгонка гамма-распределения

http://aspirantura.spb.ru/forum/pict...&pictureid=976
Рисунок 2 -- Подгонка распределения Пуассона

Рисунок 2 можно также получить не прибегая к построению модели.
Для распределения Пуассона с лямбдой равной средней длине слова это выглядит так:

Код:

> plotdist(LT,"pois",para=list(lambda=mean(LT)))

А можно легко и непринужденно посчитать статистические параметры и проверить гипотезы.

Код:

># Для гамма-распределения

>  gofstat(XX,print.test=TRUE)

Kolmogorov-Smirnov statistic:  0.09400709 

Kolmogorov-Smirnov test:  rejected 

   The result of this test may be too conservative as it  

   assumes that the distribution parameters are known

Cramer-von Mises statistic:  68.65376 

Cramer-von Mises test:  rejected 

Anderson-Darling statistic:  397.2767 

Anderson-Darling test:  rejected 



># Для Распределения Пуассона

> g2 <- gofstat(XY,print.test=TRUE)

Chi-squared statistic:  445.9628 

Degree of freedom of the Chi-squared distribution:  11 

Chi-squared p-value:  1.041315e-88 

> g2$chisqtable

      obscounts theocounts

<= 3  3119.0000  3749.2137

<= 4  5450.0000  4358.0510

<= 5  6564.0000  6275.4530

<= 6  7044.0000  7530.3751

<= 7  7518.0000  7745.3553

<= 8  7071.0000  6970.6637

<= 9  5620.0000  5576.4062

<= 10 4016.0000  4014.9226

<= 11 2545.0000  2627.8905

<= 12 1494.0000  1576.6990

<= 13  854.0000   873.2291

<= 14  416.0000   449.0792

> 14   416.0000   379.6615

>

Первоначальный выбор возможного распределения также осуществляется легко и непринужденно, посчитав моменты.

Код:

> descdist(LT)

summary statistics

------

min:  1   max:  22 

median:  7 

mean:  7.199839 

estimated sd:  2.628829 

estimated skewness:  0.519882 

estimated kurtosis:  3.143716

Вот такая красота.
http://aspirantura.spb.ru/forum/pict...&pictureid=977

Но, поскольку у нас распределение дискретное, мы нарисуем другой график.

Код:

> descdist(LT,discrete = TRUE,boot=1000)

summary statistics

------

min:  1   max:  22 

median:  7 

mean:  7.199839 

estimated sd:  2.628829 

estimated skewness:  0.519882 

estimated kurtosis:  3.143716

http://aspirantura.spb.ru/forum/pict...&pictureid=981

Почти Пуассон, красота!

В общем, пакет мне понравился. Буду пользоваться.

P.S. Если как положено считать Хи-квадрат для дискретного распределения, то видно, что и распределение Пуассона не торт.

Еще разные распределения

Код:

> XZ<-fitdist(LT,"beta")

Ошибка в mledist(data, distname, start, fix.arg, ...) : 

  values must be in [0-1] to fit a beta distribution

> XZ<-fitdist(LT/52127,"beta")

Предупреждения

1: In dbeta(x, shape1, shape2, log) : созданы NaN

2: In dbeta(x, shape1, shape2, log) : созданы NaN

3: In dbeta(x, shape1, shape2, log) : созданы NaN

4: In dbeta(x, shape1, shape2, log) : созданы NaN

5: In dbeta(x, shape1, shape2, log) : созданы NaN

6: In dbeta(x, shape1, shape2, log) : созданы NaN

7: In dbeta(x, shape1, shape2, log) : созданы NaN

8: In dbeta(x, shape1, shape2, log) : созданы NaN

9: In dbeta(x, shape1, shape2, log) : созданы NaN

10: In dbeta(x, shape1, shape2, log) : созданы NaN

> summary(XZ)

Fitting of the distribution ' beta ' by maximum likelihood 

Parameters : 

           estimate   Std. Error

shape1     7.257806   0.01867214

shape2 52538.205482 114.78284503

Loglikelihood:  443444.6   AIC:  -886885.2   BIC:  -886867.5 

Correlation matrix:

          shape1    shape2

shape1 1.0000000 0.7921102

shape2 0.7921102 1.0000000



> gofstat(XZ,print.test=TRUE)

Kolmogorov-Smirnov statistic:  0.09402943 

Kolmogorov-Smirnov test:  rejected 

   The result of this test may be too conservative as it  

   assumes that the distribution parameters are known

Cramer-von Mises statistic:  68.67007 

Crame-von Mises test: not calculated 

Anderson-Darling statistic:  397.3218 

Anderson-Darling test: not calculated 



> XZ<-fitdist(LT,"nbinom")

Предупреждение

In dnbinom_mu(x, size, mu, log) : созданы NaN



> summary(XZ)

Fitting of the distribution ' nbinom ' by maximum likelihood 

Parameters : 

         estimate Std. Error

size 1.037875e+06 8.85828908

mu   7.199210e+00 0.01175151

Loglikelihood:  -123149.4   AIC:  246302.8   BIC:  246320.5 

Correlation matrix:

              size            mu

size  1.000000e+00 -1.325475e-06

mu   -1.325475e-06  1.000000e+00



> gofstat(XZ,print.test=TRUE)

Chi-squared statistic:  445.6481 

Degree of freedom of the Chi-squared distribution:  10 

Chi-squared p-value:  1.770972e-89 



> XZ<-fitdist(LT,"geom")

Предупреждения

1: In dgeom(x, prob, log) : созданы NaN

2: In dgeom(x, prob, log) : созданы NaN

> gofstat(XZ,print.test=TRUE)

Chi-squared statistic:  62647.84 

Degree of freedom of the Chi-squared distribution:  11 

Chi-squared p-value:  0 



> (XZ<-fitdist(LT,"weibull"))

Fitting of the distribution ' weibull ' by maximum likelihood 

Parameters:

      estimate  Std. Error

shape 2.937583 0.009692365

scale 8.075648 0.012729966



> gofstat(XZ,print.test=TRUE)

Kolmogorov-Smirnov statistic:  0.08801459 

Kolmogorov-Smirnov test:  rejected 

   The result of this test may be too conservative as it  

   assumes that the distribution parameters are known

Cramer-von Mises statistic:  65.77123 

Cramer-von Mises test:  rejected 

Anderson-Darling statistic:  400.9466 

Anderson-Darling test:  rejected

Добавлено через 3 часа 40 минут

Цитата:

Сообщение от Hogfather (Сообщение 289709)

Ну, теперь сам бог велел провести тест Колмогорова-Смирнова

Цитата:

Сообщение от Вляпалась... (Сообщение 289770)

Чистенько, аккуратненько, корректненько.

Если бы. Наврал ведь, а хоть бы кто поправил. Для дискретного распределения тест Колмогорова-Смирнова не применяется, так как его предельные распределения получены в предположении о непрерывности и случайных величин, и их законов распределения . Поэтому только Хи-квадрат, либо через метод обратного преобразования.

В общем, Колмогорова-Смирнова в данном случае не трогаем. Хотя, красивый результат вышел. То-то мне он подозрительным показался.

Вот тут меня спрашивают, а как посчитать R². Не знаю, зачем, но почему бы не посчитать. Формула есть, а заодно и MAE (среднюю абсолютную ошибку) посчитаем.

Для этого сделаем по-быстрому функцию

Код:

# Функция, вычисляющая R.Sqv и MAE

# (c) Hogfather, 2012

MyInfo<-function(DF,lambda,debug=F){

MyEcdf<-ecdf(DF)

MyLen<-length(DF)

MyKnots<-1:max(knots(MyEcdf))

dfecdf <- data.frame(knots=MyKnots,Fn=MyEcdf(MyKnots))

dfecdf$Fa<-ppois(dfecdf$knots, lambda)

dfecdf$R2<-(dfecdf$Fn-dfecdf$Fa)^2

TSS<-sum(dfecdf$R2)

dfecdf$RR2<-(dfecdf$Fn-mean(dfecdf$Fn))^2

ESS<-sum(dfecdf$RR2)

R2<-1-TSS/ESS

dfecdf$Err<-dfecdf$Fn-dfecdf$Fa

MAE<-mean(abs(dfecdf$Err))*MyLen

print(data.frame(R.Sqv=R2,MAE))

if(debug) print(dfecdf)

plot(dfecdf$knots,dfecdf$Err*MyLen,col="red",xlab="Число букв в слове",ylab="Ошибка аппроксимации, слов",main="Ошибки аппроксимации")

abline(0,0)

}

Скопируем в R, запустим. Дальше достаточно натравить её на наши данные и получить не только результат, но и красивый график.

LT - у нас определено выше, 7.199839 - это полученная в результате лямбда.
Результат:

Код:

> MyInfo(LT,7.199839)

     R.Sqv      MAE

1 0.999686 191.5201

R²=0.999686
MAE=191.5201 слов. Вот тут уже именно слов ;).
График
http://aspirantura.spb.ru/forum/pict...&pictureid=985

Теперь о R². Обратите внимание, что будет если мы чуть изменим лямбду.

Код:

> MyInfo(LT,8)

      R.Sqv      MAE

1 0.9758058 1950.267

R²=0.9758058, т.е. вполне годный. А вот MAE увеличилось на порядок (!). Такие дела.

Добавлено через 58 минут
Можно и совсем облениться, если данных много надо обработать, а вводить команды одни те же лень. Пишем функцию, которая делает за нас все.

Код:

# Функция, которая только за пивом не бегает

# (c) Hogfather, 2012

MyInfoPois<-function(DF){

#Подключим библиотеку

require(fitdistrplus)



# Для начала построим красивый график

descdist(DF,discrete = TRUE)

par(ask=T)

DFPois<-fitdist(DF, "pois")

lambda<-DFPois$estimate[[1]]

print(summary(DFPois))

gofstat(DFPois,print.test=TRUE)

plot(DFPois)

# А это уже было ранее. См функцию MyInfo

MyEcdf<-ecdf(DF)

MyLen<-length(DF)

MyKnots<-1:max(knots(MyEcdf))

dfecdf <- data.frame(knots=MyKnots,Fn=MyEcdf(MyKnots))

dfecdf$Fa<-ppois(dfecdf$knots, lambda)

dfecdf$R2<-(dfecdf$Fn-dfecdf$Fa)^2

TSS<-sum(dfecdf$R2)

dfecdf$RR2<-(dfecdf$Fn-mean(dfecdf$Fn))^2

ESS<-sum(dfecdf$RR2)

R2<-1-TSS/ESS

dfecdf$Err<-dfecdf$Fn-dfecdf$Fa

MAE<-mean(abs(dfecdf$Err))*MyLen

print(data.frame(R.Sqv=R2,MAE))

plot(dfecdf$knots,dfecdf$Err*MyLen,col="red",xlab="Число букв в слове",ylab="Ошибка аппроксимации, слов",main="Ошибки аппроксимации")

par(ask=F)

abline(0,0)

}

Результат запуска.

Код:

> MyInfoPois(LT)

summary statistics

------

min:  1   max:  22 

median:  7 

mean:  7.199839 

estimated sd:  2.628829 

estimated skewness:  0.519882 

estimated kurtosis:  3.143716 

Fitting of the distribution ' pois ' by maximum likelihood 

Parameters : 

       estimate Std. Error

lambda 7.199839 0.01175249

Loglikelihood:  -123149.4   AIC:  246300.8   BIC:  246309.7 

Chi-squared statistic:  445.9628 

Degree of freedom of the Chi-squared distribution:  11 

Chi-squared p-value:  1.041315e-88 

Ожидаю подтверждения смены страницы...

     R.Sqv      MAE

1 0.999686 191.5198

Ожидаю подтверждения смены страницы...

Ожидание смены страницы, чтобы можно было сохранить график. Для перехода к следующему графику, надо кликнуть по нему мышкой. Несложно, конечно, сразу выводить его в нужный файловый формат, чуть допилить функцию и всё. Как выводить в файл я уже писал.
Картинки повторять не буду. Они все уже приведены.

Лирическое отступление для Дмитрия В. и не только.

Ежу понятно, что вышеописанное никому не нужно, разве что, продемонстрировать возможности R (я себе такую цель ставил). Вообще, прежде чем проводить научное исследование, надо себе поставить цель. Поговорим об этом. У нас есть некие эмпирические данные, в данном случае соответствие длины слов количеству букв. Какие возможны варианты.
1. Нам интересна математическая модель, которая показывает зависимость количества слов в языке данной длины в данном словаре от количества букв в слове. Звучит идиотски, согласитесь.
Во всяком случае, это легко аппроксимируется полиномом или так любимой Дмитрием гаммой (но полином лучше будет). Да, в данном случае мы можем говорить о R квадрате.
Но! В данном случае у нас данные фиксированы. Нельзя добавлять или убавлять слова, поскольку это рушит нашу модель. Случайная выборка из словаря рушит всё напрочь! И модель не выполняет своей функции -- не объясняет закономерность.

2. Нам интересна закономерность, описывающая частотное распределение слов по длине. Тогда мы говорим о дискретном стохастическом процессе, причем нас интересуют именно вероятности и мы подгоняем не только дифференциальную, но и интегральную функцию распределения. Тогда ошибки считать -- заниматься профанацией. Для каждой выборки они будут свои. Задача стоит выбрать лучшее из возможных плохих вариантов. Тут в нас начинают работать информационные критерии AIC и BIC и мы выбираем из нескольких распределений лучшее. Если бы сошелся Хи-квадрат, было бы вообще счастье. Но, к сожалению, счастье бывает только в учебниках. В жизни приходится мучатся. Где-то так.
Никто, правда нам не мешает сказать, что для всего словаря Эр. квадрат такой-то, а средняя абсолютная ошибка такая-то (причем для дифференциальной и интегральной функции они будут разные, гы). А смысл?
Другой вариант, бутстреппинг. Т.е. случайная выборка, пересчет Эр квадрат и ошибки для каждого случая и отображение этого на двумерном графике. Но это чересчур брутально.

Надеюсь, что несильно наврал, а если и наврал меня поправят.