Параметрический анализ

[Дмитрий В.]

Спасибо.

Сообщение от Hogfather

коэффициент детерминации практически единицу.

Это точность аппроксимации, или я неправильно понял?

Сообщение от Hogfather

Средняя абсолютная ошибка модели составляет 181 слово.

Для массива в > 50 000 слов это так, меньше письки таракана
*бурчит под нос* К освоению R приступить, о выполнении доложить себе лично!

[Hogfather]

Сообщение от Дмитрий В.

Это точность аппроксимации, или я неправильно понял?

Ну, грубо говоря, да.

[Дмитрий В.]

Сообщение от Hogfather

Ну, грубо говоря, да.

А, вот, он самый! Правда, R² мне почему-то точностью аппроксимации называли

[Hogfather]

Ошибка модели большая вышла. Это не слов, а процентов. С чего я про слова решил?
Да и считать здесь MAPE некорректно. Убрал.

[Дмитрий В.]

Hogfather, понятно.

[Hogfather]

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

Итак, более корректная и интересная первая подгонка, поскольку во втором случае мы просто подгоняем функцию по 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)

Рисунок 1 -- Подгонка гамма-распределения


Рисунок 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

Вот такая красота.


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

> 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

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

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

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

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

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

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

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

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

[Вляпалась...]

Сообщение от Hogfather

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

Вот так уважаемые люди дурят провинциальных дурочек

[Hogfather]

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

Вот так уважаемые люди дурят провинциальных дурочек

Дык, я второй раз накалываюсь, доверяя русским публикациям. Думаете, это я сам придумал? Вот так и учусь на ошибках. Тест Шапиро-Уилка, например, даже в учебниках описан с ошибками. Заказывал через доброго человека (не указывая пальцем на watteau) оригинал статьи, чтобы разобраться.

Надеюсь, что сюда математик нормальный забредет, подскажет идею какую-нибудь.

Добавлено через 8 часов 12 минут
Мысль для Дмитрий В. и не только.
Пусть у Вас имеются данные по 5 языкам. Чтобы на одном графике показать все распределения, лучше всего использовать диаграмму "Ящик с усами".Привожу пример. Данные выдуманы (заполняются в начале скрипта). У Дмитрия наверняка они есть.

#Инициируем данные, чтобы были, заполнив по 1000 чисел из распределения Пуассона и нормального распределения.

lgRU<-rpois(1000,5)
lgGB<-rpois(1000,5.1)
lgNL<-rpois(1000,7)
lgFR<-rnorm(1000,2,1)
lgDE<-rnorm(1000,3,1)

##### Всё счастье тут. Вот он, график!
boxplot(list("Рус"=lgRU,"Анг"=lgGB,"Чук"=lgNL,"Нен"=lgDE,"Франц"=lgFR),main="Сравнение распределения\n длины слов в языках",xlab="Язык",ylab="Длина слова",col = "lavender", notch = TRUE, varwidth = TRUE)

Результат:


Добавлено через 58 минут
А вот пример Бутстреппинга

> bw<-bootdist(XY,niter=1001)
> plot(bw)
> summary(bw)
Parametric bootstrap medians and 95% percentile CI 
  Median     2.5%    97.5% 
7.199874 7.175821 7.223915

[Дмитрий В.]

Hogfather, , выходящие за фрейм.

[Hogfather]

Вот тут меня спрашивают, а как посчитать 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 слов. Вот тут уже именно слов .
График


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

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