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

Портал аспирантов (http://www.aspirantura.spb.ru/forum/index.php)

- Филологические науки (http://www.aspirantura.spb.ru/forum/forumdisplay.php?f=136)

- - Параметрический анализ (http://www.aspirantura.spb.ru/forum/showthread.php?t=10490)

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

Итак, коллеги. Товарищ Дмитрий В. получил интересные результаты, стал мучать их в 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)

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

Что мы, собственно говоря видим. А видим, что распределение у нас вполне милое, да слегка несимметричное, но с кем не бывает.
Пытаемся натянуть сову на глобус. Для этого используем подгонку распределения методом максимального правдоподобия (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)

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

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

Код:

> 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)

Результат
http://aspirantura.spb.ru/forum/pict...&pictureid=972

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

Код:

> 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 (Сообщение 289709)

Дальше мы просто обязаны поффтыкать на картинки, иначе мы не ученые, а кот начхал.

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

Строим четыре графика для наших распределений.

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

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

Ай, да Hogfather, хочется сказать, ай, молодец!

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

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

Это 5. Даже не 5, а 7-8, где-то так. Теперь на Вам.. ой на Вас... ну неважно, за Вас... просто обязаны выйти замуж. Как честный человек :) Но соавтором - это как минимум.

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

Hogfather, выскажу огроменное спасибо и тут.
Вляпалась...,

Цитата:

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

Но соавтором - это как минимум.

Обязательно, постараемся.

Цитата:

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

Ай, да Hogfather, хочется сказать, ай, молодец!

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

(аплодисменты, подкрепленные топотом ног) Ай, да Hogfather! Ай, да молодец!!

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

Как говорится. не только методом максимального правдоподобия славен R. Ту же задачу можно попробовать решить нелинейным методом наименьших квадратов. Для этого построим кумулятивную (интегральную) функцию распределения и попробуем подогнать понравившегося нам Пуассона. В общем, сделаем примерно то, что пытался проделать Дмитрий В. в Excel.

Код:

> # Понеслась!

> # Строим кумулятивную функцию

> MyEcdf<-ecdf(LT)

># Делаем таблицу (фрейм) для аппроксимации

># Обратите внимание, поскольку я все взял в скобки, результат отображается сразу на экране

> (dfecdf <- data.frame(knots=knots(MyEcdf),Fn=MyEcdf(1:22)))

   knots           Fn

1      1 0.0001726552

2      2 0.0052947609

3      3 0.0598346346

4      4 0.1643869780

5      5 0.2903102039

6      6 0.4254417097

7      7 0.5696663917

8      8 0.7053158632

9      9 0.8131294723

10    10 0.8901720797

11    11 0.9389951465

12    12 0.9676559173

13    13 0.9840389817

14    14 0.9920194909

15    15 0.9961248489

16    16 0.9984652867

17    17 0.9994820343

18    18 0.9997889769

19    19 0.9999232643

20    20 0.9999616322

21    21 0.9999808161

22    22 1.0000000000



> # Строим модель



> mdl<-nls( Fn ~ ppois(knots,lambda), data=dfecdf,model=T)

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

In nls(Fn ~ ppois(knots, lambda), data = dfecdf, model = T) :

  Для некоторых параметров не указаны стартовые значения.

Инициализую ‘lambda’ до '1.'.

Укажите 'start' или я использую модель 'selfStart'



> # Информация о модели

> summary(mdl)



Formula: Fn ~ ppois(knots, lambda)



Parameters:

       Estimate Std. Error t value Pr(>|t|)    

lambda  7.16774    0.01924   372.5   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 



Residual standard error: 0.006275 on 21 degrees of freedom



Number of iterations to convergence: 5 

Achieved convergence tolerance: 3.776e-08

Посчитаем адекватность полученной модели.

Код:

> # Расчет адекватности модели

> (RSS.p <- sum(residuals(mdl)^2))

[1] 0.000826937

> (TSS <- sum((dfecdf$Fn - mean(dfecdf$Fn))^2))

[1] 2.981961



> # коэффициент детерминации

> 1 - (RSS.p/TSS)

[1] 0.9997227

Что мы имеем с гуся. А с гуся имеем чуть другую лямбду (7.16774) и коэффициент детерминации практически единицу.
Для лямбды можно посчитать доверительный интервал

Код:

> confint(mdl)

Waiting for profiling to be done...

    2.5%    97.5% 

7.127763 7.207781

Графически ошибки модели можно изобразить вот так.

Код:

> plot(residuals(mdl),main="Ошибки модели")

> abline(0,0)

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

Спасибо.

Цитата:

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

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

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

Цитата:

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

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

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

Цитата:

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

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

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

Цитата:

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

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

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

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

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)

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)

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

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

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

Цитата:

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

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

Цитата:

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

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

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

Результат:
http://aspirantura.spb.ru/forum/pict...&pictureid=983

Добавлено через 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

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

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

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

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