|
11.11.2012, 22:27 | #51 |
Gold Member
Регистрация: 08.04.2012
Адрес: Воронеж
Сообщений: 2,046
|
Спасибо.
Это точность аппроксимации, или я неправильно понял? Для массива в > 50 000 слов это так, меньше письки таракана *бурчит под нос* К освоению R приступить, о выполнении доложить себе лично! |
---------
Грамотей-опричникъ
Сварщик я не настоящий, а сюда просто пописать зашел |
|
Реклама | |
|
11.11.2012, 22:43 | #52 |
Platinum Member
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,304
|
|
---------
DNF is not an option
|
|
11.11.2012, 22:47 | #53 |
Gold Member
Регистрация: 08.04.2012
Адрес: Воронеж
Сообщений: 2,046
|
|
---------
Грамотей-опричникъ
Сварщик я не настоящий, а сюда просто пописать зашел |
|
11.11.2012, 23:05 | #54 |
Platinum Member
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,304
|
Ошибка модели большая вышла. Это не слов, а процентов. С чего я про слова решил?
Да и считать здесь MAPE некорректно. Убрал. |
---------
DNF is not an option
|
|
11.11.2012, 23:39 | #55 |
Gold Member
Регистрация: 08.04.2012
Адрес: Воронеж
Сообщений: 2,046
|
Hogfather, понятно.
|
---------
Грамотей-опричникъ
Сварщик я не настоящий, а сюда просто пописать зашел |
|
12.11.2012, 13:55 | #56 |
Platinum Member
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,304
|
"Нет, всё-таки откопаем..." Итак, более корректная и интересная первая подгонка, поскольку во втором случае мы просто подгоняем функцию по 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; 12.11.2012 в 22:58. |
---------
DNF is not an option
|
|
12.11.2012, 13:57 | #57 |
Silver Member
Регистрация: 31.08.2012
Адрес: Туда, вверх и налево
Сообщений: 712
|
Вот так уважаемые люди дурят провинциальных дурочек
|
---------
и чо я, дура, научнику поверила...
|
|
12.11.2012, 23:35 | #58 |
Platinum Member
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,304
|
Дык, я второй раз накалываюсь, доверяя русским публикациям. Думаете, это я сам придумал? Вот так и учусь на ошибках. Тест Шапиро-Уилка, например, даже в учебниках описан с ошибками. Заказывал через доброго человека (не указывая пальцем на 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; 13.11.2012 в 01:33. |
---------
DNF is not an option
|
|
12.11.2012, 23:42 | #59 |
Gold Member
Регистрация: 08.04.2012
Адрес: Воронеж
Сообщений: 2,046
|
Hogfather, , выходящие за фрейм.
|
---------
Грамотей-опричникъ
Сварщик я не настоящий, а сюда просто пописать зашел |
|
13.11.2012, 14:11 | #60 |
Platinum Member
Регистрация: 22.07.2010
Адрес: Санкт-Петербург
Сообщений: 3,304
|
Вот тут меня спрашивают, а как посчитать R2. Не знаю, зачем, но почему бы не посчитать. Формула есть, а заодно и 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) } LT - у нас определено выше, 7.199839 - это полученная в результате лямбда. Результат: Код:
> MyInfo(LT,7.199839) R.Sqv MAE 1 0.999686 191.5201 MAE=191.5201 слов. Вот тут уже именно слов . График Теперь о R2. Обратите внимание, что будет если мы чуть изменим лямбду. Код:
> MyInfo(LT,8) R.Sqv MAE 1 0.9758058 1950.267 Добавлено через 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 и мы выбираем из нескольких распределений лучшее. Если бы сошелся Хи-квадрат, было бы вообще счастье. Но, к сожалению, счастье бывает только в учебниках. В жизни приходится мучатся. Где-то так. Никто, правда нам не мешает сказать, что для всего словаря Эр. квадрат такой-то, а средняя абсолютная ошибка такая-то (причем для дифференциальной и интегральной функции они будут разные, гы). А смысл? Другой вариант, бутстреппинг. Т.е. случайная выборка, пересчет Эр квадрат и ошибки для каждого случая и отображение этого на двумерном графике. Но это чересчур брутально. Надеюсь, что несильно наврал, а если и наврал меня поправят. |
---------
DNF is not an option
|
|