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

# +
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.formula.api import ols


from scipy import stats
from sklearn.metrics import mean_squared_error, r2_score, explained_variance_score

from sklearn.linear_model import LinearRegression
from sklearn.feature_selection import f_regression 
from math import sqrt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
import evaluate
import math
# -

# ### 1. Load the tips dataset from either pydataset or seaborn.

# +
from pydataset import data

tips = data('tips')
tips.head()
# -

# ### 2. Fit a linear regression model (ordinary least squares) and compute yhat, predictions of tip using total_bill.
#
#

# +
# create the model object
lm = LinearRegression(normalize=True)

# fit the model to trainig data
lm.fit(tips[['total_bill']], tips.tip)

# make prediction
tips['yhat'] = lm.predict(tips[['total_bill']])
# -

# make a baseline prediction (mean of the tip)
tips['baseline'] = tips.tip.mean()

tips.head()

# +
# plot data and prediction line

sns.scatterplot(x = 'total_bill', y = 'tip', data = tips)
sns.lineplot(x = 'total_bill', y = 'baseline', data = tips)
sns.lineplot(x = 'total_bill', y = 'yhat', data = tips)
# -

# ### 3. Plot the residuals for the linear regression model that you made.
#
#

tips['residual'] = tips.tip - tips.yhat
tips['baseline_residual'] = tips.tip - tips.baseline

tips.head()

sns.relplot(x = 'total_bill', y = 'residual', data = tips)
plt.axhline(0, ls = ':')

# Heteroscedasticity: 
#   - unequal variance of errors
#   - Heteroscedasticity may also have the effect of giving too much weight to a small subset of the data (namely the subset where the error variance was largest) when estimating coefficients.
#   - possibily apply some transformations(?)
#

# ### 4. Calculate the sum of squared errors, explained sum of squares, total sum of squares, mean squared error, and root mean squared error for your model.
#
#

tips.head()

SSE = (tips.residual**2).sum()
print(f' The SSE of the OLS model is {round(SSE,1)}')

SSE_baseline = (tips.baseline_residual**2).sum()
print(f' The SSE of the baseline model is {round(SSE_baseline,1)}')

# +
#Mean squared error OLS:

MSE = SSE/len(tips)
MSE

# +
#Mean squared error of baseline model:

MSE_baseline = SSE_baseline/len(tips)
MSE_baseline
# -

# RMSE of OLS model
RMSE = mean_squared_error(tips.tip, tips.yhat, squared = False)
RMSE

# RMSE for the baseline model
RMSE_baseline = mean_squared_error(tips.tip, tips.baseline, squared = False)
RMSE_baseline

# +
# ESS = sum(tips.yhat - tips.tip.mean())**2

ESS = sum((tips.yhat - tips.baseline)**2)
ESS

# +
# Total Sum of Errors

TSS = ESS + SSE
TSS
# -



# ### 5. Write python code that compares the sum of squared errors for your model against the sum of squared errors for the baseline model and outputs whether or not your model performs better than the baseline model.

# +
df_eval = pd.DataFrame(np.array(['SSE', 'MSE','RMSE']), columns=['metric'])
df_eval['model_error'] = np.array([SSE,  MSE, RMSE])


df_eval

# +
df_eval['baseline_error'] = np.array([SSE_baseline,MSE_baseline, RMSE_baseline])

df_eval
# -

df_eval['better_than_baseline'] = df_eval.baseline_error > df_eval.model_error

df_eval

# ### 7. What is the amount of variance explained in your model?

from sklearn.metrics import r2_score
r2_score(tips.tip, tips.yhat)


# ### 9. Create a file named evaluate.py that contains the following functions.

def plot_residuals(actual, predicted):
    residuals = actual - predicted
    plt.hlines(0, actual.min(), actual.max(), ls=':')
    plt.scatter(actual, residuals)
    plt.ylabel('residual ($y - \hat{y}$)')
    plt.xlabel('actual value ($y$)')
    plt.title('Actual vs Residual')
    plt.show()


# +
def residuals(actual, predicted):
    return actual - predicted

def sse(actual, predicted):
    return (residuals(actual, predicted) **2).sum()

def mse(actual, predicted):
    n = actual.shape[0]
    return sse(actual, predicted) / n

def rmse(actual, predicted):
    return math.sqrt(mse(actual, predicted))

def ess(actual, predicted):
    return ((predicted - actual.mean()) ** 2).sum()

def tss(actual):
    return ((actual - actual.mean()) ** 2).sum()

def r2_score(actual, predicted):
    return ess(actual, predicted) / tss(actual)


# +
def regression_errors(actual, predicted):
    return pd.Series({
        'sse': sse(actual, predicted),
        'ess': ess(actual, predicted),
        'tss': tss(actual),
        'mse': mse(actual, predicted),
        'rmse': rmse(actual, predicted),
    })

def baseline_mean_errors(actual):
    predicted = actual.mean()
    return {
        'sse': sse(actual, predicted),
        'mse': mse(actual, predicted),
        'rmse': rmse(actual, predicted),
    }

def better_than_baseline(actual, predicted):
    rmse_baseline = rmse(actual, actual.mean())
    rmse_model = rmse(actual, predicted)
    return rmse_model < rmse_baseline



# -

better_than_baseline(tips.tip, tips.yhat)

# ### 10. Load the mpg dataset and fit a model that predicts highway mileage based on engine displacement. Take a look at all the regression evaluation metrics, and determine whether this model is better than the baseline model. Use the functions from your evaluate.py to help accomplish this.

mpg = data('mpg')

mpg.head()

# +
#plot displacement vs highway mpg

plt.scatter(mpg.displ, mpg.hwy)

# +
# create the model object
lm = LinearRegression(normalize=True)

# fit the model to trainig data
lm.fit(mpg[['displ']], mpg.hwy)

# make prediction
predictions = lm.predict(mpg[['displ']])

# +
# plot regression line

plt.scatter(mpg.displ, mpg.hwy)
plt.plot(mpg.displ, predictions)
plt.xlabel('displacement (lites)')
plt.ylabel('highway mpg')

# +
# plot displacement vs residuals

plt.scatter(mpg.displ, (mpg.hwy - predictions))
plt.axhline(0, ls = ':')
plt.ylabel('residuals')
plt.xlabel('displacement (liters)')
# -

mpg[mpg.displ > 5]

# calculate regressions errors
evaluate.regression_errors(mpg.hwy, predictions)

# +
# is our model better than baseline?

evaluate.better_than_baseline(mpg.hwy, predictions)

# +
# R2 score

evaluate.r2_score(mpg.hwy, predictions)
# -



StarcoderdataJupyterScriptsDedupFiltered