Выборка, выборочное распределение одномерной случайной величины. Построение эмпирической плотности и эмпирической функции распределения в пакете R.
Выборка, выборочное распределение одномерной случайной величины. Построение эмпирической плотности и эмпирической функции распределения в пакете R.
Теория:
Эмпирическую плотность можно оценить с помощью гистограммы, ядерных оценок.
K(x) – четная, интегрируемая:
Пусть h > 0, тогда:
(x) =
где , i = 1:n – значения выборки
Рассмотрим некоторые ядра:
Практика в R:
density(x, bw = "nrd0", adjust = 1,
kernel = c("gaussian", "epanechnikov", "rectangular",
"triangular", "biweight",
"cosine", "optcosine"),
weights = NULL, window = kernel, width,
give.Rkern = FALSE,
n = 512, from, to, cut = 3, na.rm = FALSE, ...)
Обязательные параметры: выборка, ядро, используя параметр kernel. Так же можно указать критерий по параметру сглаживания или его численное значение - параметр bw = "nrd0", (другие. так же указывать как строковую константу. nrd, ucv, bcv, SJ).
ecdf(x)
На вход только выборку.
Также добавив приставку d – плотность или p – функция распределения, можно получить соответственно плотность или функцию распределения указанного после приставки закона распределения:
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
dunif(x, min = 0, max = 1, log = FALSE)
Выборочные числовые характеристики случайных величин. Оценка основных выборочных числовых характеристик в пакете R.
Теория:
Практика в R:
func.mean <- function() {
#Выборочное среднее для вектора
mean(trees$Girth)
[1] 13.24839
#Выборочное среднее для кадра данных
mean(trees)
#При отсутствующих значениях (na.rm=FALSE - NA, na.rm=TRUE - пропустит NA)
mean(c(1,2,NA))
[1] NA
mean(c(1,2,NA),na.rm=TRUE)
[1] 1.5
}
func.mediana <- function() {
#Выборочная медиана
median(trees$Girth)
}
func.quantile <- function() {
#Квантиль уровней 0, 25, 50, 75, 100 процентов - вероятность того, что нет значений выше
quantile(trees$Girth)
0% 25% 50% 75% 100%
8.30 11.05 12.90 15.25 20.60
#Квантиль уровней seq(0,1,0.1)
quantile(trees$Girth, seq(0,1,0.1))
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
8.3 10.5 11.0 11.2 11.4 12.9 13.7 14.2 16.3 17.9 20.6
}
func.var <- function() {
#Вычисление несмещенной выборочной дисперсии (тупо делится на n-1, а не на n) для кадра данных
var(trees)
#Для вектора
var(trees$Girth)
}
func.sqrt_disp <- function() {
#Вычисление среднеквадратичного отклонения
sd(trees$Girth)
}
func.diff <- function() {
#Находим min и max значения с помощью range и находим размах, используя diff
min_max = range(trees$Girth)
diff(min_max)
}
func.IQR <- function() {
#Выборочный междукваРтильный размах (между 1/4 и 3/4 кваНтилями)
IQR(trees$Girth)
}
Типы статистических данных. Частотные распределения. Маргинальные частотные распределения. Проблема нахождения выборочных числовых характеристик по группированным данным. Ее решение в пакете R.
Теория:
Практика в R:
Таблица частот cтроится функцией:
table(x,…)
На вход выборка или несколько выборок если хотим найти многомерное частотное распределение.
cut(x,br)
На входе выборка и массив из точек, которыми вся область значений выборки разбивается на интервалы. Функция возвращает не список x1,…,xN – т.е. выборку, а список интервалов, которым принадлежат x1,…,xN соответственно
x <- c(0,1,2,3,4,5,1,2,3,4,12,12,31,23,12,14,13,1,1,1,2,3)
range(x)
[1] 0 31
br <- seq(-1,31,4)
x.cut <- cut(x,br)
x.cut
[1] (-1,3] (-1,3] (-1,3] (-1,3] (3,7] (3,7] (-1,3] (-1,3] (-1,3] (3,7]
[11] (11,15] (11,15] (27,31] (19,23] (11,15] (11,15] (11,15] (-1,3] (-1,3] (-1,3]
[21] (-1,3] (-1,3]
Levels: (-1,3] (3,7] (7,11] (11,15] (15,19] (19,23] (23,27] (27,31]
table(x.cut)
(-1,3] (3,7] (7,11] (11,15] (15,19] (19,23] (23,27] (27,31]
12 3 0 5 0 1 0 1
y<-c('a','v','a','s','s')
table(y)
y
A s v
2 2 1
Маргинальные частоты – для многомерных частотных распределений можно посчитать частоту по какому-то одному признаку или по нескольким. Пример:
attach(iris)
i.l <-cut(Sepal.Length,breaks=6)
i.w <-cut(Sepal.Width,breaks=6)
table(i.l,i.w)
I.w
i.l (2,2.4] (2.4,2.8] (2.8,3.2] (3.2,3.6] (3.6,4] (4,4.4]
(4.3,4.9] 2 1 14 5 0 0
(4.9,5.5] 5 4 3 15 8 2
(5.5,6.1] 2 17 13 1 2 1
(6.1,6.7] 2 10 17 6 0 0
(6.7,7.3] 0 1 11 1 0 0
(7.3,7.9] 0 3 2 0 2 0
table.i <-table(i.l,i.w)
margin.table(table.i,1)
I.l
(4.3,4.9] (4.9,5.5] (5.5,6.1] (6.1,6.7] (6.7,7.3] (7.3,7.9]
22 37 36 35 13 7
margin.table(table.i,2)
I.w
(2,2.4] (2.4,2.8] (2.8,3.2] (3.2,3.6] (3.6,4] (4,4.4]
11 36 60 28 12 3
Сначала построили таблицу частот по двум параметрам, потом посчитали маргинальные частоты по первому и второму признакам. По первому – просуммировали строки, по второму столбцы
Далее для трех признаков, и маргинальные частоты по 2-м из них
table(Species,i.l,i.w)->tb.i.sp
ftable(tb.i.sp)
i.w (2,2.4] (2.4,2.8] (2.8,3.2] (3.2,3.6] (3.6,4] (4,4.4]
Species i.l
setosa (4.3,4.9] 1 0 14 5 0 0
(4.9,5.5] 0 0 2 15 8 2
(5.5,6.1] 0 0 0 0 2 1
(6.1,6.7] 0 0 0 0 0 0
(6.7,7.3] 0 0 0 0 0 0
(7.3,7.9] 0 0 0 0 0 0
versicolor (4.3,4.9] 1 0 0 0 0 0
(4.9,5.5] 5 4 1 0 0 0
(5.5,6.1] 1 11 10 1 0 0
(6.1,6.7] 2 2 8 1 0 0
(6.7,7.3] 0 1 2 0 0 0
(7.3,7.9] 0 0 0 0 0 0
virginica (4.3,4.9] 0 1 0 0 0 0
(4.9,5.5] 0 0 0 0 0 0
(5.5,6.1] 1 6 3 0 0 0
(6.1,6.7] 0 8 9 5 0 0
(6.7,7.3] 0 0 9 1 0 0
(7.3,7.9] 0 3 2 0 2 0
margin.table(tb.i.sp,c(1,3))
I.w
Species (2,2.4] (2.4,2.8] (2.8,3.2] (3.2,3.6] (3.6,4] (4,4.4]
Setosa 1 0 16 20 10 3
Versicolor 9 18 21 2 0 0
Virginica 1 18 23 6 2 0
Проблема нахождения выборочных числовых характеристик по группированным данным.
То есть у нас вместо выборки таблица частот. Для поиска числовых характеристик требуются другие формулы они просты, но получается нужно переписывать все функции R для работы с выборками под группированные данные. Вместо этого создатели языка предложили следующее решение: по группированным данным, генерировать фиктивную выборку.
a<-rbinom(10,10,0.2)
mean(a)
[1] 1.7table(a)
a
0 1 2 3
1 3 4 2
rep(c(0,1,2,3),c(1,3,4,2))->a.fake
mean(a.fake)
[1] 1.7a.fake
[1] 0 1 1 1 2 2 2 2 3 3a
[1] 1 1 3 0 2 2 2 2 3 1Метод обратной функции
Пусть – функция распределения, число из которого мы хотим получить, – обратная функция к (она существует на интервале ) строгой монотонности ), – случайная величина, распределённая равномерно в . Заметим, что:
То есть – случайная величина, распределённая по нужному нам закону распределения.
Метод отсечения
Пусть – плотность нужного нам распределения, а – плотность какого-то популярного распределения, причём . Тогда алгоритм генерации такой:
Генерируем из равномерного на распределения и из распределения
Если , то возвращаемся на шаг 1. Иначе, – искомое число из распределения
Практика:
В языке R есть функции для генерации чисел из большинства популярных распределений, например (в дальнейшем, – число генерируемых чисел):
runif(n, min = 0, max = 1)
Равномерное непрерывное распределение
rnorm(n, mean = 0, sd = 1)
Нормальное распределение
rpois(n, lambda)
Распределение Пуассона
?distributions
Чтобы посмотреть все остальные в справке
Однако, часто возникает необходимость в генерации чисел из какого-нибудь специфического распределения, для которого в языке R функция не предусмотрена.
Метод обратной функции:
1) На бумажке (или в каком-нибудь вольфраме символьно) считаем функцию . Для этого надо выразить через из уравнения для
2) Печатаем полученную функцию в скрипте языка R:
G <- function(y)
{
…
}
3) В консоли или отдельном скрипте пишем такую штуку и получаем выборку нашего распределения из значений:
x <- apply(runif(n), G)
Lambda1 lambda2
2.149147 3.528324
Постановка задачи
Пусть случайная величина Х имеет функцию распределения .
Односторонние интервалы: – левый; – правый.
Двусторонний интервал: .
Доверительная вероятность: .
Задача: по заданной доверительной вероятности требуется вычислить .
Решение этой задачи находится через квантили распределения.
Квантиль распределения:
– неявное уравнение. Требуется найти – квантиль уровня .
Пример
, -фикс.
gm <- 0.9
f <- function(a) {
qnorm(gm + a, 0, 1) – qnorm(a, 0, 1)
}
curve(f, 0, 1 - gm)
Пусть
2*qnorm(0.95,0,1) # = 3.289707 – минимум функции на графике выше
Гипотезы
– нулевая гипотеза
– альтернативная гипотеза
Простая гипотеза – гипотеза, которой удовлетворяет только одно распределение вероятности
Сложная гипотеза – несколько распределений
Data: x.t
X-squared = 6.913, df = 10, p-value = 0.7336
chisq.test(x.t, p=pr.wrong) #отвергаем
Data: x.t
X-squared = 55.28, df = 8, p-value = 3.895e-09
15. Методы проверки нормальности данных при неизвестных параметрах распределения (вероятностная бумага, критерии Лильи – Фора).
Теория:
Практика в R:
library("nortest")
lillie.test(iris$Sepal.Width)
Z test for the mean value
data:
= -0.59725, p-value = 0.7248
Decision Accepted
Power 0.9996034
Welch Two Sample t-test
data: log(myx$MYXA30) and log(myx$MYXA60)
t = -0.37751, df = 29.213, p-value = 0.7085
Mean of x mean of y
2.807217 2.960523
Data: x by grp
Kruskal-Wallis chi-squared = 1.712, df = 2, p-value = 0.4249
Pearson's Chi-squared test
Data: tab
X-squared = 568.57, df = 12, p-value < 2.2e-16
Rho
0.9547151
Data: tab
p-value = 0.009423
0.001034782 0.656954980
sample estimates:
Odds ratio
0.05851868
Модель Гаусса-Маркова
,
(1 вход и 1 выход)
– неслучайная величина
– н.о.р.
Неизвестные параметры
– наблюдаемые погрешности
0.7143
plot(alcohol$x, alcohol$y)
abline(alcohol.lm, col="blue")
summary(alcohol.lm)
Call:
lm(formula = y ~ x - 1, data = alcohol)
Residuals:
1 2 3 4 5 6
E-14 -6.228e-16 2.552e-15 -4.459e-15 -4.270e-15 -4.081e-15
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 7.143e-01 5.524e-17 1.293e+16 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
587.320391 -3.942494 0.008115
plot (gal$h, gal$d)
lines (gal$h, gal.lm$fitted.values)
summary(gal.lm)
Call:
lm(formula = d ~ h + I(h^2), data = gal)
Residuals:
1 2 3 4 5 6 7
-9.311 19.672 3.800 -9.883 -24.198 17.992 1.928
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.873e+02 1.220e+02 4.813 0.008564 **
h -3.942e+00 6.134e-01 -6.428 0.003012 **
I(h^2) 8.115e-03 7.351e-04 11.040 0.000383 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Shapiro-Wilk normality test
data: gal.lm$residuals
W = 0.94492, p-value = 0.6834
(gal.lm3 = update(gal.lm, . ~ . + I(h ^ 3)))
Call:
lm(formula = d ~ h + I(h^2) + I(h^3), data = gal)
Coefficients:
(Intercept) h I(h^2) I(h^3)
-6.902e+01 1.260e+00 -4.983e-03 1.056e-05
plot (gal$h, gal$d)
lines (gal$h, gal.lm3$fitted.values)
summary(gal.lm3)
Call:
lm(formula = d ~ h + I(h^2) + I(h^3), data = gal)
Residuals:
1 2 3 4 5 6 7
-1.808 6.237 -1.846 -4.108 -14.097 24.453 -8.832
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.902e+01 5.020e+02 -0.137 0.899
h 1.260e+00 3.921e+00 0.321 0.769
I(h^2) -4.983e-03 9.792e-03 -0.509 0.646
I(h^3) 1.055e-05 7.872e-06 1.341 0.272
Shapiro-Wilk normality test
data: gal.lm3$residuals
W = 0.89196, p-value = 0.285
(gal.lm4 = update(gal.lm3, . ~ . - I(h ^ 3)))
Call:
lm(formula = d ~ h + I(h^2), data = gal)
Coefficients:
(Intercept) h I(h^2)
587.320391 -3.942494 0.008115
summary(gal.lm4)
Call:
lm(formula = d ~ h + I(h^2), data = gal)
Residuals:
1 2 3 4 5 6 7
-9.311 19.672 3.800 -9.883 -24.198 17.992 1.928
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.873e+02 1.220e+02 4.813 0.008564 **
h -3.942e+00 6.134e-01 -6.428 0.003012 **
I(h^2) 8.115e-03 7.351e-04 11.040 0.000383 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Min 1Q Median 3Q Max
-4.6195 -1.1002 -0.1656 1.7451 4.1976
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.977e-01 9.636e-01 -0.309 0.76
I(Height * Girth^2) 2.124e-03 5.949e-05 35.711 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
1 2
24.997246 1.576019
Выборка, выборочное распределение одномерной случайной величины. Построение эмпирической плотности и эмпирической функции распределения в пакете R.
Теория:
Эмпирическую плотность можно оценить с помощью гистограммы, ядерных оценок.
K(x) – четная, интегрируемая:
Пусть h > 0, тогда:
(x) =
где , i = 1:n – значения выборки
Рассмотрим некоторые ядра:
Практика в R:
density(x, bw = "nrd0", adjust = 1,
kernel = c("gaussian", "epanechnikov", "rectangular",
"triangular", "biweight",
"cosine", "optcosine"),
weights = NULL, window = kernel, width,
give.Rkern = FALSE,
n = 512, from, to, cut = 3, na.rm = FALSE, ...)
Обязательные параметры: выборка, ядро, используя параметр kernel. Так же можно указать критерий по параметру сглаживания или его численное значение - параметр bw = "nrd0", (другие. так же указывать как строковую константу. nrd, ucv, bcv, SJ).
ecdf(x)
На вход только выборку.
Также добавив приставку d – плотность или p – функция распределения, можно получить соответственно плотность или функцию распределения указанного после приставки закона распределения:
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
dunif(x, min = 0, max = 1, log = FALSE)