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

from datetime import timedelta, date
from collections import Counter 
import numpy as np 
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt 
from plotly.offline import init_notebook_mode, iplot 
import plotly.graph_objs as go
from plotly import tools
from sklearn.linear_model import LinearRegression
init_notebook_mode(connected=True) 
# %matplotlib inline 

# ### Approach #2: Trigonometric regression model  
#
# Part 2 continues the attempt to reduce the number of parameters in our equation, but this time using a model based approach rather than a Fourier transform.
#
# We start by reproducing our time series.

# +
def is_prime_wed(date):
    return (date.day in {2,3,5,7,11,13,17,19,23,29,31}) & (date.weekday() == 2)

def count_prime_weds(year):
    start_date = date(year, 1, 1)
    end_date = date(year, 12, 31)
    
    dates = [start_date + timedelta(n) for n in range((end_date - start_date).days + 1)] 
    
    return sum(1 for date in dates if is_prime_wed(date))


# +
start_year = 1912
end_year = 2233
years = range(start_year, end_year + 1)

prime_weds = [count_prime_weds(year) for year in range(start_year, end_year + 1)]
# -

plt.figure(figsize=(15, 1))
plt.plot(years, prime_weds) 
plt.show()

# Re-creating scaled dependant/independant variables for our model, this time we will double the size of our input data to avoid issues with n-observations ~ n-features.

y = np.tile(np.array(prime_weds), 2)
sf = end_year - start_year + 1  # sample frequency
x = np.arange(len(y)) / sf


# Currently we only have one independant variable for our model: `x` which represents the year (but has been scaled between 0-1)
#
# As the relationship between `x` and `y` is clearly non-linear we will get a terrible fit using just one variable. 
#
# To generate more independant variables we can transform `x`. We could create a polynomial series e.g. `c₁x + c₂x^2 + c₃x^3...` but for this task it makes more sense to use a trigonometric series `c₁sin(x) + c₂sin(2x) + c₃sin(3x)...`, for the reasons stated in part 1 of the notebook.
#
# Here, each `c` is a coefficient which we will try to optimise, and each transform of x is a parameter in our model. We are trying to produce the simplest model, that is, the one with the fewest parameters.

def sine_transform(x, freq):
    ''' generate sine waves of specified frequencies from time series array '''
    return np.sin(freq * x * 2*np.pi)


# Here we have defined our sine transform which we use to "generate" features for our model. Each feature will be the same transform with a different frequency, and the amplitude of each frequency will be the coefficient we will try to learn in the model. Plotting examples of different frequencies below (fixed amplitude)

# +
plt.figure(figsize=(15,2))

freqs = [25, 5, 1]
for freq in freqs:
    plt.plot(x, sine_transform(x, freq), label=str(freq))
    plt.legend(title='Frequencies', loc='lower left')
plt.show()


# -

# We can also define a cosine transformation, this is the same as sine but with a 90 degrees phase shift and should hopefully give us some more flexibility when we try to fit the model 

def cosine_transform(x, freq):
    ''' generate cosine waves of specified frequencies from time series array '''
    return np.cos(freq * x * 2*np.pi)


plt.figure(figsize=(15,2))
plt.plot(x, sine_transform(x, 5), label=str(5) + ' (sine)')
plt.plot(x, cosine_transform(x, 5), label=str(5) + ' (cosine)')
plt.legend(title='Frequencies', loc='lower left')
plt.show()

# Generating input data for model, by defining sine and cosine transformations for frequencies between 2 and 161 (nyquist). We start at 2 as the model will have a perfect fit otherwise and the t-values will be meaningless (and cannot be used to trim the model) 

freqs = range(1, 162)
transforms = []
feature_names = []
for freq in freqs:
    transforms.append(sine_transform(x, freq))
    feature_names.append('sin' + str(freq))
    transforms.append(cosine_transform(x, freq))
    feature_names.append('cos' + str(freq))

# Add "bias" (intercept) term

transforms.insert(0, np.ones(len(x)))
feature_names.insert(0, 'bias')

# Define final array containing all features

X = np.array(transforms).T
X.shape

# We now want to try fitting our model, to make sure none of our resiudal exceed the required threshold (0.5). We initially try an Ordinary Least Squares fit.

# +
df = pd.DataFrame(X, columns=feature_names)

results = sm.OLS(y, df).fit(method="qr")
abs(results.resid).max()


# -

# We can see our max residual is clearly < 0.5, will now plot the resisudals for all observations.

def plot_residuals(x, residuals):
    plt.figure(figsize=(15, 1))
    plt.plot(x, residuals)
    plt.ylim((-0.5, 0.5))
    plt.show()


plot_residuals(x, results.resid)


# It is not suprising that we have such a close fit in our model, we are using over 300 parameters!
#
# We will now try reducing the complexity of the model by fitting it in a step-wise loop, each time removing the parameter which corresponds to the lowest t-value and then checking our new model does not exceed our required threshold.

def fit_threshold(df, y, threshold=0.5):
    ''' fit model recursively until we exceed threshold, then return previous result '''
    
    results = sm.OLS(y, df).fit(method="qr") # need qr, otherwise incorrect result (see "bug" notebook) 
    max_residual = abs(results.resid).max()
    
    while max_residual < threshold:
        exclude_feature = results.tvalues[results.tvalues.index != 'bias'].abs().idxmin()
        print('Features:', len(df.columns), '| Max residual: ', max_residual)
        df = df[df.columns.difference([exclude_feature])]
        next_results = fit_threshold(df, y)
        return next_results or results
    
    return None 


# +
df = pd.DataFrame(X, columns=feature_names)

results = fit_threshold(df, y, 0.5)
# -

abs(results.resid).max()

plot_residuals(x, results.resid)

# Our new result is the model with the fewest possible features using the above rule i.e. if we were to remove one more feature then we would exceed our threshold.
#
# We now only have **150** parameters in our model, which is a big improvement, but is still not particularly impressive as we only have 161 values in our input x. This means we could just effectively create a model with a one-to-one mapping between each x to each y using only 11 more parameters.
#
# In addition we showed in Part One of the notebook that we can create a "model" with only 111 parameters and not exceed the residual threshold.
#
# One thing we could try to improve our result is using a different cost function, other than least squares. Least squares will penalise outliers more heavily, and will minimise the sum of the error over all observations. In reality we only care about the maximum error from any one point. A brief look into trying to minimize the maximum error flagged up a few things to look into: L-infinity norm, minimax, chebychev approximation theory. 
#
# #### To be continued...



StarcoderdataJupyterScriptsDedupFiltered