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

%matplotlib inline

"""
Explanation: Математическая статистика
Практическое задание 5
В данном задании предлагается провести некоторое исследование модели линейной регрессии и критериев для проверки статистических гипотез, в частности применить этим модели к реальным данным.
Правила:

Выполненную работу нужно отправить на почту probability.diht@yandex.ru, указав тему письма "[номер группы] Фамилия Имя - Задание 5". Квадратные скобки обязательны. Вместо Фамилия Имя нужно подставить свои фамилию и имя.
Прислать нужно ноутбук и его pdf-версию. Названия файлов должны быть такими: 5.N.ipynb и 5.N.pdf, где N - ваш номер из таблицы с оценками.
Никакой код из данного задания при проверке запускаться не будет.
Некоторые задачи отмечены символом <font size="5" color="red">*</font>. Эти задачи являются дополнительными. Успешное выполнение большей части таких задач (за все задания) является необходимым условием получения бонусного балла за практическую часть курса.
Баллы за каждую задачу указаны далее. Если сумма баллов за задание меньше 25% (без учета доп. задач), то все задание оценивается в 0 баллов.

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

Задача 1 - 7 баллов
Задача 2 - 2 балла
Задача 3<font size="5" color="red">*</font> - 3 балла
Задача 4 - 2 балла
Задача 5<font size="5" color="red">*</font> - 10 баллов
Задача 6 - 5 баллов
Задача 7 - 4 балла
Задача 8<font size="5" color="red">*</font> - 4 балла
Задача 9<font size="5" color="red">*</font> - 10 баллов
End of explanation
"""


from scipy.linalg import inv
from numpy.linalg import norm

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

"""
Explanation: 1. Линейная регрессия
Задача 1. По шаблону напишите класс, реализующий линейную регрессию. Интерфейс этого класса в некоторой степени соответствует классу <a href="http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html#sklearn.linear_model.LinearRegression">LinearRegression</a> из библиотеки sklearn.
End of explanation
"""


import csv

# Я загурзил файл локально
data = np.array(list(csv.reader(open('ice_cream.txt', 'r'), delimiter='\t')))
source = data[1:, :].astype(float)
print(source[:5])

# Переводим в Цельсий
source[:, 4] = (source[:, 4] - 32) / 1.8
source[:, 5] = (source[:, 5] - 32) / 1.8

# Размер выборки и число параметров
n, k = source.shape[0], 2
print(n)

# Отклик и регрессор
Y = source[:, 1]
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4]
print(X[:5])

# Обучаем модель
model = LinearRegression()
model.fit(X, Y)

# Выводим общую информацию 
model.summary()

grid = np.linspace(-5, 25, 1000)

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()

"""
Explanation: Загрузите данные о потреблении мороженного в зависимости от температуры воздуха и цены (файл 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$ --- МНК-оценки коэффициентов.
На график нанесите точки выборки.
Убедитесь, что построейнный график совпадает с графиком из презентации с первой лекции, правда, с точностью до значений температура (она была неправильно переведена из Фаренгейта в Цельсий).
End of explanation
"""


# Размер выборки и число параметров
n, k = source.shape[0], 4

# Отклик и регрессор
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])

# Обучаем модель
model = LinearRegression()
model.fit(X, Y)

# Выводим общую информацию 
model.summary()

grid = np.linspace(-5, 25, 1000)

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]

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()

"""
Explanation: Вывод. Действительно, график тот же (с поправкой на пересчет температуры). Линейная регрессия неплохо, но не идеально приближает зависимость потребления мороженого в зависимости от температуры. Для более точного вывода стоит посчитать доверительный интервал.
Теперь учтите влияние года (столбец Year) для двух случаев:
* модель $ic = \theta_1 + \theta_2\ t + \theta_3 y_1 + \theta_4 y_2$, где $y_1 = I{1\ год}, y_2 = I{2\ год}$. Поясните, почему нельзя рассмативать одну переменную $y$ --- номер года.
* для каждого года рассматривается своя линейная зависимость $ic = \theta_1 + \theta_2\ t$.
В каждом случае нарисуйте графики. Отличаются ли полученные результаты? От чего это зависит? Как зависит потребление мороженного от года? 
Первый случай.
End of explanation
"""


# Число признаков
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)

# Три подвыборки
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])

# Отклики и регрессоры
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])

# Обучаем модели
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)

# Выводим общую информацию 
model_0.summary()
model_1.summary()
model_2.summary()

grid = np.linspace(-5, 25, 1000)

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()

"""
Explanation: Вывод. В разрезе разных лет также можно уследить, что линейная зависимсоть неплохо приближает реальную. Номер года нельзя брать как признак, так это больше характеризация класса каких-то значений, в то время линейная регрессия работает со значениями, от которых зависимость линейная, численная
Второй случай.
End of explanation
"""


# Размер выборки и число параметров
n, k = source.shape[0], 8
print(n, k)

# 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]

# Обучаем модель
model = LinearRegression()
model.fit(X, Y)

# Выводим общую информацию 
model.summary()

"""
Explanation: Вывод. При разделении выборки на три части результаты в какой-то степени стали лучше. Кроме того, такое разделение позволяет заметить, что с при увеличении года в целом потребляют больше мороженого, и что прирост литров мороженого на один градус возрастает при увеличении года.
Наконец, обучите модель на предсказание потребления мороженного в зависимости от всех переменных.
Не забудьте, что для года нужно ввести две переменных.
Для полученной модели выведите summary.
End of explanation
"""


# Размер выборки и число параметров
n, k = source.shape[0], 4
print(n, k)

# Отклик и регрессор
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])

# Обучаем модель
model = LinearRegression()
model.fit(X, Y)

# Выводим общую информацию 
model.summary()

grid = np.linspace(-5, 25, 1000)

y = model.theta[0] + model.theta[1] * grid + model.theta[2] * grid ** 2 + model.theta[3] * grid ** 3

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()

"""
Explanation: Вывод. Похоже, не все признаки стоит учитывать. Некоторые значения по модулю очень маленькие по сравнению со значениями, характерными для соответствующих признаков, что их лучше не учитывать. Например, это номер измерения (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$. Хорошие ли получаются результаты?
End of explanation
"""


D = inv(X.T @ X)
print(D)

from scipy.linalg import eigvals

vals = eigvals(D)
print(vals)   # Комплексных нет, мы не налажали
vals = vals.real.astype(float)
print(vals)

CI = np.sqrt(vals.max() / vals.min())
print(CI)

"""
Explanation: Вывод. Результаты выглядят более естественно, однако теперь кажется, что при 30-40 градусах люди потребляют невероятно много мороженого, что, скорее всего, неправда (все-таки есть какой-то лимит). Кроме того, значения последних двух параметров очень малы, что говорит о малом влиянии этих параметров при малых температурах. Возможно, линейная модель и так достаточно хороша. Однако в нашем случае будет видно, что при отрицательных температурах потребление мороженого падает стремительно. Наверное, в этом есть некоторая правда. 
Чтобы понять, почему так происходит, выведите значения матрицы $(X^T X)^{-1}$ для данной матрицы и посчитайте для нее индекс обусловленности $\sqrt{\left.\lambda_{max}\right/\lambda_{min}}$, где $\lambda_{max}, \lambda_{min}$ --- максимальный и минимальный собственные значения матрицы $X^T X$. Собственные значения можно посчитать функцией <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigvals.html">scipy.linalg.eigvals</a>.
Прокомментируйте полученные результаты. Помочь в этом может следующая <a href="https://ru.wikipedia.org/wiki/%D0%A7%D0%B8%D1%81%D0%BB%D0%BE_%D0%BE%D0%B1%D1%83%D1%81%D0%BB%D0%BE%D0%B2%D0%BB%D0%B5%D0%BD%D0%BD%D0%BE%D1%81%D1%82%D0%B8">статья</a>.
End of explanation
"""


from sklearn import linear_model
from sklearn import cross_validation
from sklearn.metrics import mean_squared_error

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)

"""
Explanation: Вывод. Как говорилось на семинаре, высокий индекс обусловленности (больше 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.
Для выполнения задания воспользуйтесь следующими функциями:
* <a href="http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html#sklearn.linear_model.LinearRegression">sklearn.linear_model.LinearRegression</a>
--- реализация линейной регрессии. В данной реализации свободный параметр $\theta_1$ по умолчанию автоматически включается в модель. Отключить это можно с помощью fit_intercept=False, но это не нужно. В данной задаче требуется, чтобы вы воспользовались готовой реализацией линейной регрессии, а не своей. Ведь на практике важно уметь применять готовые реализации, а не писать их самостоятельно.


<a href="http://scikit-learn.org/0.16/modules/generated/sklearn.cross_validation.train_test_split.html">sklearn.cross_validation.train_test_split</a>
--- функция разбиения данных на train и test. Установите параметр test_size=0.3.


<a href="http://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_error.html">sklearn.metrics.mean_squared_error</a>
--- реализация MSE.


Для перебора реализуйте функцию.
End of explanation
"""


# Я загурзил файл локально и исправил проблелы на табы (там были лишние проблеы в некоторых местах)
data = np.array(list(csv.reader(open('yh.data', 'r'), delimiter='\t')))
yacht = data.astype(float)
print(yacht)

Y = yacht[:, 6]
X = yacht[:, :6]

X_train, X_test, Y_train, Y_test = cross_validation.train_test_split(X, Y, test_size=0.3)

best_features(X_train, X_test, Y_train, Y_test)

"""
Explanation: Примените реализованный отбор признаков к датасетам
* <a href="http://archive.ics.uci.edu/ml/datasets/Yacht+Hydrodynamics">Yacht Hydrodynamics</a> --- для парусных яхт нужно оценить остаточное сопротивление на единицу массы смещения (последний столбец) в зависимости от различных характеристик яхты.

<a href="http://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_boston.html#sklearn.datasets.load_boston">Boston Housing Prices</a> --- цены на дома в Бостоне в зависимости от ряда особенностей.

Yacht Hydrodynamics.
End of explanation
"""


# Я загурзил файл локально и убрал первую стрчоку с кучей запятых
data = np.array(list(csv.reader(open('bhp.csv', 'r'))))
houses = data[1:, :].astype(float)
print(houses)

# Число столбцов
cols = houses.shape[1]

Y = houses[:, 9]  # Столбец TAX
X = houses[:, np.delete(np.arange(cols), 9)]

X_train, X_test, Y_train, Y_test = cross_validation.train_test_split(X, Y, test_size=0.3)

best_features(X_train, X_test, Y_train, Y_test)

"""
Explanation: Вывод. Признак 5 (Froude number) наиболее полезный признак, остальные встречаются как-то рандомно (кроме 4, который почти не встречается, это Length-beam ratio).
Boston Housing Prices.
End of explanation
"""


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))

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()

"""
Explanation: Вывод. Первые 3 и последние 3 признака самые полезные.
Задача 3<font size="5" color="red">*</font>. Загрузите <a href="http://people.sc.fsu.edu/~jburkardt/datasets/regression/x01.txt">датасет</a>, в котором показана зависимость веса мозга от веса туловища для некоторых видов млекопитающих. Задача состоит в том, чтобы подобрать по этим данным хорошую модель регрессии. Для этого, можно попробовать взять некоторые функции от значения веса туловища, например, степенную, показательную, логарифмическую. Можно также сделать преобразование значений веса мозга, например, прологарифмировать. Кроме того, можно разбить значения веса туловища на несколько частей и на каждой части строить свою модель линейной регрессии.
Задача 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)$ .
End of explanation
"""


alpha = 0.05  # Уровень значимости

# Считаем параметры для выборок четырех разных размеров
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'])

"""
Explanation: Вывод. При увеличении размера выборки площадь доверительной область (она, очевидно, измерима даже по Жордану) сильно уменьшается. Это говорит о том, что чем больше значений получено, тем точнее с заданным уровнем доверия можно оценить параметр.
Задача 5<font size="5" color="red">*</font>. 
Пусть дана линейная гауссовская модель $Y = X\theta + \varepsilon$, где $\varepsilon \sim \mathcal{N}(0, \beta^{-1}I_n)$.
Пусть $\theta$ имеет априорное распределение $\mathcal{N}(0, \alpha^{-1}I_k)$. 
Такая постановка задачи соответствует Ridge-регрессии.
Оценкой параметров будет математическое ожидание по апостериорному распределению, аналогично можно получить доверительный интервал.
Кроме того, с помощью апостериорного распределения можно получить доверительный интервал для отклика на новом объекте, а не только точечную оценку.
Реализуйте класс RidgeRegression подобно классу LinearRegression, но добавьте в него так же возможность получения доверительного интервала для отклика на новом объекте.
Примените модель к некоторых датасетам, которые рассматривались в предыдущих задачах.
Нарисуйте графики оценки отклика на новом объекте и доверительные интервалы для него.
2. Проверка статистических гипотез
Задача 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. Внимательно проверьте правильность строгих и нестрогих знаков.
End of explanation
"""


# В двух случаях строим для 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']))

"""
Explanation: Вывод. Если $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, если отвергается).
Какие выводы можно сделать?
End of explanation
"""


p_conf = 0.75  # Если больше, то уверены в том, что действительно неудачи связаны с черной кошкой.
beta = 0.8  # Магическая константа с семинара и из условия

grid = np.linspace(0.501, 0.999, 500)

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()

"""
Explanation: Вывод. Выходит, почти всегда на малых выборках даже при большом $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$.
End of explanation
"""


size_conf = 27  # Оптимальные размер выборки

# В двух случаях строим для 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']))

"""
Explanation: Вывод. Оптимальный в некотором смысле размер выборки получается примерно 25–30, на пересечении $\theta=p^*=0.75$ и $\beta(\theta, X)=\beta=0.8$. Выберем 27.
End of explanation
"""


alpha = 0.05  # Уровень значимости

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

sps.kstest(sample, sps.weibull_min(c=k, scale=l).cdf)

"""
Explanation: Вывод. При выбранном значении $p^$ и подобранном оптималном значении выборки если $p < p^$, то гипотеза не отвергается почти всегда, а при $p > p^*$ гипотеза отвергается почти всегда.
Справка для выполнения следующих задач
Критерий согласия хи-квадрат
<a href=https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html#scipy.stats.chisquare>scipy.stats.chisquare</a>(f_obs, f_exp=None, ddof=0)
f_obs --- число элементов выборки, попавших в каждый из интервалов
f_exp --- ожидаемое число элементов выборки (по умолчанию равномерное)
ddof --- поправка на число степеней свободы. Статистика асимптотически будет иметь распределение хи-квадрат с числом степеней свободы $k - 1 - ddof$, где $k$ --- число интервалов.
Возвращает значение статистики критерия и соответствующее p-value.
Критерий согласия Колмогорова
<a href=https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kstest.html#scipy.stats.kstest>scipy.stats.kstest</a>(rvs, cdf, args=())
rvs --- выборка
cdf --- функция распределения (сама функция или ее название)
args --- параметры распределения
Возвращает значение статистики критерия и соответствующее p-value.
Задача 7.


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


Проверьте, что при больших $n$ распределение статистики из задач 3 и 4 задания 2 действительно хорошо приближают предельное распределение.


Проверьте, что остатки в регрессии из задач выше нормальны.


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


Использовать можно два описанных выше критерия, либо любой другой критерий, если будет обоснована необходимость его применения в данной задаче, а так же будет приведено краткое описание критерия.
Уровень значимости взять равным 0.05.
End of explanation
"""


# Код из задачи 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]

# Проверка критерием Колмогорова
sps.kstest(sample, 'norm')

"""
Explanation: Вывод. Выборка согласуется с распределением Вейбулла (так как 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$.
End of explanation
"""


# Код из задачи 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]

# Проверка критерием Колмогорова
sps.kstest(sample, 'norm')

"""
Explanation: Вывод. Выборка согласуется со стандартным нормальным распределением (так как pvalue больше 0.05, гипотеза не отвергается).
Пусть $X_1, ..., X_n$ --- выборка из распределения $Pois(\theta)$. Известно, что $\overline{X}$ является асимптотически нормальной оценкой параметра $\theta$.
End of explanation
"""


# Код из задачи 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]

# Проверка критерием Колмогорова
sps.kstest(sample, 'expon')

"""
Explanation: Вывод. Выборка согласуется со стандартным нормальным распределением (так как 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)$.
End of explanation
"""


# Я загурзил файл локально
data = np.array(list(csv.reader(open('ice_cream.txt', 'r'), delimiter='\t')))
source = data[1:, :].astype(float)

# Переводим в Цельсий
source[:, 4] = (source[:, 4] - 32) / 1.8
source[:, 5] = (source[:, 5] - 32) / 1.8

# Размер выборки и число параметров
n, k = source.shape[0], 2

# Отклик и регрессор
Y = source[:, 1]
X = np.zeros((n, k))
X[:, 0] = np.ones(n)
X[:, 1] = source[:, 4]

# Обучаем модель
model = LinearRegression()
model.fit(X, Y)

# Получаем остатки
errors = Y - model.predict(X)

# Проверяем критерием Колмогороова.
sps.kstest(errors, sps.norm(scale=model.sigma_sq ** 0.5).cdf)

"""
Explanation: Вывод. Выборка согласуется со экспоненциальным распределением с параметром 1 (так как pvalue больше 0.05, гипотеза не отвергается).
Мороженое.
End of explanation
"""


from statsmodels.distributions.empirical_distribution import ECDF

source = np.array(list(csv.reader(open('friends.csv', 'r'))))
sample = np.sort(source[:, 1].astype(int))

grid = np.linspace(0, sample.max() + 1, 1000)

"""
Explanation: Вывод. Распределение согласуется с нормальным.
Выбирал всех своих друзей в ВК, к которым удалось получить доступ. Это подсчитывалось в C# через HTTP-запросы. Зато прошедшее время у меня появились новые друзья, я их учитывать не буду. У моих друзей тоже появились новые друзья, они тоже не посчитаны.
End of explanation
"""


sps.kstest(sample, sps.rayleigh(scale=180).cdf)

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()

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()

"""
Explanation: В первой практике я говорил, что не могу сделать вывод о распределении на основании такого графика. Может быть, сейчас получится. Выдвинем гипотезу, что распределение числа друзей моих друзей есть распределение Рэлея с параметром масштаба 180. Почему нет, имеем право.
End of explanation
"""


