# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Пуассоновская регрессия
# ### <NAME>
# 2021

# +
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
# %matplotlib inline
plt.style.use('ggplot')
# -

# ## 1. Чтение и подготовка данных

# Рассмотрим данные о количестве велосипедистов. Количество велосипедистов зависит от погодных условий в рассматриваемый день: чем хуже погода, тем меньше желающих. В качестве признаков возьмем:
# - максимальную температуру в рассматриваемый день (F);
# - минимальную температуру в рассматриваемый день (F);
# - количество осадков.

data = pd.read_csv('data/nyc_bicyclist_counts.csv', index_col=['Date'], parse_dates=True)

data.head()

# Целевая переменная &ndash; `'BB_COUNT'` &ndash; содержит только целые положительные числа, что должно быть учтено при выборе предсказательной модели.

data['BB_COUNT'].plot(figsize=(12,5))
plt.show()

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

data['DAY_OF_WEEK'] = data.index.dayofweek
data['MONTH'] = data.index.month

data

# Данные переменные являются категориальными.

# <div class="alert alert-info">
#
# <h3> Задание 1</h3>
# <p></p>
# <ul>1. Определите функцию, которая принимает на вход исходные данные $(X,y)$ и параметры модели $\theta$. Данная функция должна возвращать среднеквадратичную ошибку модели.</ul>
#     
# <ul>2. Определите аналогичную функцию, которая возвращает значение функционала качества пуассоновской регрессии.</ul>  
#     
# <ul>3. Обучите обе модели с помощью функции minimize из SciPy. Сравните качество аппроксимации моделей. Метрикой качества выберите среднюю абсолютную ошибку.</ul>    
#  
# <ul>4. Отобразите на графике исходный ряд и результаты аппроксимации линейной и пуассоновской регрессиями.</ul>
# <p></p>   
#
# <p></p>
# </div>

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score, KFold
from sklearn.compose import TransformedTargetRegressor
from scipy.optimize import minimize
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error


def mse(X,y,theta):
    return ((y-np.dot(X,theta))**2).mean()


X=data.drop('BB_COUNT', axis=1)
y = data['BB_COUNT']

X['const'] = 1

theta0 = np.ones(X.shape[1])

lin_reg = minimize(lambda theta: mse(X,y,theta), tuple(theta0))
lin_reg_params = lin_reg.x

y_pred_lin = np.dot(X, lin_reg_params)

mean_absolute_error(y,y_pred_lin)


def pois(X,y,theta):
    mu = np.dot(X,theta)
    return (np.exp(mu)-y*mu).mean()


theta0 = np.zeros(X.shape[1])

pois_reg = minimize(lambda theta1: pois(X,y,theta1), tuple(theta0))
pois_reg_params = pois_reg.x

y_pred_pois = np.dot(X, pois_reg_params)
data['pois_approx'] = y_pred_pois

mean_absolute_error(y,data['pois_approx'])

data['lin_approx'] = y_pred_lin
data['pois_approx'] = np.exp(y_pred_pois)

a = data['BB_COUNT'].plot(figsize=(15,7), label = 'initial')
b = data['lin_approx'].plot(label = 'lin')
c = data['pois_approx'].plot(label = 'pois')
a.legend()
b.legend()
c.legend()
plt.show()

# <div class="alert alert-info">
#
# <h3> Задание 2</h3>
# <p></p>
# <ul>Линейные модели чувствительны к виду категориальных признаков. Преобразуйте категориальные признаки с помощью One Hot Encoding и повторите шаги 3-4 из задания 1. Как изменилось качество моделей?</ul>
# <p></p>
# </div>

data = pd.read_csv('data/nyc_bicyclist_counts.csv', index_col=['Date'], parse_dates=True)

data['DAY_OF_WEEK'] = data.index.dayofweek
data['MONTH'] = data.index.month

X=data.drop(['BB_COUNT','DAY_OF_WEEK','MONTH'], axis = 1)
y = data['BB_COUNT']

from sklearn.preprocessing import OneHotEncoder

enc = OneHotEncoder(handle_unknown='ignore')

enc.fit(data[['DAY_OF_WEEK','MONTH']])

X[['SUN','MON','TUES','WEN','THUR','FRI','SAT','APR','JUN','JUL','AUG','SEP','OCT','NOV']] = enc.transform(data[['DAY_OF_WEEK','MONTH']]).toarray()

X['const'] = 1

theta0 = np.ones(X.shape[1])

lin_reg = minimize(lambda theta: mse(X,y,theta), tuple(theta0))
lin_reg_params = lin_reg.x

y_pred_lin = np.dot(X, lin_reg_params)

mean_absolute_error(y,y_pred_lin)

theta0 = np.zeros(X.shape[1])

pois_reg = minimize(lambda theta1: pois(X,y,theta1), tuple(theta0))
pois_reg_params = pois_reg.x

y_pred_pois = np.dot(X, pois_reg_params)
data['pois_approx'] = y_pred_pois

mean_absolute_error(y,data['pois_approx'])

data['lin_approx'] = y_pred_lin
data['pois_approx'] = np.exp(y_pred_pois)

a = data['BB_COUNT'].plot(figsize=(15,7), label = 'initial')
b = data['lin_approx'].plot(label = 'lin')
c = data['pois_approx'].plot(label = 'pois')
a.legend()
b.legend()
c.legend()
plt.show()

# <div class="alert alert-info">
#
# <h3> Задание 3</h3>
# <p></p>
# <ul>Преобразуйте категориальные признаки с помощью Фурье-разложения и повторите шаги 3-4 из задания 1. Какого качества моделей удалось достичь?</ul>
# <p></p>
# </div>

data = pd.read_csv('data/nyc_bicyclist_counts.csv', index_col=['Date'], parse_dates=True)

data['DAY_OF_WEEK'] = data.index.dayofweek
data['MONTH'] = data.index.month

X=data.drop(['BB_COUNT','DAY_OF_WEEK','MONTH'], axis = 1)
y = data['BB_COUNT']

X[['SIN_M','COS_M','SIN_W','COS_W']] = np.array([np.sin(2*np.pi/7*data['MONTH']),np.cos(2*np.pi/7*data['MONTH']),np.sin(2*np.pi/12*data['DAY_OF_WEEK']),np.cos(2*np.pi/12*data['DAY_OF_WEEK'])]).transpose()

X['const'] = 1

theta0 = np.zeros(X.shape[1])

lin_reg = minimize(lambda theta: mse(X,y,theta), tuple(theta0))
lin_reg_params = lin_reg.x

lin_reg_params

y_pred_lin = np.dot(X, lin_reg_params)

mean_absolute_error(y,y_pred_lin)

pois_reg = minimize(lambda theta1: pois(X,y,theta1), tuple(theta0))
pois_reg_params = pois_reg.x

y_pred_pois = np.dot(X, pois_reg_params)

data['lin_approx'] = y_pred_lin
data['pois_approx'] = np.exp(y_pred_pois)

mean_absolute_error(y,data['pois_approx'])

a = data['BB_COUNT'].plot(figsize=(15,7), label = 'initial')
b = data['lin_approx'].plot(label = 'lin')
c = data['pois_approx'].plot(label = 'pois')
a.legend()
b.legend()
c.legend()
plt.show()



StarcoderdataJupyterScriptsDedupFiltered