#!/usr/bin/env python
# coding: utf-8

# # Задача 1
# ## Условие
# Загрузите данные из набора Forest Fires (файл forestfires.csv) о лесных пожарах в Португалии. Задача состоит в том, чтобы с помощью линейной регрессии научиться предсказывать координату area (площадь пожара) в виде линейной комбинации других данных.
#
# Чтобы работать с числовыми координатами, нечисловые координаты (month, day) нужно перевести в числовые. Для простоты можно заменить координату month на индикатор летнего сезона, а координату day не использовать вообще. По желанию можно сделать преобразование другим способом. Также добавьте координату, тождественно равную единице (вес при этой координате  интерпретируется как сдвиг).
#
# Разбейте выборку на две части в соотношении 7:3 (перемешав её с помощью random.shuffle). По первой части постройте регрессионную модель. Примените модель ко второй части выборки и посчитайте по ней среднеквадратичную ошибку.
#
# Для переменной area выполните преобразование $f(x) = ln(x+c)$ и постройте для нее новую регрессионную модель. Посчитайте среднеквадратичную ошибку для преобразованных значений. При каком $c$ предсказания получаются лучше всего?
#
# При выбранном $c$ сделайте разбиение выборки в соотношении 7:3 разными способами (перемешивая каждый раз). Найдите способ оценить разброс качества от разбиения. Сильно ли меняется качество? Сделайте выводы.

# In[7]:


import scipy.stats as sps
import matplotlib.pyplot as plt
import numpy as np
from sortedcontainers import SortedList
from IPython.display import display, Markdown
import pandas as pd
params = {'legend.fontsize': 15,
          'legend.handlelength': 2,
           'figure.figsize': (15, 10)}
plt.rcParams.update(params)


# In[71]:


df = pd.read_csv('./forestfires.csv')


# In[72]:


df = df[df['area'] < df['area'].quantile(0.9)]
summer_season = ['jun', 'jul', 'aug']
conditions = [(df['month'] == 'jun') | (df['month'] == 'jul') | (df['month'] == 'aug')]
choices = [1]
df['month'] = np.select(conditions, choices)
df['shift'] = 1
del df['day']


# In[73]:


df = df.sample(frac=1)
quant = int(.7*len(df))
training_df, testing_df = df[:quant], df[quant:]


# In[74]:


def linear_regression(X, Z):
    Zt = np.transpose(Z)
    theta = np.dot(np.dot(np.linalg.inv(np.dot(Zt, Z)), Zt), X)
    est = X - np.dot(Z, theta)
    sigma = np.dot(est, est)/(len(Z) - len(Z[0]))
    return theta, sigma

def getXZ(df):
    return df['area'].to_numpy(), df.loc[:, df.columns != 'area'].to_numpy()

def compute_estimation(training_df):
    theta, sigma = linear_regression(*getXZ(training_df))
    return theta, sigma

def apply_regression(training_df, testing_df):
    theta, sigma = compute_estimation(training_df)
    X, Z = getXZ(testing_df)
    mse = ((np.dot(Z, theta) - X)**2).mean()
    return mse


# In[75]:


x_axis = np.arange(1, len(testing_df) + 1)
plt.plot(x_axis, testing_df['area'], label='Area')
_, Z = getXZ(testing_df)
theta, _= linear_regression(*getXZ(training_df))
plt.plot(x_axis, np.dot(Z, theta), label='Estimation of area')
plt.xlabel('Number')
plt.ylabel('Area')
plt.title('Линейная регрессия для площади')
plt.legend(loc='best')
plt.show()
mse = apply_regression(training_df, testing_df)
print(mse)


# In[80]:


def get_log_mse(training_df, testing_df, c):
    _, Z_testing = getXZ(testing_df)
    _, Z_training = getXZ(training_df)
    testing_df_log_area =  testing_df['area'].apply(lambda x: np.log(x + c))
    training_df_log_area =  training_df['area'].apply(lambda x: np.log(x + c))
    theta, sigma = linear_regression(training_df_log_area, Z_training)
    X, _ = getXZ(testing_df)
    mse = ((np.exp(np.dot(Z_testing, theta)) - c - X)**2).mean()
    return mse


# In[77]:


min_mse = mse
min_c = 0
_, Z_testing = getXZ(testing_df)
_, Z_training = getXZ(training_df)
mse_list = np.array([])
x_axis = np.arange(0, 100, 0.1)
for c in x_axis:
    testing_df_log_area =  testing_df['area'].apply(lambda x: np.log(x + c))
    training_df_log_area =  training_df['area'].apply(lambda x: np.log(x + c))
    theta, sigma = linear_regression(training_df_log_area, Z_training)
    X, _ = getXZ(testing_df)
    mse = ((np.exp(np.dot(Z_testing, theta)) - c - X)**2).mean()
    mse_list = np.append(mse_list, mse)
    if mse < min_mse:
        min_mse = mse
        min_c = c
        theta_min = theta
print(min_mse, min_c)
plt.plot(x_axis, mse_list, label='MSE')
plt.scatter(min_c, min_mse, color='red', label='minimum')
plt.xlabel('c')
plt.ylabel('MSE')
plt.legend(loc='best')
plt.title('Зависимость MSE от c')
plt.show()


# In[ ]:





# In[84]:


x_axis = np.arange(1, len(testing_df) + 1)
plt.plot(x_axis, testing_df['area'], label='Area')
plt.plot(x_axis, np.exp(np.dot(Z_testing, theta_min)) - min_c, label='Estimation')
plt.xlabel('Number')
plt.ylabel('Area')
plt.title('Линейная регрессия для площади. Экспонециальная зависимость')
plt.legend(loc='best')
plt.show()


# In[82]:


for tries in range(10):
    df = df.sample(frac=1)
    quant = int(.7*len(df))
    training_df, testing_df = df[:quant], df[quant:]
    print(get_log_mse(training_df, testing_df, min_c))


# # Вывод
# Мы получили оценку на величину, сначала в предположении ее линейной зависимости от остальных, а потом в предположении сдвинутой экспонециальной зависимости. Я убрал из рассмотрения очень сильно отклоняющиеся значения приближаемого поля(оставил квантиль 0.9), чтобы общий тренд увидеть более четче. Мы видим, что mse для второго предположения получается примерно такой же, и ошибка для оптимального c сильно зависит от разбиения. Из этого можно сделать вывод, что обе гипотезы не очень правдоподобны, то есть, вероятно, area не выражается подобным способом.

# # Задача 2
# ## Условие
# Пусть $$X_i = \beta_1+i\beta_2+\varepsilon_0+\ldots+\varepsilon_i,\, i= 0, 1, \ldots, n$$ --- расстояния, которое проехал трамвай за $i$ секунд по показанию датчика. Здесь $\beta_1$~--- начальное расстояние, $\beta_2$ --- скорость трамвая, $\varepsilon_0$ --- ошибка начального показания датчика. Трамвай едет с постоянной скоростью, и через каждую секунду датчик фиксирует расстояние, которое проехал трамвай. Отсчет времени идет от предыдущего замера, причем отсчет происходит с ошибкой. Для $i = 1, \ldots, n$ величина $\varepsilon_i$ есть ошибка приращения расстояния, то есть $\varepsilon_i = \varepsilon_i^t \beta_2,$ где $\varepsilon_i^t$ --- ошибка отсчета времени. Все ошибки $\varepsilon_i$ независимы и распределены по закону $N(0, \sigma^2)$.
#
# Сведите задачу к линейной модели и найдите оценки наименьших квадратов для начального расстояния $\beta_1$ и скорости $\beta_2,$ а также несмещенную оценку для $\sigma^2,$ из которой выразите оценку дисперсии отсчета времени.
#
# Данные возьмите из файла Regression.csv. Сделайте выводы.
#
# ## Решение
# В домашнем задании мы показали, что при сведении к линейной регрессии($Y_0 = X_0$, $Y_i = X_i - X_{i-1}$), получаются оценки
# * $\hat{\beta_1} = Y_0$
# * $\hat{\beta_2} = \frac{\sum\limits_{i = 1}^{n}{Y_i}}{n}$
# * $\hat{\sigma^2} = \frac{n}{n-1}s^2(Y_1, \ldots, Y_n)$
#
# Т.к. $\varepsilon_i \sim N(0, \sigma^2), \varepsilon_i^t = \frac{\varepsilon_i}{\beta_2}$, то
# $\varepsilon_i^t \sim N(0, \frac{\sigma^2}{\beta_2^2})$.
# Возьмем тогда $\hat{\varepsilon_i^t} = \frac{\hat{\sigma^2}}{\hat{\beta_2^2}}$

# In[89]:


df = pd.read_csv('./regression.csv')


# In[92]:


def get_estimations(Xs):
    Ys = np.zeros(len(Xs))
    Ys[0] = Xs[0]
    for i in range(1, len(Xs)):
        Ys[i] = Xs[i] - Xs[i-1]
    beta_1 = Ys[0]
    beta_2 = Ys[1:].mean()
    n = len(Xs) - 1
    sigma = n/(n - 1)*np.var(Ys[1:])
    time_error = sigma / beta_2 ** 2
    return beta_1, beta_2, sigma, time_error


# In[120]:


Xs = df.to_numpy()
beta_1, beta_2, sigma, time_error = get_estimations(Xs)
print('beta_1: ', beta_1)
print('beta_2: ', beta_2)
print('sigma: ', sigma)
print('time_error: ', time_error)


# In[123]:


estims = np.array([])
Xs = Xs.reshape(1, len(Xs))[0]
for i in range(3, len(Xs) + 1):
    estims = np.append(estims, get_estimations(Xs[:i])[3])
plt.xlabel('Номер')
plt.ylabel('Дисперсия ошибки отсчета времени')
plt.title('Дисперсия ошибки отсчета времени')
plt.plot(np.arange(1, len(estims) + 1), estims, label='Дисперсия ошибки отсчета времени')
plt.legend(loc='best')
plt.show()


# # Вывод
# Мы получили, что