Phystech@DataScience¶

Домашнее задание 10¶

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

  • Выполненную работу нужно отправить телеграм-боту @thetahat_pds_bot. Для начала работы с ботом каждый раз отправляйте /start. Дождитесь подтверждения от бота, что он принял файл. Если подтверждения нет, то что-то не так. Работы, присланные иным способом, не принимаются.
  • Дедлайн см. в боте. После дедлайна работы не принимаются кроме случаев наличия уважительной причины.
  • Прислать нужно ноутбук в формате ipynb. Если вы строите интерактивные графики, их стоит прислать в формате html.
  • Следите за размером файлов. Бот не может принимать файлы весом более 20 Мб. Если файл получается больше, заранее разделите его на несколько.
  • Выполнять задание необходимо полностью самостоятельно. При обнаружении списывания всем участникам списывания дается штраф -2 балла к итоговой оценке за семестр.
  • Решения, размещенные на каких-либо интернет-ресурсах, не принимаются. Кроме того, публикация решения в открытом доступе может быть приравнена к предоставлении возможности списать.
  • Обратите внимание на правила использования ИИ-инструментов при решении домашнего задания.
  • Код из рассказанных на занятиях ноутбуков можно использовать без ограничений.
  • Для выполнения задания используйте этот ноутбук в качестве основы, ничего не удаляя из него. Можно добавлять необходимое количество ячеек.
  • Комментарии к решению пишите в markdown-ячейках.
  • Выполнение задания (ход решения, выводы и пр.) должно быть осуществлено на русском языке.
  • Решение проверяется системой ИИ-проверки No description has been provided for this image ThetaGrader. Результат проверки валидируется и исправляется человеком, после чего комментарии отправляются студентам.
  • Если код будет не понятен проверяющему, оценка может быть снижена.
  • Никакой код из данного задания при проверке запускаться не будет. Если код студента не выполнен, недописан и т.д., то он не оценивается.

Важно!!! Правила заполнения ноутбука:

  • Запрещается удалять имеющиеся в ноутбуке ячейки, менять местами положения задач.
  • Сохраняйте естественный линейный порядок повествования в ноутбуке сверху-вниз.
  • Отвечайте на вопросы, а также добавляйте новые ячейки в предложенных местах, которые обозначены <...>.
  • В markdown-ячейках, содержащих описание задачи, находятся специальные отметки, которые запрещается модифицировать.
  • При нарушении данных правил работа может получить 0 баллов.

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

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

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

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

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

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

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

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

  • Задача 5 — 25 баллов
In [ ]:
# Bot check

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

# Status: not final
# Перед отправкой в финальном решении удали "not" в строчке выше.
# Так бот проверит, что ты отправляешь финальную версию, а не промежуточную.
# Никакие значения в этой ячейке не влияют на факт сдачи работы.
In [1]:
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¶

Пусть $x_1, \dots, x_n$ — реализация выборки, $X^{*}_1, X^{*}_2, \dots, X^{*}_n$ — построенная по ней бутстрепная выборка. С какой вероятностью элемент $x_i$ исходной выборки попадёт в бутстрепную? Посчитайте среднее число уникальных (попарно различных) элементов в бутстрепной выборке, если в исходной выборке все элементы различны? К какой функции от $n$ стремится это количество при $n \rightarrow \infty$?

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

Задание 2¶

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

$$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{МНК}}$. Учитывая свойства, данные в домашнем задании №6, получите распределение для $x_{new}^T\widehat{\theta}_{\text{МНК}}$. Считая величину $\sigma^2$ известной, запишите интервал. Далее замените дисперсию $\sigma^2$ на ее состоятельную оценку.

  4. Также постройте точный доверительный интервал для ожидаемого значения отклика на том же новом объекте $x_{new}$.

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

Легкая часть¶

Задание 3¶

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

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

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

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

In [35]:
# ваш код

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

Вывод:

Задание 4¶

Профиль Физика

Скачайте данные столкновениях частиц и оставьте следующие признаки:

  • E1, E2 — полная энергия электронов (ГэВ);
  • M — инвариантная масса двух электронов (ГэВ).

Датасет большой, поэтому для ускорения расчетов можно взять подвыборку размера ~ 1000

Профиль биология

Скачайте данные экспрессии белков у белых мышей и оставьте следующие белки для анализа:

  • ITSN1_N
  • DYRK1A_N
  • ELK_N

1. Постройте для каждого из них гистограммы. Что можно сказать о характере распределения?

In [ ]:
df = pd.read_csv("your/data/path")

Ответ:

2. Пусть есть выборка $X_1, ..., X_n$. Опишите процедуру бутстрепа построения доверительного интервала для величины $\theta = \mathsf{E} X_1$. Рассмотрите три типа бутстрепных доверительных интервалов. Выпишите соответствующие формулы.

Описание:

3. Реализуйте функцию вычисления выборки оценок $\theta^*_1, \ldots \theta^*_B$ для оценки $\widehat \theta$ методом бутстрепа. Функция должна принимать на вход выборку и количество бутстрепных выборок $B$.

Для ускорения расчетов можете использовать broadcasting (не является обязательным требованием). Учите, что в таком случае внутри функции у вас могут получиться очень большие матрицы, из-за которых может кончиться оперативная память.

In [ ]:
def bootstrap(sample, B=100_000): 
    '''
    Считает бутстрепные оценки для исходной выборки

    :param sample: исходная выборка
    :param B: количество бутстрепных 
    :return bootstrap_estimations: оценки по бутстрепным выборкам
    '''   
    
    return

4. Реализуйте три типа бутстрепных доверительных интервалы в виде функций, принимающих на вход выборку оценок, полученных с помощью бутстрепа. Реализовывать вычисление бутстрепных интервалов для подвыборок размера от 1 до n не нужно, только для всей выборки.

In [ ]:
def bootstrap_normal_confidence_interval(theta, bootstrap_estimations, alpha=0.95):
    '''
    Считает левую и правую границу нормального бутстрепного интервала

    :param theta: оценка параметра
    :param bootstrap_estimations: массив бутстрепных оценок
    :return left: левая граница бутстрепного интервала
    :return right: правая граница бутстрепного интервала
    '''   

    return

    
def bootstrap_central_confidence_interval(theta, bootstrap_estimations, alpha=0.95):
    '''
    Считает левую и правую границу центрального бутстрепного интервала

    :param theta: оценка параметра
    :param bootstrap_estimations: массив бутстрепных оценок
    :return left: левая граница бутстрепного интервала
    :return right: правая граница бутстрепного интервала
    '''   

    return

    
def bootstrap_quantile_confidence_interval(theta, bootstrap_estimations, alpha=0.95):
    '''
    Считает левую и правую границу квантильного бутстрепного интервала

    :param theta: оценка параметра
    :param bootstrap_estimations: массив бутстрепных оценок
    :return left: левая граница бутстрепного интервала
    :return right: правая граница бутстрепного интервала
    '''   
    
    return

5. Для каждого признака постройте бутстрепные доверительные интервалы для $\theta = \mathsf{E} X_1$ и сравните их.

In [ ]:
data = pd.DataFrame(index = ['normal', 'central', 'quantile', 'theta', 'interval_length'], columns=df.columns)

Вывод:

Визуализируйте бутстрепные интервалы для каждого признака. Для этого сгенерируйте выборку $X_1, ... X_{N}, N = 100$ и постройте график доверительных интервалов уровня доверия $0.95$, вычисленных для всех подвыборок размера $n$ вида $X_1, ... X_n$, $1 \le n \le 100$, используя написанную ниже функцию

In [ ]:
def draw_confidence_interval(
    left,
    right,
    estimation=None,
    sample=None, 
    ylim=(None, None), 
    estim_label = '',
    sample_label='',
    color=None,
    interval_label=None
):
    '''
    Рисует доверительный интервал и оценку в зависимости от размера выборки.
    
    :param left: левые границы интервалов (в зависимости от n)
    :param right: правые границы интервалов (в зависимости от n)
    :param estimation: оценки (в зависимости от n)
    :param sample: выборка
    :param ylim: ограничение вертикальной оси
    :param estim_label: подпись к оценке
    :param sample_label: подпись к выборке
    :param color: цвет, которым будет отображен доверительный интервал
    '''
    
    return

Решение:

In [ ]:

Вывод:

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

Задача 5¶

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

Допустим, вы хотите оценить реальный уровень доверия интервала для $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$)?

Ваш ответ:

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

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

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

2. Построение ТДИ

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

Точный доверительный интервал: $\theta \in $ <...>

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

Ваш ответ:

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

In [ ]:

Вывод: