В данном задании предлагается провести некоторое исследование модели линейной регрессии и критериев для проверки статистических гипотез, в частности применить этим модели к реальным данным.
Правила:
probability.diht@yandex.ru, указав тему письма "[номер группы] Фамилия Имя - Задание 5". Квадратные скобки обязательны. Вместо Фамилия Имя нужно подставить свои фамилию и имя.5.N.ipynb и 5.N.pdf, где N - ваш номер из таблицы с оценками.Баллы за задание:
In [1]:
import numpy as np
import scipy.stats as sps
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
Задача 1. По шаблону напишите класс, реализующий линейную регрессию. Интерфейс этого класса в некоторой степени соответствует классу `LinearRegression` из библиотеки sklearn.
In [2]:
from scipy.linalg import inv
from numpy.linalg import norm
In [3]:
class LinearRegression:
def __init__(self):
super()
def fit(self, X, Y, alpha=0.95):
''' Обучение модели. Предполагается модель Y = X * theta + epsilon,
где X --- регрессор, Y --- отклик,
а epsilon имеет нормальное распределение с параметрами (0, sigma^2 * I_n).
alpha --- уровень доверия для доверительного интервала.
'''
# Размер выборки и число признаков
self.n, self.k = X.shape
# Оценки на параметры
self.theta = inv(X.T @ X) @ X.T @ Y
self.sigma_sq = norm(Y - X @ self.theta) ** 2 / (self.n - self.k)
# Считаем доверительные интервалы
l_quant = sps.t.ppf((1 - alpha) / 2, df=self.n - self.k)
r_quant = sps.t.ppf((1 + alpha) / 2, df=self.n - self.k)
diag = inv(X.T @ X).diagonal()
coeff = np.sqrt(self.sigma_sq * diag)
self.conf_int = np.array([self.theta + l_quant * coeff, self.theta + r_quant * coeff]).T
return self
def summary(self):
print('Linear regression on %d features and %d examples' % (self.k, self.n))
print('Sigma: %.6f' % self.sigma_sq)
print('\t\tLower\t\tEstimation\tUpper')
for j in range(self.k):
print('theta_%d:\t%.6f\t%.6f\t%.6f' % (j, self.conf_int[j, 0],
self.theta[j], self.conf_int[j, 1]))
def predict(self, X):
''' Возвращает предсказание отклика на новых объектах X. '''
Y_pred = X @ self.theta
return Y_pred
Загрузите данные о потреблении мороженного в зависимости от температуры воздуха и цены (файл ice_cream.txt).
Примените реализованный выше класс линейной регрессии к этим данным предполагая, что модель имеет вид $ic = \theta_1 + \theta_2\ t$, где $t$ --- температура воздуха (столбец temp), $ic$ --- постребление мороженного в литрах на человека (столбец IC).
Значения температуры предварительно переведите из Фаренгейта в Цельсий [(Фаренгейт — 32) / 1,8 = Цельсий].
К обученной модели примените фунцию summary и постройте график регрессии, то есть график прямой $ic = \widehat{\theta}_1 + \widehat{\theta}_2\ t$, где $\widehat{\theta}_1, \widehat{\theta}_2$ --- МНК-оценки коэффициентов.
На график нанесите точки выборки.
Убедитесь, что построейнный график совпадает с графиком из презентации с первой лекции, правда, с точностью до значений температура (она была неправильно переведена из Фаренгейта в Цельсий).
In [4]:
import csv
In [5]:
# Я загурзил файл локально
data = np.array(list(csv.reader(open('ice_cream.txt', 'r'), delimiter='\t')))
source = data[1:, :].astype(float)
print(source[:5])
In [6]:
# Переводим в Цельсий
source[:, 4] = (source[:, 4] - 32) / 1.8
source[:, 5] = (source[:, 5] - 32) / 1.8
In [7]:
# Размер выборки и число параметров
n, k = source.shape[0], 2
print(n)
In [8]:
# Отклик и регрессор
Y = source[:, 1]
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4]
print(X[:5])
In [9]:
# Обучаем модель
model = LinearRegression()
model.fit(X, Y)
Out[9]:
In [10]:
# Выводим общую информацию
model.summary()
In [11]:
grid = np.linspace(-5, 25, 1000)
In [12]:
plt.figure(figsize=(20, 8))
plt.plot(grid, model.theta[0] + grid * model.theta[1], color='brown', label='Предсказание', linewidth=2.5)
plt.scatter(X[:, 1], Y, s=40.0, label='Выборка', color='red', alpha=0.5)
plt.legend()
plt.ylabel('Литров на человека')
plt.xlabel('Температура')
plt.title('Потребление мороженого')
plt.grid()
plt.show()
Вывод. Действительно, график тот же (с поправкой на пересчет температуры). Линейная регрессия неплохо, но не идеально приближает зависимость потребления мороженого в зависимости от температуры. Для более точного вывода стоит посчитать доверительный интервал.
Теперь учтите влияние года (столбец Year) для двух случаев:
В каждом случае нарисуйте графики. Отличаются ли полученные результаты? От чего это зависит? Как зависит потребление мороженного от года?
Первый случай.
In [13]:
# Размер выборки и число параметров
n, k = source.shape[0], 4
In [14]:
# Отклик и регрессор
Y = source[:, 1]
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4]
X[:, 2] = (source[:, 6] == 1).astype(int)
X[:, 3] = (source[:, 6] == 2).astype(int)
print(X[:5])
In [15]:
# Обучаем модель
model = LinearRegression()
model.fit(X, Y)
Out[15]:
In [16]:
# Выводим общую информацию
model.summary()
In [17]:
grid = np.linspace(-5, 25, 1000)
In [18]:
y_0 = model.theta[0] + grid * model.theta[1]
y_1 = model.theta[0] + grid * model.theta[1] + model.theta[2]
y_2 = model.theta[0] + grid * model.theta[1] + model.theta[3]
In [19]:
plt.figure(figsize=(20, 8))
plt.plot(grid, y_0, color='gold', label='Предсказание (год 0)', linewidth=2.5)
plt.plot(grid, y_1, color='turquoise', label='Предсказание (год 1)', linewidth=2.5)
plt.plot(grid, y_2, color='springgreen', label='Предсказание (год 2)', linewidth=2.5)
plt.scatter(X[source[:, 6] == 0, 1], Y[source[:, 6] == 0], s=40.0, label='Выборка (год 0)', color='yellow', alpha=0.5)
plt.scatter(X[source[:, 6] == 1, 1], Y[source[:, 6] == 1], s=40.0, label='Выборка (год 1)', color='blue', alpha=0.5)
plt.scatter(X[source[:, 6] == 2, 1], Y[source[:, 6] == 2], s=40.0, label='Выборка (год 2)', color='green', alpha=0.5)
plt.legend()
plt.ylabel('Литров на человека')
plt.xlabel('Температура')
plt.title('Потребление мороженого')
plt.grid()
plt.show()
Вывод. В разрезе разных лет также можно уследить, что линейная зависимсоть неплохо приближает реальную. Номер года нельзя брать как признак, так это больше характеризация класса каких-то значений, в то время линейная регрессия работает со значениями, от которых зависимость линейная, численная
Второй случай.
In [20]:
# Число признаков
k = 2
# Размеры подвыборок
n_0 = (source[:, 6] == 0).sum()
n_1 = (source[:, 6] == 1).sum()
n_2 = (source[:, 6] == 2).sum()
print(n_0, n_1, n_2)
In [21]:
# Три подвыборки
source_0, source_1, source_2 = np.vsplit(source, [n_0, n_0 + n_1])
print(source_0.astype(int)[:3]) # Кастим к инт, чтобы при выводе не размазывались по всей строке
print(source_1.astype(int)[:3])
print(source_2.astype(int)[:3])
In [22]:
# Отклики и регрессоры
Y_0 = source_0[:, 1]
X_0 = np.zeros((n_0, k))
X_0[:, 0] = np.ones(n_0)
X_0[:, 1] = source_0[:, 4]
print(X_0[:3])
Y_1 = source_1[:, 1]
X_1 = np.zeros((n_1, k))
X_1[:, 0] = np.ones(n_1)
X_1[:, 1] = source_1[:, 4]
print(X_1[:3])
Y_2 = source_2[:, 1]
X_2 = np.zeros((n_2, k))
X_2[:, 0] = np.ones(n_2)
X_2[:, 1] = source_2[:, 4]
print(X_2[:3])
In [23]:
# Обучаем модели
model_0 = LinearRegression()
model_0.fit(X_0, Y_0)
model_1 = LinearRegression()
model_1.fit(X_1, Y_1)
model_2 = LinearRegression()
model_2.fit(X_2, Y_2)
Out[23]:
In [24]:
# Выводим общую информацию
model_0.summary()
model_1.summary()
model_2.summary()
In [25]:
grid = np.linspace(-5, 25, 1000)
In [26]:
plt.figure(figsize=(20, 8))
plt.plot(grid, model_0.theta[0] + grid * model_0.theta[1], color='gold', label='Предсказание (год 0)', linewidth=2.5)
plt.plot(grid, model_1.theta[0] + grid * model_1.theta[1], color='turquoise', label='Предсказание (год 1)', linewidth=2.5)
plt.plot(grid, model_2.theta[0] + grid * model_2.theta[1], color='springgreen', label='Предсказание (год 2)', linewidth=2.5)
plt.scatter(X[source[:, 6] == 0, 1], Y[source[:, 6] == 0], s=40.0, label='Выборка (год 0)', color='yellow', alpha=0.5)
plt.scatter(X[source[:, 6] == 1, 1], Y[source[:, 6] == 1], s=40.0, label='Выборка (год 1)', color='blue', alpha=0.5)
plt.scatter(X[source[:, 6] == 2, 1], Y[source[:, 6] == 2], s=40.0, label='Выборка (год 2)', color='green', alpha=0.5)
plt.legend()
plt.ylabel('Литров на человека')
plt.xlabel('Температура')
plt.title('Потребление мороженого')
plt.grid()
plt.show()
Вывод. При разделении выборки на три части результаты в какой-то степени стали лучше. Кроме того, такое разделение позволяет заметить, что с при увеличении года в целом потребляют больше мороженого, и что прирост литров мороженого на один градус возрастает при увеличении года.
Наконец, обучите модель на предсказание потребления мороженного в зависимости от всех переменных.
Не забудьте, что для года нужно ввести две переменных.
Для полученной модели выведите summary.
In [27]:
# Размер выборки и число параметров
n, k = source.shape[0], 8
print(n, k)
In [28]:
# Cтроим регрессор
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4] # Температура
X[:, 2] = source[:, 0] # Дата
X[:, 3:5] = source[:, 2:4] # Пропускаем IC
X[:, 5] = source[:, 5]
X[:, 6] = (source[:, 6] == 1).astype(int) # Индикатор год 1
X[:, 7] = (source[:, 6] == 2).astype(int) # Индикатор год 2
print(X.astype(int)[:5])
# Отлкик
Y = source[:, 1]
In [29]:
# Обучаем модель
model = LinearRegression()
model.fit(X, Y)
Out[29]:
In [30]:
# Выводим общую информацию
model.summary()
Вывод. Похоже, не все признаки стоит учитывать. Некоторые значения по модулю очень маленькие по сравнению со значениями, характерными для соответствующих признаков, что их лучше не учитывать. Например, это номер измерения (theta_2) или температура в будущем месяце.(theta_5).
Но это еще не все. Постройте теперь линейную регрессию для модели $ic = \theta_1 + \theta_2\ t + \theta_3\ t^2 + \theta_4\ t^3$. Выведите для нее summary и постройте график предсказания, то есть график кривой $ic = \widehat{\theta}_1 + \widehat{\theta}_2\ t + \widehat{\theta}_3\ t^2 + \widehat{\theta}_4\ t^3$. Хорошие ли получаются результаты?
In [31]:
# Размер выборки и число параметров
n, k = source.shape[0], 4
print(n, k)
In [32]:
# Отклик и регрессор
Y = source[:, 1]
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4]
X[:, 2] = X[:, 1] ** 2
X[:, 3] = X[:, 1] ** 3
print(X[:5])
In [33]:
# Обучаем модель
model = LinearRegression()
model.fit(X, Y)
Out[33]:
In [34]:
# Выводим общую информацию
model.summary()
In [35]:
grid = np.linspace(-5, 25, 1000)
In [36]:
y = model.theta[0] + model.theta[1] * grid + model.theta[2] * grid ** 2 + model.theta[3] * grid ** 3
In [37]:
plt.figure(figsize=(20, 8))
plt.plot(grid, y, color='brown', label='Предсказание (год 0)', linewidth=2.5)
plt.scatter(X[:, 1], Y, s=40.0, label='Выборка', color='red', alpha=0.5)
plt.legend()
plt.ylabel('Литров на человека (первый случай)')
plt.xlabel('Температура')
plt.title('Потребление мороженого')
plt.grid()
plt.show()
Вывод. Результаты выглядят более естественно, однако теперь кажется, что при 30-40 градусах люди потребляют невероятно много мороженого, что, скорее всего, неправда (все-таки есть какой-то лимит). Кроме того, значения последних двух параметров очень малы, что говорит о малом влиянии этих параметров при малых температурах. Возможно, линейная модель и так достаточно хороша. Однако в нашем случае будет видно, что при отрицательных температурах потребление мороженого падает стремительно. Наверное, в этом есть некоторая правда.
Чтобы понять, почему так происходит, выведите значения матрицы $(X^T X)^{-1}$ для данной матрицы и посчитайте для нее индекс обусловленности $\sqrt{\left.\lambda_{max}\right/\lambda_{min}}$, где $\lambda_{max}, \lambda_{min}$ --- максимальный и минимальный собственные значения матрицы $X^T X$. Собственные значения можно посчитать функцией `scipy.linalg.eigvals`.
Прокомментируйте полученные результаты. Помочь в этом может следующая статья.
In [38]:
D = inv(X.T @ X)
print(D)
In [39]:
from scipy.linalg import eigvals
In [40]:
vals = eigvals(D)
print(vals) # Комплексных нет, мы не налажали
vals = vals.real.astype(float)
print(vals)
In [41]:
CI = np.sqrt(vals.max() / vals.min())
print(CI)
Вывод. Как говорилось на семинаре, высокий индекс обусловленности (больше 30) говорит о мультиколлинеарности, которая ведет к переобучению. Видимо, мы перестарались.
Задача 2. В данной задаче нужно реализовать функцию отбора признаков для линейной регрессии. Иначе говоря, пусть есть модель $y = \theta_1 x_1 + ... + \theta_k x_k$. Нужно определить, какие $\theta_j$ нужно положить равными нулю, чтобы качество полученной модели было максимальным.
Для этого имеющиеся данные нужно случайно разделить на две части --- обучение и тест (train и test). На первой части нужно обучить модель регресии, взяв некоторые из признаков, то есть рассмотреть модель $y = \theta_{j_1} x_{j_1} + ... + \theta_{j_s} x_{j_s}$. По второй части нужно посчитать ее качество --- среднеквадратичное отклонение (mean squared error) предсказания от истинного значения отклика, то есть величину $$MSE = \sum\limits_{i \in test} \left(\widehat{y}(x_i) - Y_i\right)^2,$$ где $x_i = (x_{i,1}, ..., x_{i,k})$, $Y_i$ --- отклик на объекте $x_i$, а $\widehat{y}(x)$ --- оценка отклика на объекте $x$.
Если $k$ невелико, то подобным образом можно перебрать все поднаборы признаков и выбрать наилучший по значению MSE.
Для выполнения задания воспользуйтесь следующими функциями:
`sklearn.linear_model.LinearRegression`
--- реализация линейной регрессии. В данной реализации свободный параметр $\theta_1$ по умолчанию автоматически включается в модель. Отключить это можно с помощью fit_intercept=False, но это не нужно. В данной задаче требуется, чтобы вы воспользовались готовой реализацией линейной регрессии, а не своей. Ведь на практике важно уметь применять готовые реализации, а не писать их самостоятельно.
`sklearn.cross_validation.train_test_split`
--- функция разбиения данных на train и test. Установите параметр test_size=0.3.
`sklearn.metrics.mean_squared_error` --- реализация MSE.
Для перебора реализуйте функцию.
In [42]:
from sklearn import linear_model
from sklearn import cross_validation
from sklearn.metrics import mean_squared_error
In [43]:
def best_features(X_train, X_test, Y_train, Y_test):
mses = [] # сюда записывайте значения MSE
k = X_train.shape[1]
for j in range(1, 2 ** k): # номер набора признаков
mask = np.array([j & (1 << s) for s in range(k)], dtype=bool)
features_numbers = np.arange(k)[mask] # набор признаков
model = linear_model.LinearRegression()
model.fit(X_train[:, features_numbers], Y_train)
mse = mean_squared_error(Y_test, model.predict(X_test[:, features_numbers])) # MSE для данного набора признаков
mses.append(mse)
# Печать 10 лучших наборов
print('mse\t features')
mses = np.array(mses)
best_numbres = np.argsort(mses)[:10]
for j in best_numbres:
mask = np.array([j & (1 << s) for s in range(k)], dtype=bool)
features_numbers = np.arange(k)[mask]
print('%.3f\t' % mses[j], features_numbers)
Примените реализованный отбор признаков к датасетам
Yacht Hydrodynamics --- для парусных яхт нужно оценить остаточное сопротивление на единицу массы смещения (последний столбец) в зависимости от различных характеристик яхты.
Boston Housing Prices --- цены на дома в Бостоне в зависимости от ряда особенностей.
Yacht Hydrodynamics.
In [44]:
# Я загурзил файл локально и исправил проблелы на табы (там были лишние проблеы в некоторых местах)
data = np.array(list(csv.reader(open('yh.data', 'r'), delimiter='\t')))
yacht = data.astype(float)
print(yacht)
In [45]:
Y = yacht[:, 6]
X = yacht[:, :6]
In [46]:
X_train, X_test, Y_train, Y_test = cross_validation.train_test_split(X, Y, test_size=0.3)
In [47]:
best_features(X_train, X_test, Y_train, Y_test)
Вывод. Признак 5 (Froude number) наиболее полезный признак, остальные встречаются как-то рандомно (кроме 4, который почти не встречается, это Length-beam ratio).
Boston Housing Prices.
In [48]:
# Я загурзил файл локально и убрал первую стрчоку с кучей запятых
data = np.array(list(csv.reader(open('bhp.csv', 'r'))))
houses = data[1:, :].astype(float)
print(houses)
In [49]:
# Число столбцов
cols = houses.shape[1]
In [50]:
Y = houses[:, 9] # Столбец TAX
X = houses[:, np.delete(np.arange(cols), 9)]
In [51]:
X_train, X_test, Y_train, Y_test = cross_validation.train_test_split(X, Y, test_size=0.3)
In [52]:
best_features(X_train, X_test, Y_train, Y_test)
Вывод. Первые 3 и последние 3 признака самые полезные.
Задача 3\*. Загрузите датасет, в котором показана зависимость веса мозга от веса туловища для некоторых видов млекопитающих. Задача состоит в том, чтобы подобрать по этим данным хорошую модель регрессии. Для этого, можно попробовать взять некоторые функции от значения веса туловища, например, степенную, показательную, логарифмическую. Можно также сделать преобразование значений веса мозга, например, прологарифмировать. Кроме того, можно разбить значения веса туловища на несколько частей и на каждой части строить свою модель линейной регрессии.
Задача 4. Пусть $X_1, ..., X_n$ --- выборка из распределения $\mathcal{N}(a, \sigma^2)$. Постройте точную доверительную область для параметра $\theta = (a, \sigma^2)$ уровня доверия $\alpha=0.95$ для сгенерированной выборки размера $n \in \{5, 20, 50\}$ из стандартного нормального распределения. Какой вывод можно сделать?
Вспомним, что $T=\frac{\overline{X}-a}{\frac{S}{\sqrt{n}}} \sim T_{n-1}$. Положим $u_{\frac{1\pm\sqrt{\alpha}}{2}}$ - квантили $T_{n-1}$. Тогда $\sqrt{\alpha}=P\left(u_{\frac{1-\sqrt{\alpha}}{2}} < T < u_{\frac{1+\sqrt{\alpha}}{2}}\right)=P\left(\overline{X} - \frac{S}{\sqrt{n}} u_{\frac{1+\sqrt{\alpha}}{2}} < a < \overline{X} - \frac{S}{\sqrt{n}} u_{\frac{1-\sqrt{\alpha}}{2}}\right)$, то есть доверительный интервал уровня доверия $\sqrt{\alpha}$ для $a$ есть $$\left(\overline{X} - \frac{S}{\sqrt{n}} u_{\frac{1+\sqrt{\alpha}}{2}}, \overline{X} - \frac{S}{\sqrt{n}} u_{\frac{1-\sqrt{\alpha}}{2}}\right).$$
Вспомним, что $R=\frac{n-1}{\sigma^2}S^2 \sim \chi^2_{n-1}$. Пусть $v_{\frac{1\pm\sqrt{\alpha}}{2}}$ - квантили $\chi^2_{n-1}$. Тогда $\sqrt{\alpha}=P\left(v_{\frac{1-\sqrt{\alpha}}{2}} < R < v_{\frac{1+\sqrt{\alpha}}{2}}\right)=P\left(\frac{(n-1)S^2}{v_{\frac{1+\sqrt{\alpha}}{2}}} < \sigma^2 < \frac{(n-1)S^2}{v_{\frac{1-\sqrt{\alpha}}{2}}}\right)$, то есть доверительный интервал уровня доверия $\sqrt{\alpha}$ для $\sigma^2$ есть $$\left(\frac{(n-1)S^2}{v_{\frac{1+\sqrt{\alpha}}{2}}}, \frac{(n-1)S^2}{v_{\frac{1-\sqrt{\alpha}}{2}}}\right).$$
Из домашнего задания $\overline{X}$ и $S^2$ независимы, значит, $T$ и $R$ тоже независимы. Таким образом, $\alpha=P(u_{\ldots} < T < u_{\ldots}, v_{\ldots} < R < v_{\ldots})$, откуда декартово произведение полученных доверительных интервалов уровней доверия $\sqrt{\alpha}$ есть доверительная область уровня доверия $\alpha$ для $(a, \sigma^2)$ .
In [53]:
def find_conf_reg(sample):
size = sample.size
alpha_r = np.sqrt(0.95)
u_1 = sps.t.ppf((1 - alpha_r) / 2, df=size - 1)
u_2 = sps.t.ppf((1 + alpha_r) / 2, df=size - 1)
v_1 = sps.chi2.ppf((1 - alpha_r) / 2, df=size - 1)
v_2 = sps.chi2.ppf((1 + alpha_r) / 2, df=size - 1)
mean = sample.mean()
std = sample.std()
a_low = mean - u_2 * (std / size) ** 0.5
a_high = mean - u_1 * (std / size) ** 0.5
s_low = (size - 1) * std / v_2
s_hight = (size - 1) * std / v_1
return ((a_low, a_high), (s_low, s_hight))
In [54]:
plt.figure(figsize=(20, 8))
for size in [5, 20, 50]:
a_conf, s_conf = find_conf_reg(sps.norm.rvs(size=size))
plt.hlines(s_conf[0], a_conf[0], a_conf[1], linewidth=2.5, color='tomato')
plt.hlines(s_conf[1], a_conf[0], a_conf[1], linewidth=2.5, color='tomato')
plt.vlines(a_conf[0], s_conf[0], s_conf[1], linewidth=2.5, color='tomato')
plt.vlines(a_conf[1], s_conf[0], s_conf[1], linewidth=2.5, color='tomato')
plt.ylabel('Выборочная дисперсия')
plt.xlabel('Выборочное среднее')
plt.title('Доверительная область')
plt.grid()
plt.show()
Вывод. При увеличении размера выборки площадь доверительной область (она, очевидно, измерима даже по Жордану) сильно уменьшается. Это говорит о том, что чем больше значений получено, тем точнее с заданным уровнем доверия можно оценить параметр.
Задача 5\*. Пусть дана линейная гауссовская модель $Y = X\theta + \varepsilon$, где $\varepsilon \sim \mathcal{N}(0, \beta^{-1}I_n)$. Пусть $\theta$ имеет априорное распределение $\mathcal{N}(0, \alpha^{-1}I_k)$. Такая постановка задачи соответствует Ridge-регрессии. Оценкой параметров будет математическое ожидание по апостериорному распределению, аналогично можно получить доверительный интервал. Кроме того, с помощью апостериорного распределения можно получить доверительный интервал для отклика на новом объекте, а не только точечную оценку.
Реализуйте класс RidgeRegression подобно классу LinearRegression, но добавьте в него так же возможность получения доверительного интервала для отклика на новом объекте. Примените модель к некоторых датасетам, которые рассматривались в предыдущих задачах. Нарисуйте графики оценки отклика на новом объекте и доверительные интервалы для него.
Задача 6. Существует примета, что если перед вам дорогу перебегает черный кот, то скоро случится неудача. Вы же уже достаточно хорошо знаете статистику и хотите проверить данную примету. Сформулируем задачу на математическом языке. Пусть $X_1, ..., X_n \sim Bern(p)$ --- проведенные наблюдения, где $X_i = 1$, если в $i$-м испытании случилась неудача после того, как черный кот перебежал дорогу, а $p$ --- неизвестная вероятность такого события. Нужно проверить гипотезу $H_0: p=1/2$ (отсутствие связи между черным котом и неудачей) против альтернативы $H_1: p>1/2$ (неудача происходит чаще если черный кот перебегает дорогу).
Известно, что $S = \left\{T(X) > c_\alpha\right\}$, где $T(X) = \sum X_i$, является равномерно наиболее мощным критерием для данной задачи. Чему при этом равно $c_\alpha$? При этом p-value в данной задаче определяется как $p(t) = \mathsf{P}_{0.5}(T(X) > t)$, где $t = \sum x_i$ --- реализация статистики $T(X)$.
Для начала проверьте, что критерий работает. Возьмите несколько значений $n$ и реализаций статистики $T(X)$. В каждом случае найдите значение $c_\alpha$ и p-value. Оформите это в виде таблицы.
Пользуйтесь функциями из scipy.stats, про которые подробно написано в файле python_5. Внимательно проверьте правильность строгих и нестрогих знаков.
In [55]:
alpha = 0.05 # Уровень значимости
In [56]:
# Считаем параметры для выборок четырех разных размеров
stats = np.zeros((4, 4))
for i, size in enumerate([5, 15, 30, 50]):
t = sps.bernoulli(p=0.5).rvs(size=size).sum() # Сумма бернуллиевских случайных величин
pvalue = sps.binom(n=size, p=0.5).sf(t) # Статистика T имеет биномиальное распределение
c_alpha = sps.binom(n=size, p=0.5).ppf(1 - alpha)
stats[i, :] = np.array([size, t, c_alpha, pvalue])
pd.DataFrame(data=stats, columns=['$n$', '$t$', '$c_{\\alpha}$', 'p-value'])
Out[56]:
Вывод. Если $t/n$ сильно больше 0.5 (скажем, 0.67), то p-value выходит меньше 0.05, и гипотеза отвергается. Критерий работает.
Для каких истинных значений $p$ с точки зрения практики можно считать, что связь между черным котом и неудачей есть? Теперь сгенерируйте 10 выборок для двух случаев: 1). $n=5, p=0.75$; 2). $n=10^5, p=0.51$. В каждом случае в виде таблицы выведите реализацию статистики $T(X)$, соответствующее p-value и 0/1 - отвергается ли $H_0$ (выводите 1, если отвергается). Какие выводы можно сделать?
In [57]:
# В двух случаях строим для 10 выборок табличку
for size, p in [(5, 0.75), (100000, 0.51)]:
stats = np.zeros((10, 5))
for i in np.arange(10):
t = sps.bernoulli(p=p).rvs(size=size).sum()
pvalue = sps.binom(n=size, p=0.5).sf(t)
c_alpha = sps.binom(n=size, p=0.5).ppf(1 - alpha)
rejected = int(t > c_alpha)
stats[i, :] = np.array([size, p, t, pvalue, rejected])
print(pd.DataFrame(data=stats, columns=['size', 'prob', 'stat', 'p-value', 'rej']))
Вывод. Выходит, почти всегда на малых выборках даже при большом $p$ мы гипотезу не отвергнем, а при больших даже при $p$, близких к 0.5, гипотеза будет отвергнута. Похоже на ошибки II и I рода соответственно.
Возникает задача подбора оптимального размера выборки.
Для этого сначала зафиксируйте значение $p^* > 1/2$, которое будет обладать следующим свойством. Если истинное $p > p^*$, то такое отклонение от $1/2$ с практической точки зрения признается существенным, то есть действительно чаще случается неудача после того, как черный кот перебегает дорогу. В противном случае отклонение с практической точки зрения признается несущественным.
Теперь для некоторых $n$ постройте графики функции мощности критерия при $1/2 < p < 1$ и уровне значимости 0.05. Выберите такое $n^*$, для которого функция мощности дает значение 0.8 при $p^*$. Для выбранного $n^*$ проведите эксперимент, аналогичный проведенным ранее экспериментам, сгенерировав выборки для следующих истинных значений $p$: 1). $1/2 < p < p^*$; 2). $p > p^*$. Сделайте вывод.
Выберем $p^*=0.75$.
In [58]:
p_conf = 0.75 # Если больше, то уверены в том, что действительно неудачи связаны с черной кошкой.
beta = 0.8 # Магическая константа с семинара и из условия
In [59]:
grid = np.linspace(0.501, 0.999, 500)
In [60]:
plt.figure(figsize=(20, 8))
for size in [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]:
power = sps.binom(n=size, p=grid).sf(sps.binom(n=size, p=0.5).ppf(1 - alpha))
plt.plot(grid, power, label='n={}'.format(size))
plt.legend()
plt.vlines(p_conf, 0, 1)
plt.hlines(beta, 0.5, 1)
plt.title('Мощность критериев')
plt.ylabel('$\\beta(\\theta, X)$')
plt.xlabel('$\\theta$')
plt.grid()
plt.show()
Вывод. Оптимальный в некотором смысле размер выборки получается примерно 25–30, на пересечении $\theta=p^*=0.75$ и $\beta(\theta, X)=\beta=0.8$. Выберем 27.
In [61]:
size_conf = 27 # Оптимальные размер выборки
In [62]:
# В двух случаях строим для 10 выборок табличку
for p in [0.6, 0.85]:
stats = np.zeros((10, 6))
for i in np.arange(10):
t = sps.bernoulli(p=p).rvs(size=size_conf).sum()
pvalue = sps.binom(n=size_conf, p=0.5).sf(t)
c_alpha = sps.binom(n=size_conf, p=0.5).ppf(1 - alpha)
rejected = int(t > c_alpha)
stats[i, :] = np.array([size_conf, p, t, c_alpha, pvalue, rejected])
print(pd.DataFrame(data=stats, columns=['size', 'prob', 'stat', 'c_alpha', 'p-value', 'rej']))
Вывод. При выбранном значении $p^*$ и подобранном оптималном значении выборки если $p < p^*$, то гипотеза не отвергается почти всегда, а при $p > p^*$ гипотеза отвергается почти всегда.
Критерий согласия хи-квадрат
`scipy.stats.chisquare`(f_obs, f_exp=None, ddof=0)
f_obs --- число элементов выборки, попавших в каждый из интервалов
f_exp --- ожидаемое число элементов выборки (по умолчанию равномерное)
ddof --- поправка на число степеней свободы. Статистика асимптотически будет иметь распределение хи-квадрат с числом степеней свободы $k - 1 - ddof$, где $k$ --- число интервалов.
Возвращает значение статистики критерия и соответствующее p-value.
Критерий согласия Колмогорова
`scipy.stats.kstest`(rvs, cdf, args=())
rvs --- выборка
cdf --- функция распределения (сама функция или ее название)
args --- параметры распределения
Возвращает значение статистики критерия и соответствующее p-value.
Задача 7.
Проверьте, что ваша выборка значений скорости ветра из задания 2 действительно согласуется с распределением Вейбулла.
Проверьте, что при больших $n$ распределение статистики из задач 3 и 4 задания 2 действительно хорошо приближают предельное распределение.
Проверьте, что остатки в регрессии из задач выше нормальны.
Подберите класс распределений для выборки количества друзей из задания 1.
Использовать можно два описанных выше критерия, либо любой другой критерий, если будет обоснована необходимость его применения в данной задаче, а так же будет приведено краткое описание критерия. Уровень значимости взять равным 0.05.
In [63]:
alpha = 0.05 # Уровень значимости
In [64]:
sample = [3.4, 0.5, 0.2, 1, 1.7, 1, 1, 3.9, 4.1, 3.6, 0.5, 0.7,
1.2, 0.5, 0.5, 1.9, 0.5, 0.3, 1.5, 1.9, 1.9, 2.4, 1.2,
2.9,3.2, 1.2, 1.7, 2.9, 1.5, 2.4, 3.4, 0.7, 1.2, 1.7,
1.5, 3.2, 3.9, 1.7, 2.7, 1, 1.5, 1.5, 2.9, 0.7, 2.2, 2.2,
1.9, 1.7, 1.7, 1.9, 1.9, 3.9, 1.2, 1.5, 2.4, 3.3, 2.9,
2.2, 4.6, 3.9, 2.2, 1.2, 3.6, 3.2, 2.2, 2.9, 3.4, 2.4,
2.9, 3.2, 1.7, 1.7, 2.2, 2.7, 3.2, 3.2, 2.9, 1.9, 1.7,
2.2, 1.7, 1.2, 1.2, 1.9, 0.7, 2.2, 1.5, 1.5, 2.7, 4.9,
3.2, 0.7, 2.2, 3.6, 3.6, 1.7, 3.2, 3.4, 1, 0.5, 3.4, 5.3,
4.4, 6.8, 4.6, 3.4, 2.2, 2.2, 2.7, 2.2, 1.2, 1.7, 1.9,
1.2, 1.2, 3.6, 2.4, 1, 2.9, 3.6, 1.7, 2.8] # Выборка из второго практикума
k = 2.00307 # Параметры из второго практикума
l = 2.53379
In [65]:
sps.kstest(sample, sps.weibull_min(c=k, scale=l).cdf)
Out[65]:
Вывод. Выборка согласуется с распределением Вейбулла (так как pvalue больше 0.05, гипотеза не отвергается).
Пусть $X_1, ..., X_n$ --- выборка из распределения $\mathcal{N}(\theta, 1)$. Известно, что $\overline{X}$ является асимптотически нормальной оценкой параметра $\theta$. Вам нужно убедиться в этом, сгенерировав множество выборок и посчитав по каждой из них оценку параметра в зависимости от размера выборки.
Сгенерируйте 200 выборок $X_1^j, ..., X_{300}^j$ из распределения $\mathcal{N}(0, 1)$. По каждой из них посчитайте оценки $\widehat{\theta}_{jn} = \frac{1}{n}\sum\limits_{i=1}^n X_i^j$ для $1 \leqslant n \leqslant 300$, то есть оценка параметра по первым $n$ наблюдениям $j$-й выборки. Для этой оценки посчитайте статистику $T_{jn} = \sqrt{n} \left( \widehat{\theta}_{jn} - \theta \right)$, где $\theta = 0$.
In [66]:
# Код из задачи 3а прака 2
sample = sps.norm.rvs(size=(200, 300))
estimator = sample.cumsum(axis=1) / (np.ones(200).reshape(200, 1) @ np.linspace(1, 300, 300).reshape(1, 300))
stat = estimator * np.linspace(1, 300, 300) ** 0.5
sample = stat[:, -1]
In [67]:
# Проверка критерием Колмогорова
sps.kstest(sample, 'norm')
Out[67]:
Вывод. Выборка согласуется со стандартным нормальным распределением (так как pvalue больше 0.05, гипотеза не отвергается).
Пусть $X_1, ..., X_n$ --- выборка из распределения $Pois(\theta)$. Известно, что $\overline{X}$ является асимптотически нормальной оценкой параметра $\theta$.
In [68]:
# Код из задачи 3б прака 2
sample = sps.poisson.rvs(mu=1, size=(200, 300))
estimator = sample.cumsum(axis=1) / (np.ones(200).reshape(200, 1) @ np.linspace(1, 300, 300).reshape(1, 300))
stat = (estimator - 1) * np.linspace(1, 300, 300) ** 0.5
sample = stat[:, -1]
In [69]:
# Проверка критерием Колмогорова
sps.kstest(sample, 'norm')
Out[69]:
Вывод. Выборка согласуется со стандартным нормальным распределением (так как pvalue больше 0.05, гипотеза не отвергается).
Пусть $X_1, ..., X_n$ --- выборка из распределения $U[0, \theta]$. Из домашнего задания известно, что $n\left(\theta - X_{(n)}\right) \stackrel{d_\theta}{\longrightarrow} Exp\left(1/\theta\right)$.
In [70]:
# Код из задачи 4 прака 2
sample = sps.uniform.rvs(size=(200, 300))
estimator = np.maximum.accumulate(sample, axis=1)
stat = (1 - estimator) * np.linspace(1, 300, 300)
sample = stat[:, -1]
In [71]:
# Проверка критерием Колмогорова
sps.kstest(sample, 'expon')
Out[71]:
Вывод. Выборка согласуется со экспоненциальным распределением с параметром 1 (так как pvalue больше 0.05, гипотеза не отвергается).
Мороженое.
In [72]:
# Я загурзил файл локально
data = np.array(list(csv.reader(open('ice_cream.txt', 'r'), delimiter='\t')))
source = data[1:, :].astype(float)
In [73]:
# Переводим в Цельсий
source[:, 4] = (source[:, 4] - 32) / 1.8
source[:, 5] = (source[:, 5] - 32) / 1.8
In [74]:
# Размер выборки и число параметров
n, k = source.shape[0], 2
In [75]:
# Отклик и регрессор
Y = source[:, 1]
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4]
In [76]:
# Обучаем модель
model = LinearRegression()
model.fit(X, Y)
Out[76]:
In [77]:
# Получаем остатки
errors = Y - model.predict(X)
In [78]:
# Проверяем критерием Колмогороова.
sps.kstest(errors, sps.norm(scale=model.sigma_sq ** 0.5).cdf)
Out[78]:
Вывод. Распределение согласуется с нормальным.
Выбирал всех своих друзей в ВК, к которым удалось получить доступ. Это подсчитывалось в C# через HTTP-запросы. Зато прошедшее время у меня появились новые друзья, я их учитывать не буду. У моих друзей тоже появились новые друзья, они тоже не посчитаны.
In [79]:
from statsmodels.distributions.empirical_distribution import ECDF
In [80]:
source = np.array(list(csv.reader(open('friends.csv', 'r'))))
sample = np.sort(source[:, 1].astype(int))
In [81]:
grid = np.linspace(0, sample.max() + 1, 1000)
В первой практике я говорил, что не могу сделать вывод о распределении на основании такого графика. Может быть, сейчас получится. Выдвинем гипотезу, что распределение числа друзей моих друзей есть распределение Рэлея с параметром масштаба 180. Почему нет, имеем право.
In [82]:
sps.kstest(sample, sps.rayleigh(scale=180).cdf)
Out[82]:
In [83]:
plt.figure(figsize=(20, 5))
plt.plot(grid, ECDF(sample)(grid), color='red', label='ЭФР')
plt.plot(grid, sps.rayleigh(scale=180).cdf(grid), color='blue', label='ФР')
plt.legend()
plt.title('Число друзей моих друзей в ВК')
plt.grid()
plt.show()
In [84]:
plt.figure(figsize=(20, 5))
plt.plot(grid, sps.rayleigh(scale=180).pdf(grid), color='blue', label='Плотность')
plt.legend()
plt.title('Плотность распределения Рэлея')
plt.grid()
plt.show()
Вывод. Не отвергается гипотеза, что распределения числа друзей моих друзей имеет распределение Рэлея с параметром масштаба 180. График функции распределения примерно совпадает с графиком ЭФР для выборки. Это очень здорово. С помощью критерия Колмогорова получилось найти неплохо приближающее распределение. Однако если почитать о том, где данное распреление встречается, оказывается, что оно находит применение в физике, так что встретить его в такой стезе весьма оказалось весьма неожиданно.
Задача 8\*. Проведите исследование согласно примеру 2 параграфа 2 главы 18 книги М.Б. Лагутина "Наглядная математическая статистика".
Задача 9\*. Изучите Q-Q plot и критерий Шапиро-Уилка для проверки нормальности, напишите их теоретическое пояснение. В изучении могут помочь материалы курса ПСАД.
Постройте графики Q-Q plot для различных распределений и дайте к ним пояснение. Проверьте различные данные на нормальность с помощью различных критериев и Q-Q plot. Данные можно использовать из задачи 7 или какие-либо еще, например, отдельные компоненты из Ирисов Фишера. Постарайтесь так же правильно контролировать вероятность общей ошибки первого рода.