R:Статистическая проверка принадлежности нормальному распределения — различия между версиями
м |
м |
||
| Строка 13: | Строка 13: | ||
Для демонстрации работы функций, реализующий различные критерий проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение: | Для демонстрации работы функций, реализующий различные критерий проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение: | ||
| − | + | {{r-code|code= | |
> x <- rnorm(n = 100) | > x <- rnorm(n = 100) | ||
| − | + | }} | |
=== Статистические критерии === | === Статистические критерии === | ||
| Строка 64: | Строка 64: | ||
Пример вызова данной функции: | Пример вызова данной функции: | ||
| − | + | {{r-code|code= | |
> normalTest(x, method = "sw") | > normalTest(x, method = "sw") | ||
| Строка 78: | Строка 78: | ||
Description: | Description: | ||
Fri Feb 14 19:59:59 2014 by user: | Fri Feb 14 19:59:59 2014 by user: | ||
| − | + | }} | |
Помимо функции <code>normalTest()</code> данный пакет включает в себя следующие функции: | Помимо функции <code>normalTest()</code> данный пакет включает в себя следующие функции: | ||
| Строка 115: | Строка 115: | ||
С помощью <code>apply</code>-функций можно последовательно применить функцию к вектору, списку или массиву. Прежде чем всего нам необходимо сформировать таблицу данных. С помощью функции <code>replicate()</code> сгенерируем 10 переменных, имеющих стандартное нормальное распределение, которые объединяются в класс <code>data.frame</code>. | С помощью <code>apply</code>-функций можно последовательно применить функцию к вектору, списку или массиву. Прежде чем всего нам необходимо сформировать таблицу данных. С помощью функции <code>replicate()</code> сгенерируем 10 переменных, имеющих стандартное нормальное распределение, которые объединяются в класс <code>data.frame</code>. | ||
| − | + | {{r-code|code= | |
> DF <- data.frame(replicate(n = 10, rnorm(n = 100))) | > DF <- data.frame(replicate(n = 10, rnorm(n = 100))) | ||
| − | + | }} | |
Структура сгенерированной таблицы выглядит следующим образом: | Структура сгенерированной таблицы выглядит следующим образом: | ||
| − | + | {{r-code|code= | |
> str(DF) | > str(DF) | ||
'data.frame': 100 obs. of 10 variables: | 'data.frame': 100 obs. of 10 variables: | ||
| Строка 134: | Строка 134: | ||
$ X9 : num -0.8508 0.4481 -0.2828 -0.5464 0.0605 ... | $ X9 : num -0.8508 0.4481 -0.2828 -0.5464 0.0605 ... | ||
$ X10: num 1.421 0.408 1.254 -0.956 -1.91 ... | $ X10: num 1.421 0.408 1.254 -0.956 -1.91 ... | ||
| − | + | }} | |
Для решения поставленной задачи можно воспользоваться функцией <code>sapply()</code>. Но прежде, нам необходимо немного отформатировать формат вывода результатов нашей функции: нам нужно извлечь значения критерия и его уровень значимости, т.к. результат функции <code>shapiro.test()</code> содержит также информацию, которая не подлежит включению в таблицу. | Для решения поставленной задачи можно воспользоваться функцией <code>sapply()</code>. Но прежде, нам необходимо немного отформатировать формат вывода результатов нашей функции: нам нужно извлечь значения критерия и его уровень значимости, т.к. результат функции <code>shapiro.test()</code> содержит также информацию, которая не подлежит включению в таблицу. | ||
| − | + | {{r-code|code= | |
> shapiro.test(x) | > shapiro.test(x) | ||
| Строка 145: | Строка 145: | ||
data: x | data: x | ||
W = 0.9903, p-value = 0.6882 | W = 0.9903, p-value = 0.6882 | ||
| − | + | }} | |
Структура результата применения функции <code>shapiro.test()</code> представлена ниже: | Структура результата применения функции <code>shapiro.test()</code> представлена ниже: | ||
| − | + | {{r-code|code= | |
> str(shapiro.test(x)) | > str(shapiro.test(x)) | ||
List of 4 | List of 4 | ||
| Строка 158: | Строка 158: | ||
$ data.name: chr "x" | $ data.name: chr "x" | ||
- attr(*, "class")= chr "htest" | - attr(*, "class")= chr "htest" | ||
| − | + | }} | |
Как видим, помимо значений критерия и уровня значимости здесь содержится информация о применяемом методе. Мы можем отфильтровать вывод следующим образом: | Как видим, помимо значений критерия и уровня значимости здесь содержится информация о применяемом методе. Мы можем отфильтровать вывод следующим образом: | ||
| − | + | {{r-code|code= | |
> normTest <- function (x) { | > normTest <- function (x) { | ||
+ res <- shapiro.test(x) | + res <- shapiro.test(x) | ||
+ return(c(res$statistic, p.value = res$p.value)) | + return(c(res$statistic, p.value = res$p.value)) | ||
+ } | + } | ||
| − | + | }} | |
Результат теперь будет выглядеть следующим образом: | Результат теперь будет выглядеть следующим образом: | ||
| − | + | {{r-code|code= | |
> normTest(x) | > normTest(x) | ||
W p.value | W p.value | ||
0.9903 0.6882 | 0.9903 0.6882 | ||
| − | + | }} | |
Теперь можно использовать данную функцию при обработке столбцов нашей таблицы. | Теперь можно использовать данную функцию при обработке столбцов нашей таблицы. | ||
| − | + | {{r-code|code= | |
> t(sapply(DF, normTest)) | > t(sapply(DF, normTest)) | ||
W p.value | W p.value | ||
| Строка 192: | Строка 192: | ||
X9 0.9915 0.7834 | X9 0.9915 0.7834 | ||
X10 0.9808 0.1531 | X10 0.9808 0.1531 | ||
| − | + | }} | |
Того же результата можно добиться и с помощью функции <code>lapply()</code><ref>По результатам сравнения производительности, данный вариант оказался чуть быстрее предыдущего.</ref>: | Того же результата можно добиться и с помощью функции <code>lapply()</code><ref>По результатам сравнения производительности, данный вариант оказался чуть быстрее предыдущего.</ref>: | ||
| − | + | {{r-code|code= | |
> do.call(rbind, lapply(DF, normTest)) | > do.call(rbind, lapply(DF, normTest)) | ||
W p.value | W p.value | ||
| Строка 209: | Строка 209: | ||
X9 0.9915 0.7834 | X9 0.9915 0.7834 | ||
X10 0.9808 0.1531 | X10 0.9808 0.1531 | ||
| − | + | }} | |
===== Применение функций к нескольким по группам ===== | ===== Применение функций к нескольким по группам ===== | ||
| Строка 215: | Строка 215: | ||
Добавим к нашей таблице группы испытуемых: | Добавим к нашей таблице группы испытуемых: | ||
| − | + | {{r-code|code= | |
> DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE)) | > DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE)) | ||
| − | + | }} | |
Состав групп получился следующим: | Состав групп получился следующим: | ||
| − | + | {{r-code|code= | |
> table(DF$GRP) | > table(DF$GRP) | ||
A B C | A B C | ||
38 25 37 | 38 25 37 | ||
| − | + | }} | |
Рассчитаем значения критерия Шапиро - Уилка для первого столбца для каждоый группы испытуемых: | Рассчитаем значения критерия Шапиро - Уилка для первого столбца для каждоый группы испытуемых: | ||
| − | + | {{r-code|code= | |
> do.call(rbind, tapply(DF$X1, DF$GRP, normTest)) | > do.call(rbind, tapply(DF$X1, DF$GRP, normTest)) | ||
W p.value | W p.value | ||
| Строка 235: | Строка 235: | ||
B 0.9607 0.28697 | B 0.9607 0.28697 | ||
C 0.9410 0.07256 | C 0.9410 0.07256 | ||
| − | + | }} | |
=== Графические методы === | === Графические методы === | ||
| Строка 245: | Строка 245: | ||
Построение Q–Q plot с помощью пакета <code>stats</code> выглядит следующим образом: | Построение Q–Q plot с помощью пакета <code>stats</code> выглядит следующим образом: | ||
| − | + | {{r-code|code= | |
> qqnorm(x) | > qqnorm(x) | ||
> qqline(x) | > qqline(x) | ||
| − | + | }} | |
| + | |||
[[Файл:Stats-qqnorm.svg|400px|центр]] | [[Файл:Stats-qqnorm.svg|400px|центр]] | ||
| Строка 255: | Строка 256: | ||
Построение Q–Q plot с помощью пакета <code>QTLRel</code> выглядит следующим образом: | Построение Q–Q plot с помощью пакета <code>QTLRel</code> выглядит следующим образом: | ||
| − | + | {{r-code|code= | |
> qqPlot(x, x = "norm") | > qqPlot(x, x = "norm") | ||
| − | + | }} | |
| + | |||
[[Файл:Qtlrel-qqplot.svg|400px|центр]] | [[Файл:Qtlrel-qqplot.svg|400px|центр]] | ||
| Строка 264: | Строка 266: | ||
Альтернативный вариант реализован в функции <code>qqPlot()</code> из пакета <code>car</code>: | Альтернативный вариант реализован в функции <code>qqPlot()</code> из пакета <code>car</code>: | ||
| − | + | {{r-code|code= | |
> qqPlot(x, distribution = "norm") | > qqPlot(x, distribution = "norm") | ||
| − | + | }} | |
| + | |||
[[Файл:Car-qqplot.svg|400px|центр]] | [[Файл:Car-qqplot.svg|400px|центр]] | ||
| Строка 273: | Строка 276: | ||
Построение P-P plot можно осуществить с помощью функции <code>probplot</code> из пакета <code>e1071</code>: | Построение P-P plot можно осуществить с помощью функции <code>probplot</code> из пакета <code>e1071</code>: | ||
| − | + | {{r-code|code= | |
> probplot(x, qdist = qnorm) | > probplot(x, qdist = qnorm) | ||
| − | + | }} | |
| + | |||
[[Файл:E1071-probplot.svg|400px|центр]] | [[Файл:E1071-probplot.svg|400px|центр]] | ||
| Строка 282: | Строка 286: | ||
Ещё один интересный способ графического анализа представлен функцией <code>histDist</code> из пакета <code>gamlss</code>: | Ещё один интересный способ графического анализа представлен функцией <code>histDist</code> из пакета <code>gamlss</code>: | ||
| − | + | {{r-code|code= | |
> histDist(x, family = "NO", density = TRUE) | > histDist(x, family = "NO", density = TRUE) | ||
| Строка 299: | Строка 303: | ||
AIC: 307.414 | AIC: 307.414 | ||
SBC: 312.624 | SBC: 312.624 | ||
| − | + | }} | |
| + | |||
[[Файл:gamlss-histdist.svg|400px|центр]] | [[Файл:gamlss-histdist.svg|400px|центр]] | ||
Версия 16:24, 17 февраля 2014
 
|
Материал «R:Статистическая проверка принадлежности нормальному распределения», созданный автором Артём Клевцов, публикуется на условиях лицензии Creative Commons «Attribution» («Атрибуция») 4.0 Всемирная. | |
|
Предположение о принадлежности случайной величины нормальному закону распределения лежит в основе многих статистических методов и критериев. В ряде случаев соблюдение данного требования является критичным для применения того или иного метода.
На практике мы встречаемся с двумя вариантами задач по проверке принадлежности распределения нормальному закону: для одномерного и многомерного распределения.
 |
Перед использованием функций из пакетов их необходимо предварительно установить и загрузить: КодR <syntaxhighlight lang="r">> install.packages(pkgs = "pkgname") > library(package = "pkgname")</syntaxhighlight> |
Содержание
Одномерное нормальное распределение
В качестве [math]H_0[/math] для всех нижеприведённых критериев является предположение, что «случайная величина [math]X[/math] распределена нормально».
Для демонстрации работы функций, реализующий различные критерий проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение:
<syntaxhighlight lang="r">> x <- rnorm(n = 100)</syntaxhighlight>
Статистические критерии
Пакет stats
В данном пакете реализованы две функции, которые позволяют осуществить проверку принадлежности распределения нормальному закону.
-
shapiro.test- критерий Шапиро - Уилка -
ks.test- критерий Колмогорова - Смирнова[1]
Данные функции возвращают результат в виде S3-класса - htest.
Пакет nortest
В данный пакет входят следующие функции:
-
ad.test- критерий Андерсона - Дарлинга -
cvm.test- критерий Крамера - фон Мизеса -
lillie.test- критерий Лиллиефорса -
pearson.test- критерий [math]\chi^2[/math] Пирсона -
sf.test- критерий Шапиро - Франчия
Данные функции возвращают результат в виде S3-класса - htest.
Пакет moments
В данный пакет входят следующие функции:
-
agostino.test- критерий Д'Агостино -
bonett.test- критерий Бонетта – Сайера -
jarque.test- критерий Жарка-Бера
Данные функции также возвращают результат в виде S3-класса - htest.
Пакет fBasics
В данном пакете не предлагается никакой оригинальной реализации критериев - код в основном заимствован из пакетов stats, nortest, moments. Данный пакет предлагает альтернативный вывод результатов в виде объекта S4-класса fHTEST, в том время как все предыдущие функции использовали S3-класс htest.
Функция normalTest() является «обёрктой» для ряда функций из того же пакета - fBasics. Задать необходимый критерий можно задать с помощью аргумента method. Доступны следующие критерии:
-
sw- критерий Шапиро - Уилка -
jb- критерий Жарка-Бера -
ks- критерий Колмогорова - Смирнова -
da- критерий Д'Агостино -
ad- критерий Андерсона - Дарлинга.
Пример вызова данной функции:
<syntaxhighlight lang="r">> normalTest(x, method = "sw")
Title:
Shapiro - Wilk Normality Test
Test Results:
STATISTIC: W: 0.9831 P VALUE: 0.2301
Description:
Fri Feb 14 19:59:59 2014 by user:</syntaxhighlight>
Помимо функции normalTest() данный пакет включает в себя следующие функции:
-
shapiroTest- критерий Шапиро - Уилка -
ksnormTest- критерий Колмогорова - Смирнова[2] -
jarqueberaTest- критерий Жарка-Бера -
dagoTest- критерий Д'Агостино -
adTest- критерий Андерсона - Дарлинга -
cvmTest- критерий Крамера - фон Мизеса -
lillieTest- критерий Лиллиефорса -
pchiTest- критерий Пирсона -
sfTest- критерий Шапиро - Франчия
Данные функции также возвращают результат в виде S4-класса - fHTEST.
Пакет TeachingDemos
Данные пакет содержит только одну функцию, имеющую отношение к критериям проверки принадлежности распределения нормальному закону - SnowsPenultimateNormalityTest(). Данная функция возвращают результат в виде S3-класса - htest.
Пакет tseries
Данный пакет содержит только одну функцию, имеющую отношение к критериям проверки принадлежности распределения нормальному закону - jarque.bera.test, которая является реализацией критерия Жарка-Бера. Данная функция возвращают результат в виде S3-класса - htest.
Пакет lawstat
В данный пакет входят следующие функции:
-
rjb.test- критерий Жарка-Бера -
sj.test- SJ-критерий
Маленькие хитрости
Применение функций к нескольким переменным
С помощью apply-функций можно последовательно применить функцию к вектору, списку или массиву. Прежде чем всего нам необходимо сформировать таблицу данных. С помощью функции replicate() сгенерируем 10 переменных, имеющих стандартное нормальное распределение, которые объединяются в класс data.frame.
<syntaxhighlight lang="r">> DF <- data.frame(replicate(n = 10, rnorm(n = 100)))</syntaxhighlight>
Структура сгенерированной таблицы выглядит следующим образом:
<syntaxhighlight lang="r">> str(DF) 'data.frame': 100 obs. of 10 variables:
$ X1 : num 1.051 1.08 -0.477 -1.396 3.423 ... $ X2 : num -0.602 2.29 -0.758 -1.615 -0.364 ... $ X3 : num 0.0559 -1.0117 0.5242 0.4105 -0.3191 ... $ X4 : num -0.0965 0.2006 0.29 0.7702 -0.0182 ... $ X5 : num -0.7074 -1.6111 0.3478 0.2504 0.0609 ... $ X6 : num -1.432 0.535 -0.932 0.581 -1.606 ... $ X7 : num -1.42407 -0.31827 -2.04648 -0.19856 0.00301 ... $ X8 : num 0.511 0.192 0.467 -1.308 2.496 ... $ X9 : num -0.8508 0.4481 -0.2828 -0.5464 0.0605 ... $ X10: num 1.421 0.408 1.254 -0.956 -1.91 ...</syntaxhighlight>
Для решения поставленной задачи можно воспользоваться функцией sapply(). Но прежде, нам необходимо немного отформатировать формат вывода результатов нашей функции: нам нужно извлечь значения критерия и его уровень значимости, т.к. результат функции shapiro.test() содержит также информацию, которая не подлежит включению в таблицу.
<syntaxhighlight lang="r">> shapiro.test(x)
Shapiro-Wilk normality test
data: x W = 0.9903, p-value = 0.6882</syntaxhighlight>
Структура результата применения функции shapiro.test() представлена ниже:
<syntaxhighlight lang="r">> str(shapiro.test(x)) List of 4
$ statistic: Named num 0.99 ..- attr(*, "names")= chr "W" $ p.value : num 0.688 $ method : chr "Shapiro-Wilk normality test" $ data.name: chr "x" - attr(*, "class")= chr "htest"</syntaxhighlight>
Как видим, помимо значений критерия и уровня значимости здесь содержится информация о применяемом методе. Мы можем отфильтровать вывод следующим образом:
<syntaxhighlight lang="r">> normTest <- function (x) { + res <- shapiro.test(x) + return(c(res$statistic, p.value = res$p.value)) + }</syntaxhighlight>
Результат теперь будет выглядеть следующим образом:
<syntaxhighlight lang="r">> normTest(x)
W p.value 0.9903 0.6882</syntaxhighlight>
Теперь можно использовать данную функцию при обработке столбцов нашей таблицы.
<syntaxhighlight lang="r">> t(sapply(DF, normTest))
W p.value
X1 0.9831 0.2301 X2 0.9936 0.9213 X3 0.9800 0.1333 X4 0.9829 0.2219 X5 0.9874 0.4625 X6 0.9862 0.3874 X7 0.9839 0.2617 X8 0.9833 0.2360 X9 0.9915 0.7834 X10 0.9808 0.1531</syntaxhighlight>
Того же результата можно добиться и с помощью функции lapply()[3]:
<syntaxhighlight lang="r">> do.call(rbind, lapply(DF, normTest))
W p.value
X1 0.9831 0.2301 X2 0.9936 0.9213 X3 0.9800 0.1333 X4 0.9829 0.2219 X5 0.9874 0.4625 X6 0.9862 0.3874 X7 0.9839 0.2617 X8 0.9833 0.2360 X9 0.9915 0.7834 X10 0.9808 0.1531</syntaxhighlight>
Применение функций к нескольким по группам
Добавим к нашей таблице группы испытуемых:
<syntaxhighlight lang="r">> DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE))</syntaxhighlight>
Состав групп получился следующим:
<syntaxhighlight lang="r">> table(DF$GRP)
A B C
38 25 37</syntaxhighlight>
Рассчитаем значения критерия Шапиро - Уилка для первого столбца для каждоый группы испытуемых:
<syntaxhighlight lang="r">> do.call(rbind, tapply(DF$X1, DF$GRP, normTest))
W p.value
A 0.9522 0.13281 B 0.9607 0.28697 C 0.9410 0.07256</syntaxhighlight>
Графические методы
Многие исследователи также используют графические методы для определения степени отклонения распределения от нормального закона. В R реализована возможность построения Q-Q и P-P графиков, гистограмм и кривых распределения плотности вероятностей.
Пакет stats
Построение Q–Q plot с помощью пакета stats выглядит следующим образом:
<syntaxhighlight lang="r">> qqnorm(x) > qqline(x)</syntaxhighlight>
Пакет QTLRel
Построение Q–Q plot с помощью пакета QTLRel выглядит следующим образом:
<syntaxhighlight lang="r">> qqPlot(x, x = "norm")</syntaxhighlight>
Пакет car
Альтернативный вариант реализован в функции qqPlot() из пакета car:
<syntaxhighlight lang="r">> qqPlot(x, distribution = "norm")</syntaxhighlight>
Пакет e1071
Построение P-P plot можно осуществить с помощью функции probplot из пакета e1071:
<syntaxhighlight lang="r">> probplot(x, qdist = qnorm)</syntaxhighlight>
Пакет gamlss
Ещё один интересный способ графического анализа представлен функцией histDist из пакета gamlss:
<syntaxhighlight lang="r">> histDist(x, family = "NO", density = TRUE)
Family: c("NO", "Normal") Fitting method: "nlminb"
Call: gamlssML(y = y, family = "NO", formula = x)
Mu Coefficients: [1] -0.2273 Sigma Coefficients: [1] 0.09813
Degrees of Freedom for the fit: 2 Residual Deg. of Freedom 98
Global Deviance: 303.414
AIC: 307.414
SBC: 312.624</syntaxhighlight>
С помощью аргумента family можно задать семейство распределений для подгонки и сравнения[4].
Многомерное нормальное распределение
Примечания
- ↑ Для оценки нормальности вызов выглядит следующим образом:
ks.test(x, y = "pnorm") - ↑ Данная функция вызывает
ks.test(x, "pnorm")для трёх альтернативных гипотез - двусторонней и двух односторонних. - ↑ По результатам сравнения производительности, данный вариант оказался чуть быстрее предыдущего.
- ↑ Более подробную информацию о доступных семействах распределений можно получить с помощью команды
help("gamlss.family").
