ЛР1_Вар26_ШабуняАлексей_МС4.Rmd

---
title: "ЛР1_Вар26_ШабуняАлексей_БСТ214"
author: "Шабуня Алексей Андреевич"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
  pdf_document:
    latex_engine: xelatex
    df_print: paged
  word_document: default
lang: "ru-russian"
---

## 1.1. Импорт данных для работы

Подключаем необходимый для импорта данных пакет командой `library`:

```{r}
library('rio')
```

Присваиваем импортированным данным имя `dataframe`:

```{r}
dataframe <- import('Данные по вариантам для выполнения типового расчета.xlsx')
```

Сохраним теперь отдельным объектом данные по исследуемому показателю, которые сохранены в файле в столбце v26:

```{r}
data <- dataframe['v26']
```

Посмотрим на начало таблицы, используя функцию `head`:

```{r}
head(data) # показывает первые шесть (по умолчанию) строк каждого столбца
```

...и на конец таблицы, используя функцию `tail`:

```{r}
tail(data) # показывает последние шесть (по умолчанию) строк каждого столбца
```

Для начального знакомства с данными будет полезен результат вызова функции `str`, которая позволяет узнать класс объекта, в котором содержатся данные (в нашем случае `data.frame`), число наблюдений (`observations` -- 101) и столбцов-признаков (`variable` -- 1), а также тип данных каждого признака (указан после имени -- `chr`) и первые четыре наблюдения:

```{r}
str(data)
```

Данные успешно импортированы. Теперь необходимо привести их к виду, при котором мы сможем использовать функции статистического анализа R.

## 1.2. Подготовка данных для статистического анализа

Можно заметить, что к нашим 100 наблюдениям добавилось одно лишнее, которое было внизу столбца, и тип данных не числовой (`numeric`), а текстовый (`character`). Исправить все можно так:

```{r}
data <- data[1:100,] # отрезает в data первые 100 наблюдений
data <- as.data.frame(as.numeric(data)) # изменяет формат данных на числовой
```

Переименуем столбец данных в соответствии с исследуемой величиной:

```{r}
names(data) <- 'Производственная площадь ста\n предприятий машиностроения региона'
```

Проверим структуру исправленных данных, теперь все верно:

```{r}
str(data)
```

Наконец можем перейти к статистическому анализу.

## 1.3. Дескриптивная статистика

Самое начальное представление об исследуемой величине по имеющейся выборке можно получить, вызвав команду `summary` базового пакета:

```{r}
summary(data)
```

Здесь мы видим минимальное значение в выборке `Min.` равное 1.2, выборочные нижний квартиль `1st Qu.` равный 1.7, медиану `Median` -- 2.2, среднюю `Mean` -- 2.270, верхний квартиль `3rd Qu.` -- 2.625, максимальное значение в выборке `Max` равное 3.8.

Теперь подключаем пакет `DescTools` (Tools for Descriptive Statistics), используемый для более подробной информации о распределении изучаемого признака:

```{r}
library('DescTools')
```

Применяем функцию `Desc` (Describe data):

```{r, fig.cap='*Рис. 1. Гистограмма, ящичковая диаграмма и кумулята исследуемой переменной, выдаваемые функцией `Desc` пакета `DescTools`*'}
Desc(data)
```

В результате мы получаем таблицу с кратким описанием признака. Также мы получаем расширенную дескриптивную статистику для каждого признака в датафрейме.

Также функция выдает три графика (рис. 1): а) гистограмму плотности и подогнанную с помощью непараметрических методов кривую плотности распределения, б) ящичковую диаграмму (`boxplot`) -- по ней также можно судить о распределении признака, схожести его с нормальным распределением, но главное ее применение -- выявление наличия экстремальных значений, в) кумуляту относительных частот.

Для удобства дальнейшего анализа запишем данные в вектор под именем `volume`, а также создадим переменную с объёмом выборки `n`:

```{r}
volume <- data[[1]] 
n <- length(volume)
```

## 1.4. Гистограмма эмпирических частот

Построим сначала гистограмму частот. Для этого воспользуемся функцией `hist`:

```{r, fig.cap='*Гистограмма эмпирических частот производственной площади ста предприятий машиностроения региона, построенная с помощью базовой функции `hist`*'}
hist = hist(volume, freq = TRUE, col = 'red',
            breaks = 'sturges', 
            xlab = 'Нижняя граница интервала',
            ylab = 'Частота', ylim = c(0, 40),
            main = 'Гистограмма эмпирических частот')
```

Для нанесения на гистограмму частот (рис. 2) теоретической кривой нормального распределения воспользуемся функцией `lines`, которая соединяет линией точки, координаты которых переданы в функцию аргументами. Где вектор из абсцисс мы сгенерируем функцией `seq`:

```{r}
x_values <- seq(from = min(volume)-0.1,
                to = max(volume)+0.1,length = 1000)
```

Значения теоретической кривой нормального распределения определим по формуле: $$m_{iT} = P(a_i \le X \le b_i) \cdot n = (F(b_i, \hat{\mu}, \hat{\sigma}) - F(a_i, \hat{\mu}, \hat{\sigma})) \cdot n \approx f(x_i, \hat{\mu}, \hat{\sigma}) \cdot h \cdot n,$$ где $\hat{\mu}, \hat{\sigma}$ -- оценки параметров, $h$ -- ширина интервала гистограммы, $f(\cdot)$ -- функция плотности нормального распределения, $F(\cdot)$ -- функция распределения.

Ширину интервалов определим из полученного нами ранее результата вызова функции `hist`. Объект `hist` представляет из себя список, содержащий всю основную информацию о построенной гистограмме частот, первый элемент списка --- список `breaks`, который содержит границы интервалов. Шаг мы определим как разницу между вторым и первым элементами этого списка (т.е. верхней и нижней границы первого интервала):

```{r}
h = unlist(hist['breaks'])[2]-unlist(hist['breaks'])[1]
```

Вектор из значений функции плотности нормального распределения в сгенерированных точках получим, воспользовавшись функцией `dnorm`, аргументами которой являются вектор из координат точек, значение функции плотности в которых нужно определить, средняя и стандартное отклонение совокупности:

```{r}
y_values <- dnorm(x_values, mean = mean(volume), 
                  sd = sd(volume))*h*n
```

Итак, нанесем теоретическую кривую нормального закона на гистограмму частот (рис. 3):

```{r, fig.cap='*Рис. 3. Гистограмма эмпирических частот производственной площади ста предприятий машиностроения региона и теоретическая кривая нормального распределения, построенные с помощью базовой функции `hist`*'}
hist = hist(volume, freq = TRUE, col = 'red', 
            breaks = 'sturges', 
            xlab = 'Нижняя граница интервала', 
            ylab = 'Частота', ylim = c(0, 40), 
            main = 'Гистограмма частот и теоретическая кривая Гаусса') 

lines(x_values, y_values, col = 'blue')
```

## 1.5. Интервальные оценки параметров нормального распределения

Для построения доверительных интервалов используем функции пакета `DescTools`.

Зададим сперва надежность, с которой будем строить все интервальные оценки параметров:

```{r}
gamma = 0.99 #надежность
```

Доверительный интервал для средней вычисляется функцией `MeanCI` (CI от Confidence Interval), аргументами которой необходимо передать вектор значений наблюдений выборки, по которой строится оценка, уровень надежности (`conf.level`, по умолчанию 0.95), если известно генеральное среднеквадратическое отклонение, то указать его значение, в противном случае, по умолчанию используется его выборочная оценка.

Предположим сначала, что нам известно значение генерального среднеквадратического отклонения $\sigma(X)=0.2$ , тогда доверительный интервал для средней строится следующим образом:

```{r}
MeanCI(volume, sd = 0.2, conf.level = gamma)
```

Первое значение результата вызова функции соответствует точечной оценке средней (`mean`) -- среднему арифметическому $\bar{x} = 11.8682$, второе -- нижней границе доверительного интервала (`lwr.ci` от lower bound of the confidence interval), третье -- верхней границе интервала (`upr.ci` от upper bound of the confidence interval).

Таким образом, при известной генеральной дисперсии интервальная оценка генеральной средней: $$P(2.218283 \le \mu \le 2.321317) = 0.99$$

В случае, если генеральная дисперсия неизвестна, доверительный интервал определяется следующим образом:

```{r}
MeanCI(volume, conf.level = gamma)
```

$$P(2.103465 \le \mu \le 2.436135) = 0.99$$

Доверительный интервал для дисперсии можно найти, воспользовавшись функцией `VarCI`:

```{r}
VarCI(volume, conf.level = gamma)
```

Получаем точечную оценку дисперсии $\hat{S}^2=0.4010929$, нижнюю и верхнюю границу доверительного интервала для дисперсии с надежностью 0.99: $$P(0.2856976 \le \sigma^2 \le 0.5970250) = 0.99$$

Для среднеквадратического отклонения определим границы асимптотического доверительного интервала, написав код самостоятельно: $$P\left(\frac{S \cdot \sqrt{2n}}{\sqrt{2n-3} + t} \le \sigma \le \frac{S \cdot \sqrt{2n}}{\sqrt{2n-3} - t}\right) = \gamma = \Phi(t)$$

Квантиль нормального распределения определим функцией `qnorm`, аргументом которой передадим $1-\alpha/2$, где $\alpha=1-\gamma$.

```{r}
sdCI <- c(sd(volume)*sqrt(2*n)/(sqrt(2*n - 3) + qnorm(1-(1-gamma)/2)), 
          sd(volume)*sqrt(2*n)/(sqrt(2*n - 3) - qnorm(1-(1-gamma)/2))) 
names(sdCI) <- c('sd_min', 'sd_max') 
sdCI
```

Таким образом, получаем асимптотический доверительный интервал для генерального среднего квадратического отклонения: $$P(0.5391737 \le \sigma \le 0.7815539 ) = 0.99$$

## 1.6. Интервальные оценки вероятности (генеральной доли) биномиального распределения

### Задача 8 из ДКР1

##### Органы статистики региона хотят оценить долю населения, которая работает в частном секторе, во всем занятом населении. Предполагая, что число работников частных предприятий может быть описано биномиальным законом:

а) найти с надежностью 0.991 границы, в которых будет лежать доля жителей региона, работающих на частных предприятиях, если из 10 случайно отобранных занятых 5 работают на частных предприятиях;

В этом случае имеем дело с *выборкой небольшого объема*, поэтому используем *интервал Клоппера-Пирсона (Clopper-Pearson interval)* построения интервальной оценки вероятности.

Запишем данные условием задачи значения числа испытаний, числа успехов, надежности, с которой надо определить границы доверительного интервала;

Используем функцию `BinomCI` пакета `DescTools`, указав обязательные аргументы, требуемую надежность и метод *Клоппера-Пирсона*:

```{r}
n1 = 10 # число испытаний
m1 = 5 # число успехов
gamma1 = 0.991 # надежность
BinomCI(m1, n1, conf.level=gamma1, method = 'clopper-pearson')
```

Получаем точечную оценку вероятности (`est`), нижнюю границу (`lwr.ci`) и верхнюю границу (`upr.ci`) доверительного интервала, построенного с надежностью 0.991, для доли жителей региона, работающих в частном секторе: $$P(0.1252871 < p < 0.8747129) = 0.991$$

Стоит отметить, что интервальная оценка, полученная *методом Клоппера-Пирсона*, обладает наибольшей шириной. Сравним с интервалами, полученными другими методами.

В случае малого числа испытаний рекомендуют также использовать *интервалы Вилсона (Wilson) и Джеффриса*.

```{r}
BinomCI(m1, n1, conf.level=gamma1, method = 'wilson')
```

Интервал Вилсона: $$P(0.1815783 < p < 0.8184217) = 0.991$$

```{r}
BinomCI(m1, n1, conf.level=gamma1, method = 'jeffreys')
```

Интервал Джеффриса: $$P(0.1548174 < p < 0.8451826) = 0.991$$

Видно, что полученные такими методами интервалы уже, чем построенный нами ранее *интервал Клоппера-Пирсона*.

б) найти с надежностью 0.991 границы, в которых будет лежать доля жителей региона, работающих на частных предприятиях, если из 600 случайно отобранных занятых 384 работают на частных предприятиях;

Запишем данные условием задачи значения числа испытаний, числа успехов, надежности, с которой надо определить границы доверительного интервала;

Используем функцию `BinomCI` пакета `DescTools`, указав обязательные аргументы, требуемую надежность. В случае *выборки большого объема* для построения интервальной оценки вероятности часто используют нормальную аппроксимацию *(интервал Вальда)*:

```{r}
n2 = 600 # число испытаний
m2 = 384 # число успехов
gamma2 = 0.991 # надежность
BinomCI(m2, n2, conf.level=gamma2, method = 'wald')
```

Получаем точечную оценку вероятности (`est`), нижнюю границу (`lwr.ci`) и верхнюю границу (`upr.ci`) доверительного интервала, построенного с надежностью 0.991: $$P(0.5888144 < p < 0.6911856) = 0.991$$

Также при большом числе испытаний можно рассмотреть *интервал Агрести-Коула*, и в нашем случае он дал интервал меньшей ширины:

```{r}
BinomCI(m2, n2, conf.level=gamma1, method = 'agresti-coull')
```

$$P(0.5874805 < p < 0.6893713) = 0.991$$