R Notebook

Методы статистической обработки информации

---

Практика 1. Моделирование и описательная статистика

Выполнил Глушков Егор Александрович, гр. 20М.04-мм

Варианты №4, 7

№4. Отрицательно биномиальное распределение

\(\beta_{neg}(k=3, p=0.4);\;P\{\xi=j\}=\frac{\Gamma(k+j)}{\Gamma(k)j!}p^k(1-p)^j,\;j=0,1,...,\infty\)

Моделируем данные по заданным параметрам \(k=3,\;p=0.4\)

library(moments)

N <- 1500
sz <- 3; p <- 0.4
X <- rnbinom(n=N, size=sz, prob=p)

На модельных данных строим гистограмму

hist(X)

Описательная статистика

Находим матожидание, дисперсию, стандартное отклонение, ошибку среднего и сравниваем с гипотетическими параметрами

expectation <- c(mean(X), sz * (1 - p) / p)
variance <- c(sd(X)^2, sz * (1 - p) / p^2)
st_dev <- c(sd(X), sqrt(variance[2]))
err <- c(st_dev[1] / sqrt(N), st_dev[2] / sqrt(N))

df1 <- data.frame(expectation, variance, st_dev, err, row.names = c("simulated", "hypothetical")); df1

Находим медиану, минимум, максимум, размах

c(median=median(X), min=min(X), max=max(X), R=max(X)-min(X))

median    min    max      R 
     4      0     19     19 

Находим квартили и интерквартильный размах

q <- quantile(X, probs=seq(0, 1, 0.25)); q

  0%  25%  50%  75% 100% 
   0    2    4    6   19 

Q <- q[4] - q[2]; unname(Q)

[1] 4

c(q25<-qnbinom(0.25, size=sz, prob=p), q75<-qnbinom(0.75, size=sz, prob=p), intqR=q75-q25)

            intqR 
    2     6     4 

Центральные моменты 2-4 порядков по формуле и через библиотеку [совпадают]

m2 <- moment(X, order=2, central=TRUE)
m3 <- moment(X, order=3, central=TRUE)
m4 <- moment(X, order=4, central=TRUE)
c(mean((X-mean(X))^3), m3)

[1] 37.27772 37.27772

c(mean((X-mean(X))^4), m4)

[1] 507.9319 507.9319

Подсчет коэффициентов асимметрии и эксцесса по формулам и через библиотеку [совпадают]

print(paste("skewness: ", m3 / m2 ^ 1.5, " (formula), ", skewness(X), " (lib)"))

[1] "skewness:  1.04018762905602  (formula),  1.04018762905602  (lib)"

print(paste("kurtosis: ", m4 / m2 ^ 2 - 3, " (formula), ", kurtosis(X) - 3, " (lib; corrected)"))

[1] "kurtosis:  1.2988789558756  (formula),  1.2988789558756  (lib; corrected)"

gamma1 <- (2 - p) / sqrt(sz * (1 - p))
gamma2 <- 6 / sz + p ^ 2 / (sz * (1 - p))
print(paste("excess sk and kur: ", gamma1, gamma2, " (by formulae)"))

[1] "excess sk and kur:  1.19256958799989 2.08888888888889  (by formulae)"

Характеристики асимметрии и эксцесса в целом соответствуют форме заданного распределения.

Оценка параметров

Метод моментов

Приравниваем выборочные моменты к соответствующим моментам распределения (при 2 неизвестных параметрах, т.е. используем 2 момента): \(\alpha_1 = \frac{k(1-p)}{p}\) приравниваем к \(\overline{x}\), \(\mu_2 = \frac{k(1-p)}{p^2} = \frac{\overline{x}}{p}\) приравниваем к \(m_2\). Получаем \(p_{mm} = \frac{\overline{x}}{m_2}\), \(k_{mm} = \frac{\overline{x}^2}{m_2 - \overline{x}}\)

p. <- mean(X) / m2
k. <- mean(X)^2 / (m2 - mean(X))
c(p.=p., k.=k.)

       p.        k. 
0.4104918 3.1070206 

Моделируем исследуемое распределение с теоретическими и с оценочными (по методу моментов) параметрами. Строим график

xx. <- rnbinom(1500, size=k., prob=p.)

hist(X, freq=FALSE, xlim=c(-1, 25), main="Method of moments")
lines(density(X), col=2)

lines(density(xx.), col=3)
legend('topright',c("hyp","mm"),pch=20,col=c(2,3))

Метод максимального правдоподобия

Оценки параметров по методу максимального правдоподобия

Func.prob.log <- function(x) - sum(dnbinom(X, size = x[1], prob = x[2], log = TRUE))
res <- optim(c(k., p.), Func.prob.log) # k. and p. as initial values
k.. <- res$par[1]; p.. <- res$par[2];
c(k..=k.., p..=p..)

      k..       p.. 
2.9438163 0.3974773 

Построим график со сравнением теоретического распределения с распределениями, построенными по параметрам, полученным в методах моментов (mm) и максимального правдоподобия (mle)

xx.. <- rnbinom(1500, size=k.., prob=p..)

hist(X, freq=FALSE, xlim=c(-1, 25), main="MM, MLE")
lines(density(X), col=2)

lines(density(xx.), col=3)
lines(density(xx..), col=4)

legend('topright',c("hyp","mm", "mle"),pch=20,col=c(2,3,4))

В случае, когда параметр k (т.е. size) известен и требуется найти p, то оценка параметра p по методу максимального правдоподобия, будет иметь вид: \(p_{mle} = \frac {k} {k + \overline{x}}\)

k <- sz
p_mle <- k/(k+mean(X)); p_mle

[1] 0.402037

Заметим, что оценки, полученные как по методу моментов, так и по методу максимального правдоподобия, достаточно близки к теоретическим. Для получения более точных оценок следует увеличить объем моделируемых данных

Проверка согласия эмпирического и теоретического распределения по критерию хм-квадрат Пирсона

# эмпирические частоты
h <- hist(X, plot=FALSE)

# гипотетические частоты
p.i <- sapply(seq(length(h$breaks)-1) + 1, function(i) pnbinom(h$breaks[i], size=k, prob=p) - pnbinom(h$breaks[i-1], size=k, prob=p))
p.i[1] <- pnbinom(h$breaks[2], size=k, prob=p)
p.i[length(p.i)] <- 1 - pnbinom(h$breaks[length(h$breaks) - 1], size=k, prob=p)
sum(p.i) # должно быть равно 1

[1] 1

Проверка условия \(np_i \geq 5\) (для применения критерия согласия Пирсона)

tab <- cbind(h$counts, length(X) * p.i); tab

      [,1]       [,2]
 [1,]  487 476.160000
 [2,]  374 393.984000
 [3,]  295 282.175488
 [4,]  164 169.305293
 [5,]  100  91.521604
 [6,]   45  46.182614
 [7,]   20  22.192731
 [8,]    7  10.282381
 [9,]    6   4.631134
[10,]    2   3.564754

Объединяем строки с малыми значениями

Tab <- rbind(tab[1:8, ], c(sum(tab[9:nrow(tab), 1]), sum(tab[9:nrow(tab), 2]))); Tab

      [,1]       [,2]
 [1,]  487 476.160000
 [2,]  374 393.984000
 [3,]  295 282.175488
 [4,]  164 169.305293
 [5,]  100  91.521604
 [6,]   45  46.182614
 [7,]   20  22.192731
 [8,]    7  10.282381
 [9,]    8   8.195889

Наблюдаемое значение статистики \(\chi^2\)

chi2 <- sum(apply(Tab, 1, function(x) (x[1] - x[2])^2 / x[2])); chi2

[1] 4.09438

alpha <- 0.05
p_value <- 1 - pchisq(chi2, nrow(Tab) - 2 - 1); 
print(paste("p_value", round(p_value, 3), sep=" = "))

[1] "p_value = 0.664"

if (p_value > alpha) print(paste("The hypothesis of the agreement of the empirical distribution with the theoretical one is not rejected with the significance level alpha", alpha, sep="=")) else print("The hypothesis of the agreement of the empirical distribution with the theoretical one is rejected")

[1] "The hypothesis of the agreement of the empirical distribution with the theoretical one is not rejected with the significance level alpha=0.05"

\(p.value = 0.664\), поэтому гипотеза о согласии эмпирического и теоретического распределения НЕ отвергается при заданном \(\alpha = 0.05\).