Физтех.Статистика

Скачать ipynb

Phystech@DataScience

Домашнее задание 8

Правила, прочитайте внимательно:

• Выполненную работу нужно отправить телеграм-боту @miptstats_pds_bot. Для начала работы с ботом каждый раз отправляйте /start. Работы, присланные иным способом, не принимаются.
• Дедлайн см. в боте. После дедлайна работы не принимаются кроме случаев наличия уважительной причины.
• Прислать нужно ноутбук в формате ipynb.
• Выполнять задание необходимо полностью самостоятельно. При обнаружении списывания все участники списывания будут сдавать устный зачет.
• Решения, размещенные на каких-либо интернет-ресурсах, не принимаются. Кроме того, публикация решения в открытом доступе может быть приравнена к предоставлении возможности списать.
• Для выполнения задания используйте этот ноутбук в качестве основы, ничего не удаляя из него. Можно добавлять необходимое количество ячеек.
• Комментарии к решению пишите в markdown-ячейках.
• Выполнение задания (ход решения, выводы и пр.) должно быть осуществлено на русском языке.
• Если код будет не понятен проверяющему, оценка может быть снижена.
• Никакой код из данного задания при проверке запускаться не будет. Если код студента не выполнен, недописан и т.д., то он не оценивается.
• Код из рассказанных на занятиях ноутбуков можно использовать без ограничений.

Правила оформления теоретических задач:

• Решения необходимо прислать одним из следующих способов:
  • фотографией в правильной ориентации, где все четко видно, а почерк разборчив,
    • отправив ее как файл боту вместе с ноутбуком или
    • вставив ее в ноутбук посредством Edit -> Insert Image (фото, вставленные ссылкой, не принимаются);
  • в виде $LaTeX$ в markdown-ячейках.
• Решения не проверяются, если какое-то требование не выполнено. Особенно внимательно все проверьте в случае выбора второго пункта (вставки фото в ноутбук). Неправильно вставленные фотографии могут не передаться при отправке. Для проверки попробуйте переместить ipynb в другую папку и открыть его там.
• В решениях поясняйте, чем вы пользуетесь, хотя бы кратко. Например, если пользуетесь независимостью, то достаточно подписи вида "X и Y незав."
• Решение, в котором есть только ответ, и отсутствуют вычисления, оценивается в 0 баллов.

Баллы за задание:

Легкая часть (достаточно на "хор"):

• Задача 1 — 20 баллов
• Задача 2 — 30 баллов

Сложная часть (необходимо на "отл"):

• Задача 3 — 30 баллов
• Задача 4 — 35 баллов

Легкая часть (достаточно на "хор"):

• Задача 5 — 35 баллов

Сложная часть (необходимо на "отл"):

• Задача 6 — 35 баллов

---

# Bot check

# HW_ID: phds_hw8
# Бот проверит этот ID и предупредит, если случайно сдать что-то не то.

# Status: not final
# Перед отправкой в финальном решении удали "not" в строчке выше.
# Так бот проверит, что ты отправляешь финальную версию, а не промежуточную.
# Никакие значения в этой ячейке не влияют на факт сдачи работы.

import pandas as pd
import numpy as np
import scipy.stats as sps
import matplotlib.pyplot as plt
import seaborn as sns

from scipy.special import gamma, loggamma

sns.set_style("whitegrid")

%matplotlib inline

Теоретическая часть

Легкая часть

Задача 1

Найдите оценку параметра $\theta$ методом максимального правдоподобия по выборке размера $n$ из распределения:

1. Для $\mathcal{N}(a,\sigma^2)$ найдите ОМП в следующих случаях:
   a. $\theta = (a, \sigma^2)$
   b. $\theta = \sigma^2$, $a$ известно
   c. $\theta = a$, $\sigma^2$ известно (для этого случая посчитайте также асимптотическую дисперсию оценки по теореме из лекции)

2. $\mathrm{Pois}(\theta)$. Посчитайте асимптотическую дисперсию, если оценка является асимптотически нормальной по теореме из лекции.

Задача 2

Рассмотрим задачу классификации в случае наличия более двух классов.

Пусть $X_1,...,X_n$ — выборка из категориального распределения, то есть $P_\theta(X_1 = j) = \theta_j$ для $j \in \{1, ..., k\}$, причем $\theta = (\theta_1, ..., \theta_k), \theta_j \geqslant 0$ и $\theta_1 + ... + \theta_k = 1$.

1. Получите функцию правдоподобия для данной задачи. Используя результат, выпишите функцию потерь, которую вы бы минимизировали при обучении модели классификации. Проверьте, что в случае $k=2$ результат совпадает с logloss'ом.
2. Найдите оценку максимального правдоподобия параметра $\theta$ и проверьте ее на состоятельность.

Сложная часть

Задача 3

Рассмотрим модель линейной регрессии

$$Y = X \theta + ɛ,$$

где $Y \in \mathbb{R}^n$ - отклик, $X \in \mathbb{R}^{n \times d}$, $\theta \in \mathbb{R}^d$, $ɛ \in \mathbb{R}^n$ - шум, $n > d$, $rk X = d$.

Будем рассматривать нормальный и гомоскедастичный шум, т.е. $ɛ \sim \mathcal{N}(0, 1)$

1. Получите выражение для функции правдоподобия в данной модели. Минимизации какой функции потерь эквивалента максимизация правдоподобия в данной задаче?

2. Найдите оценку максимального правдоподобия для $(\theta, \sigma^2)$.

3. Пусть $x_{new} \in \mathbb{R}^{d}$ - новый объект. Постройте асимптотический доверительный интервал для ожидаемого значения отклика на этом объекте $y_{new} = x_{new}^T \theta$.
   Указание: используя модель регрессии, получите распределение, которое имеет величина $\widehat{\theta}_{\text{МНК}}$. Учитывая свойства, данные в домашнем задании №1, получите распределение для $x_{new}^T\widehat{\theta}_{\text{МНК}}$. Считая величину $\sigma^2$ известной, запишите интервал. Далее замените дисперсию $\sigma^2$ на ее состоятельную оценку.

Практическая часть

Легкая часть

Задача 4

Cкачайте данные wine dataset и рассмотрите столбцы Alcalinity of ash, Nonflavanoid phenols, Proanthocyanins и Hue для вина первого типа. Тип вина указан в первом столбце. Для работы с табличными данными используйте библиотеку pandas.

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

• асимтотические доверительные интервалы при помощи центральной предельной теоремы;
• точные неасимптотические при помощи распределений хи-квадрат, Стьюдента.

Запишите их в виде таблицы.

# ваш код

Сделайте выводы по полученной таблице.

Вывод:

Задача 5

Функция правдоподобия

Дана параметрическая модель и 3 выборки, состоящие из 2-3 наблюдений. Для удобства, выборки представлены в виде python-кода — каждая выборка записана как список ее элементов; множество выборок представлено как список списков, соответствующих выборкам из множества. Нужно для каждой выборки построить график функции правдоподобия.

a). Параметрическая модель $\mathcal{N}(\theta, 1)$, выборки: [[-1, 1], [-5, 5], [-1, 5]]

b). Параметрическая модель $Exp(\theta)$, выборки: [[1, 2], [0.1, 1], [1, 10]]

c). Параметрическая модель $U[0, \theta]$, выборки: [[0.2, 0.8], [0.5, 1], [0.5, 1.3]]

d). Параметрическая модель $Bin(5, \theta)$, выборки: [[0, 1], [5, 5], [0, 5]]

e). Параметрическая модель $Pois(\theta)$, выборки: [[0, 1], [0, 10], [5, 10]]

f). Параметрическая модель $Сauchy(\theta)$, где $\theta$ — параметр сдвига, выборки: [[-0.5, 0.5], [-2, 2], [-4, 0, 4]]

Выполнить задание, не создавая много кода, поможет следующая функция.

Допишите ее.

def draw_likelihood(density_function, grid, samples, label):
    """Изображает график функции правдоподобия для каждой из 3 выборок.

    Аргументы:
    density_function --- функция, считающая плотность (обычную или дискретную).
        На вход данная функция должна принимать массив размера (1, len_sample)
        и возвращать массив размера (len_grid, len_sample).
    grid --- массив размера (len_grid, 1) --- сетка для построения графика;
    samples --- три выборки;
    label --- latex-код параметрической модели.
    """
    assert len(samples) == 3, "Число выборок не равно 3."

    plt.figure(figsize=(18, 5))
    for i, sample in enumerate(samples):
        sample = np.array(sample)[np.newaxis, :]
        likelihood = <подсчитайте значение правдоподобия, используя density_function и sample>

        plt.subplot(1, 3, i+1)
        plt.plot(grid, likelihood)
        plt.xlabel('$\\theta$', fontsize=16)
        plt.grid(ls=':')
        plt.title(label + ', sample=' + str(sample), fontsize=16)
    plt.show()

Первый пункт можно выполнить с помощью следующего кода:

grid = np.linspace(-5, 5, 1000).reshape((-1, 1))
draw_likelihood(sps.norm(loc=grid).pdf, grid,
                [[-1, 1], [-5, 5], [-1, 5]], '$\\mathcal{N}(\\theta, 1)$')

Выполните остальные:

# ваш код

Прокомментируйте полученные графики. Сделайте вывод о том, как функция правдоподобия для каждой модели зависит от выборки. Является ли функция правдоподобия плотностью?

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

Вывод:

Сгенерируем выборку большого размера из стандартного нормального распределения и посчитаем ее функцию правдоподобия в модели $\mathcal{N}(\theta, 1)$. Выполните код ниже:

sample = sps.norm.rvs(size=10**5)
likelihood = sps.norm.pdf(sample).prod()
print(likelihood)

Почему результат отличается от ожидаемого? Как обойти эту неприятность для подсчета оценки максимального правдоподобия? Реализуйте это.

Подсказка: нужно использовать некоторый метод класса, реализующий это распределение.

Ответ на вопрос и описание метода решения проблемы:

# ваш код

Сложная часть

Задача 6

Асимпотические доверительные интервалы ничего не могут гарантировать на малых размерах выборки. В этой задаче вам предстоит иллюстрировать этот факт, посчитав реальный уровень доверия для интервалов.

Реальный уровень доверия (оценка доли покрытия интервалом) - доля случаев попадания истинного значения параметра в доверительный интервал.

Пример подсчета:

Допустим, вы хотите оценить реальный уровень доверия интервала для $a$, если $X \sim \mathcal{N}(a, \sigma^2)$

• Фиксируете истинные $(a, \sigma^2)$, для которых будете делать оценку
• Генерируете $B$ выборок из $\mathcal{N}(a, \sigma^2)$ с зафиксированными параметрами
• По каждой выборке получаете АДИ для $a$
• Считаете долю случаев, когда истинное $a$ попадает в интервал

Эта доля и будет оценкой доли покрытия интервала.

Важно: вы симулируете реальную ситуацию, когда вы не знаете ни $a$, ни $\sigma$, поэтому в формулах для АДИ их использовать нельзя!

Важно: при такой оценке реального уровня доверия вы используете метод Монте-Карло. Погрешность этого метода составляет $\sim \frac{1}{\sqrt{n}}$, где $n$ - количество выборок, по которым осуществляется оценка.

Вопрос: какое $n$ нужно брать, если вы хотите оценить реальный уровень доверия с точностью до 2 знаков ($\delta = 0.01$)?

Ваш ответ:

Генерация выборок для оценки

Сгенерируйте набор выборок из нормального распределения $\mathcal{N}(\theta, 1)$ при $\theta=0$

theta = 0  # истинное значение параметра
sample_size = 300
sample_count = <...>
X = <...>

Случай АДИ

На лекции вы получали формулу для асимптотического доверительного интервала для $a$

Асимптотический доверительный интервал: $\theta \in \bigg(\overline{X} - \frac{S}{\sqrt{n}}z_{(1+\alpha)/2}, \overline{X} + \frac{S}{\sqrt{n}}z_{(1+\alpha)/2} \bigg)$

Посчитайте $z$ (используйте функцию .ppf)

Постройте график зависимости реального уровня доверия интервала от размера выборки. График начинайте с $n=2$.

Для этого вам нужно провести описанную выше процедуру для подвыборок разных размеров. Для подсчета интервалов по префиксам используйте функцию np.cumsum. Для того чтобы не запутаться, разносите подсчет разных величин, входящих в формулу, по разным строчкам.

Сделайте выводы.

Случай ТДИ

На лекции вы получали формулу для точного доверительного интервала для $a$ в нормальной модели

Точный доверительный интервал: $\theta \in \bigg(\overline{X} - \frac{S}{\sqrt{n-1}}T_{n-1,(1+\alpha)/2}, \overline{X} + \frac{S}{\sqrt{n-1}}T_{n-1,(1+\alpha)/2} \bigg)$

Вопрос: чем этот интервал лучше предыдущего?

Ваш ответ:

Постройте график реального уровня доверия интервала от размера выборки для этого вида интервала. График начинайте с $n=2$. Сравните его с предыдущим.

Вывод:

Меню

• Главная
• Команда
• Учебные курсы
  • Введение в анализ данных
  • DS-поток
  • Phystech@DataScience
  • Семинары по математической статистике на ФБМФ, ФМХФ
• Видеолекции
  • Python для анализа данных
  • Лекции по статистике на DS-потоке

Контакты

• mipt.stats@yandex.ru

При использовании материалов с данного сайта обязательна ссылка на первоисточник.

Design: HTML5 UP.