# Пререквизиты ### Математический анализ 1. Производная функции 2. Нахождение экстремума функции 3. Предел 4. Обозначение суммы ### Теория вероятности 1. Закон распределения дискретной и непрерывной случайной величины 2. Функция плотности и функция вероятности непрерывной случайной величины 3. Теоретический момент распределения. Определение, виды 4. Математическое ожидание и дисперсия. 5. Основные распределения ### Метаматическая статистика 1. Определения выборки и генеральной совокупности 2. Статистический критерий 3. Тестирование гипотез # Расширенный план урока *Если мы предполагаем, что студент совсем не знаком с математической статистикой, то я начну урок с ее основ. Дам базовые определения и обозначения. Подкрепляя определения примерами из жизни.* ## Выборка и Генеральная совокупность. Статистика - это искусство делать выводы о многих объектах, имея в распоряжении лишь их малую часть. **Генеральная совокупность** - это все объекты, которые мы хотим изучить. Число объектов генеральной совокупности обычно обозначается заглавной буквой $N$. **Выборка** - это все объекты, доступные нам для изучения. Выборка является частью генеральной совокупности. Число объектов выборки обозначается маленькой буквой $n$. > Предположим, мы исследуем цены на рынке подержанных автомобилей. Для изучения мы выбрали 100 первых попавшихся объявлений о продажах подержанных автомобилей. Все подержанные автомобили - это наша ***генеральная совокупность***. Выбранные 100 автомобилей из объявлений - это наша ***выборка***. > Размер нашей выборки $n = 100$, а вот размер генеральной совокупности $N$ нам не известен. **Наблюдение** - это имеющаяся у нас информация об одном из объектов выборки. Наблюдение обозначается $x_i$, где $i$ - это порядковый номер наблюдения, а $x_i$ - само значение нашего наблюдения. > Допустим, в нашем примере с автомобилями, цена автомобиля из первого объявления равна 1 млн. рублей, а цена из последнего объявления оказалась 650 тыс. рублей. Математически мы обозначим это $x_1 = 1000000$, $x_{100} = 650000$. Далее я перехожу к определению точечных оценок на примере оценок математического ожидания и дисперсии. Раскрываю тему смещенности и несмещенности оценок ## Выборочные аналоги (точечные оценки) математического ожидания и дисперсии Тут я рассказываю о том, что среднее является несмещенной и состоятельной оценкой математического ожидания. Ввожу определение несмещенности оценки, состоятельности и эффективности. Понятие асимптотической несмещенности Рассказываю про эти три оценки и прохожусь по их свойствам (выполняются или нет): * несмещенность * состоятельность * эффективность | Характеристика | Теоретический момент | Точечная оценка | Свойства | | ----------------------- | ------------------------------------------- | -------------------------------------------------------------- | ------------------ | | Математическое ожидание | $M[X] = \sum_{i=1}^n x_ip_i$ | $$\bar x = \frac{1}{n}\sum_{i=1}^n x_i $$ | Несмещенная оценка | | Дисперсия | $$Var[X] = \sum_{i=1}^n (x_i- M[X])^2p_i$$ | $$\hat\sigma^2_n = \frac{1}{n}\sum_{i=1}^n (x_i - \bar x)^2 $$ | Смещенная оценка, асимптотически несмещенная | | Дисперсия | $$Var[X] = \sum_{i=1}^n (x_i- M[X])^2p_i$$ | $$\hat s^2_n = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar x)^2 $$ | Несмещенная оценка | Рассказать кратко про закон больших чисел и как с ростом выборки среднее стремится к математическому ожиданию ### Упражнение 1 Следующий код в python генерирует 5 случайных чисел из нормального распределения $N(5, 1)$ ``` import numpy as np np.random.seed(2222) mu, sigma = 5, 1 # мат ожидание и стандартное отклонение sample = np.random.normal(mu, sigma, 5) sample = np.around(sample, 2) print(sample) # [4.06 6.01 5.46 5.72 5.55] ``` Вывод: ![](https://hackmd.io/_uploads/ByQAarDOn.png) Найти эмпирическую оценку математического ожидания, смещенную и несмещенную оценки дисперсии и сравнить их с теоретическими моментами. В поле с ответом ввести 3 ответа, округлив до 2 знаков после запятой Ответы: 5.36; 0.46; 0.57 Решение: Эмпирическая оценка мат ожидания: $\hat \mu = \bar x = \frac{1}{n}\sum_{i=1}^n x_i = \frac{1}{5}\sum_{i=1}^5 x_i = \frac{x_1 + x_2 + x_3 + x_4 + x_5}{5} = \frac{4,06 + 6,01 + 5,46 + 5,72 + 5,55}{5} = 5,36$ Эмпирическая смещенная оценка дисперсии: $\hat\sigma^2_n = \frac{1}{n}\sum_{i=1}^n (x_i - \bar x)^2 = \frac{1}{5}\sum_{i=1}^5 (x_i - 5,36)^2 =$ $= \frac{(4,06-5,36)^2 + (6,01-5,36)^2 + (5,46-5,36)^2 + (5,72-5,36)^2 + (5,55-5,36)^2}{5} = 0.46$ Эмпирическая несмещенная оценка дисперсии: $\hat s^2_n = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar x)^2 = \frac{1}{4}\sum_{i=1}^5 (x_i - 5,36)^2 =$ $= \frac{(4,06-5,36)^2 + (6,01-5,36)^2 + (5,46-5,36)^2 + (5,72-5,36)^2 + (5,55-5,36)^2}{4} = 0.57$ > Я долго думала, включать Метод Моментов или нет, и решила, что не стоит перегружать информацией. ММ реже используется и часто дает смещенные оценки, поэтому сразу к ММП ## Оценки параметров распределений методом Максимального Правдоподобия Мы уже умеем находить точечные оценки для математического ожидания и дисперсии. В этой секции мы научимся находить различные параметры распределений при помощи метода максимального правдоподобия ММП (Maximum Likelihood Estimator - MLE) Определение функции правдоподобия - как она составляется и из чего состоит Описание алгоритма подбора параметров методом МП Вывод параметров для нормального распределения Свойства ММП оценок: * несмещенная или асимптотически несмещенная оценка * состоятельная оценка * эффективная или асимптотически эффективная оценка ### Упражнение 2 Перед вами выборка размера 10 из распределения Пуассона. Найти оценку параметра $\hat \lambda$ этого распределения Методом Максимального Правдоподобия. | x1 | x 2 | x3 | x4 | x5 | x6 | x7 | x8 | x9 | x10| | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | | 4 | 7 | 6 | 11 | 10 | 7 | 5 | 3 | 12 | 4 | Ответ: $\hat \lambda = \bar x = 6.9$ Код, сгенеривший это: ``` import numpy as np np.random.seed(1111) lam = 7 sample = np.random.poisson(lam, 10) print(sample) # [ 4 7 6 11 10 7 5 3 12 4] ``` Решение: Построим функцию правдоподобия $$L(\theta) = \prod_{i=1}^np(x_i, \theta) = \prod_{i=1}^{10}\exp(-\theta)\frac{1}{x_i!}\theta^{x_i} $$ Возьмем логарифм $$ l(\theta) = \ln \prod_{i=1}^{10}\exp(-\theta)\frac{1}{x_i!}\theta^{x_i}$$ Упростим выражение $$ l(\theta) = \sum_{i=1}^{10}\ln[\exp(-\theta)\frac{1}{x_i!}\theta^{x_i}] = -\sum_{i=1}^{10}\theta + \sum_{i=1}^{10} \ln \frac{1}{x_i!} + \sum_{i=1}^{10} \ln \theta^{x_i}$$ $$ l(\theta) = -10*\theta - \sum_{i=1}^{10} \ln x_i! + \ln \theta \sum_{i=1}^{10} x_i $$ Возьмем производную по $\theta$ $$ l'(\theta) = -10 + \frac{1}{\theta} \sum_{i=1}^{10} x_i$$ Возьмем вторую производную по $\theta$ $$ l''(\theta) = -\frac{1}{\theta^2} \sum_{i=1}^{10} x_i$$ Вторая производная < 0, значит мы находим максимум Приравняем первую производную к нулю $$ l'(\theta) = 0 $$ Отсюда: $$-10 + \frac{1}{\theta} \sum_{i=1}^{10} x_i = 0$$ $$\frac{1}{\theta} \sum_{i=1}^{10} x_i = 10$$ $$\theta = \frac{1}{10} \sum_{i=1}^{10} x_i = \bar x$$ $$\hat \theta = \bar x = \frac{4+7+6+11+10+7+5+3+12+4}{10} = 6,9$$ Так как распределение Пуассона - дискретное распределение, мы смело можем округлить нашу оценку и использовать $\hat \lambda = 7$ ## Как подобрать распределение Теперь, когда мы знаем, как под каждое конкретное распределение подбирать параметры, мы фокусируемся на том, а как же выбрать нужное распределение? ### Возможный алгоритм 1) Визуально выбрать по гистограмме "кандидатов" 3) Рассчитать MLE оценки 4) Проверить критерием Пирсона или Колмогорова-Смирнова/Андерсона (последние два только для непрерывных распределений) Тут я более подробно описываю алгоритм по шагам. Даю теорию по основным критериям Хи-квадрат Пирсона Критерий Колмогорова-Смирнова (только для непрерывных распределений) Андерсона ### Пример 1 В этом примере я делаю сэмпл из нормального распределения ``` import numpy as np import matplotlib.pyplot as plt np.random.seed(1111) sample = np.random.normal(10, 1, 1000) plt.hist(sample, bins = 10) plt.show() ``` ![](https://hackmd.io/_uploads/SknMvqsu2.png) Вычисляем ММП оценки ``` mu, sigma = scipy.stats.norm.fit(sample) print(mu, sigma) #10.03065504293586 0.9735639311032354 ``` Проводим тест на соответсвие распределению с данными параметрами ``` # Проверяем с помощью теста Колмогорова-Смирнова res = stats.goodness_of_fit(dist = scipy.stats.norm, known_params = {"loc": mu, "scale": sigma}, data = sample, statistic = "ks" ) print(res) ``` Результат: ``` GoodnessOfFitResult(fit_result= params: FitParams(loc=10.03065504293586, scale=0.9735639311032354) success: True message: 'The fit was performed successfully.', statistic=0.019875455307503898, pvalue=0.8148, null_distribution=array([0.03361172, 0.01984966, 0.04192832, ..., 0.02066031, 0.02106213, 0.02675832])) ``` ### Упражнение 3 Мы провели тест Колмогорова-Смирнова для наших данных. Мы предположили что наши данные Имеют Нормальное распределение. Мы приняли уровень значимости $\alpha = 0.05$ и провели тест Колмогорова-Смирнова с помощью библиотеки python scipy, которая в качестве вывода выдала нам следующий результат результат: ``` statistic=0.019875455307503898, pvalue=0.8148 ``` #### Какой вывод мы можем сделать? Выберите правильный ответ. 1. Наши данные описываются Нормальным распределением с предполагаемыми параметрами. 2. Наши данные НЕ описываются Нормальным распределением предполагаемыми параметрами. 3. Этот тест не применим в нашем случае. #### Ответ: (2) Наши данные описываются Нормальным распределением с предполагаемыми параметрами. p-value больше чем уровень значимости (0.8148>0.05), а это значит что мы не отвергаем нулевую гипотезу о том, что данные распределены согласно Нормальному закону с установленными нами ранее ММП параметрами 10.03 и 0.97. ### Пример 2 Показываю датасет с реальными данными о времени заполнения адреса юзерами одного сервиса доставки еды. Предлагаю предположить какие распределения могут подойти под эти данные. ![](https://hackmd.io/_uploads/H1J8DdsO3.png) Мое предположение по визуальному анализу: гамма-распределение. Первым делом я с помощью библиотеки вычисляю ММП параметры гамма-распределения ``` from scipy import stats # Выделяем нашу выборку data = d["user_time"] # Вычисляем ММП оценки параметров a, loc, scale = scipy.stats.gamma.fit(data) # Выводим ММП оценки параметров print(a, loc, scale, sep = '\n') ``` ``` # Вывод: 2.2409926143142442 1.9836253584605417 10.22655447976155 ``` После того как у нас готовы ММП оценки мы можем провести тест на то, насколько хорошо наши данные Гамма описываются распределением с данными параметрами ``` # Проверяем с помощью теста Колмогорова-Смирнова res = stats.goodness_of_fit(dist = scipy.stats.gamma, known_params = {"a": a, "loc": loc, "scale": scale}, data = data, n_mc_samples = 1000, statistic = "ks" ) print(res) ``` ``` # Вывод: GoodnessOfFitResult(fit_result= params: FitParams(a=2.2409926143142442, loc=1.9836253584605417, scale=10.22655447976155) success: True message: 'The fit was performed successfully.', statistic=0.08745078149029262, pvalue=0.000999000999000999 ``` Построим результат ``` import matplotlib.pyplot as plt x = np.arange(d.user_time.min(), d.user_time.max()) pdf = scipy.stats.gamma.pdf(x, a=a, scale=scale, loc=loc) plt.figure plt.hist(d.user_time, density = True) plt.plot(x, pdf) plt.show() ``` ![](https://hackmd.io/_uploads/Bk5OC_sdh.png) ### Упражнение 4 Мы провели тест Колмогорова-Смирнова для наших данных. Мы предположили что наши данные Имеют Гамма распределение. Мы приняли уровень значимости $\alpha = 0.01$ и провели тест Колмогорова-Смирнова с помощью библиотеки python scipy, которая в качестве вывода выдала нам следующий результат результат: ``` statistic=0.08745078149029262, pvalue=0.000999000999000999 ``` #### Какой вывод мы можем сделать? Выберите правильный ответ. 1. Наши данные описываются Гамма распределением с предполагаемыми параметрами. 2. Наши данные не описываются Гамма распределением предполагаемыми параметрами. 3. Этот тест не применим в нашем случае. #### Ответ: (2) Наши данные не описываются Гамма распределением с предполагаемыми параметрами. p-value меньше чем уровень значимости (0.0009<0.01), а это значит что мы отвергаем нулевую гипотезу о том, что данные распределены согласно Гамма закону с установленными нами ранее ММП параметрами. > После этого я прогнала всевозможные распределения из scipy но ни ондно не зафитилось до сотояния соответсвия критерию. Об этом я уже не стала писать, звучит как задел на практический урок. Почему так произошло - работа с реальными данным/ несовершенство реальных данных Как тогда выбирать - рассказать про SSE и выбор среди всех распределений по этому критерию # Конспект с важной теорией **Генеральная совокупность** - это все объекты, которые мы хотим изучить. Число объектов генеральной совокупности обычно обозначается заглавной буквой $N$. **Выборка** - это все объекты, доступные нам для изучения. Выборка является частью генеральной совокупности. Число объектов выборки обозначается маленькой буквой $n$. Пусть у нас есть $n$ наблюдений $x_1, x_2,...,x_n$ из генеральной совокупности. **Среднее значение** - это эмпирическая оценка математического ожидания. Вычисляется по формуле $$ \bar x = \frac{1}{n}\sum_{i=1}^n x_i $$ Среднее значение является несмещенной оценкой математического ожидания. $$ \hat\mu = \bar x$$ | Характеристика | Теоретический момент (дискретные распределения) | Точечная оценка | Свойства | | ----------------------- | ------------------------------------------- | -------------------------------------------------------------- | ------------------ | | Математическое ожидание | $M[X] = \sum_{i=1}^n x_ip_i$ | $$\bar x = \frac{1}{n}\sum_{i=1}^n x_i $$ | Несмещенная оценка | | Дисперсия | $$Var[X] = \sum_{i=1}^n (x_i- M[X])^2p_i$$ | $$\hat\sigma^2_n = \frac{1}{n}\sum_{i=1}^n (x_i - \bar x)^2 $$ | Смещенная оценка | | Дисперсия | $$Var[X] = \sum_{i=1}^n (x_i- M[X])^2p_i$$ | $$\hat s^2_n = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar x)^2 $$ | Несмещенная оценка | ## ## Оценки параметров распределений методом максимального правдоподобия #### Функция правдоподобия $$L(\theta) = p(x_1, \theta) * p(x_2, \theta) * ... * p(x_n, \theta) = \prod_{i=1}^np(x_i, \theta)$$ Где $L(\theta)$ - функция правдоподобия, $p(x_i, \theta)$ - функция плотности распределения c пока еще неизвестным параметром $\theta$ в точке $x_i$. Для упрощения процесса вычисления производной функции правдоподобия обычно вычисляется производная логарифма данной функции. Экстремум будет тот же самый. ### Алгоритм вычисления оценки Методом Максимального Правдоподобия 1. Построить функцию правдоподобия $L(\theta)$ 2. Взять логарифм этой функции $lnL(\theta)$ 3. Взять производную $lnL(\theta)'_{\theta}$ 4. Взять вторую производную $lnL(\theta)''_{\theta}$ и проверить достаточное условие максимума 5. Приравнять к нулю производную $lnL(\theta)'_{\theta} = 0$ 7. Найти точечную оценку $\hat \theta = f(x_1,...x_n)$ как функцию от $x_1,...,x_2$ ## MLE оценки для некоторых распределений ### Нормальное распределение $\xi \sim N(\mu, \sigma^2)$ Параметр | MLE-оценка | ---------- | ------ | $\mu$ | $$\hat \mu = \frac{1}{n}\sum_{i=1}^n x_i$$ | | $\sigma^2$ | $$\hat\sigma^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar x)^2 $$| ## Критерии для подбора распределения ### Критерий Хи-квадрат (Chi-square goodness of fit test) - универсальный: для дискретных и непрерывных распределений - применяется к данным, сгруппированным по интервалам (bins) $H_0$: наши данные подчиняются заданному распределению $H_1$: наши данные не подчиняются заданному распределению При условии, что наши данные разбиты на k интервалов (bins), статистика вычисляется #### Статистика: $$ \chi^2 = \sum_{i=1}^k(O_i - E_i)^2/E_i $$ Где $O_i$ эмпирическая частота для интервала i, а $E_i$ is теоретическая (ожидаемая) частота для интервала i. Теоретическая частота вычисляется как: $$E_i=n*(F(Y_u)−F(Y_l))$$ Где $Y_u$ - верхний предел для интервала i (upper), $Y_l$ - нижний предел для интервала i (lower), n - размер выборки (sample size). #### Вывод: Тестовая статистика распределена по $\chi_{k-c}^2$ Хи-квадрат с (k - c) степенями свободы. Где k - число интервалов, c = (число параметров тестируемого распределения + 1). #### Интерпретация: Если p-value больше, чем уровень значимости $\alpha$ - значит принимаем $H_0$, что наши распределены так, как мы предполагали Если p-value меньше, чем уровень значимости $\alpha$ - значит отвергаем $H_0$, принимаем $H_1$, что наши не распределены так, как мы предполагали ### Критерий Колмогорова-Смирнова #### Непараметрический критерий $H_0$: наши данные подчиняются заданному распределению $H_1$: наши данные не подчиняются заданному распределению #### Статистика: $$ D = \max_{i=1..N} ( F(Y_i) - \frac{i-1}{N}, \frac{i}{N} - F(Y_i)) $$ Где $F(Y)$ - это функция предполагаемого нами распределения #### Интерпретация: Если p-value больше, чем уровень значимости $\alpha$ - значит принимаем $H_0$, что наши распределены так, как мы предполагали Если p-value меньше, чем уровень значимости $\alpha$ - значит отвергаем $H_0$, принимаем $H_1$, что наши не распределены так, как мы предполагали #### Ограничения: * только для непрерывных распределений #### SSE Это сумма квадратов ошибок $$ SSE = \sum_{i=1}^ne_i^2$$ Где $e_i$ - это разница между наблюдаемым исходом и модельным/предсказанным исходом. Один из возможножных способов оценки, насколько хорошо распределение описывает данные. Мы можем использовать: $$e_i = F(Y_i) - \hat F (Y_i)$$ где $F(Y_i)$ - теоретическая функция распределения, $\hat F (Y_i)$ - эмпирическая функция распределения Чем меньше SSE, тем лучше распределение описывает наши эмпирические данные.