Modelling Lottery Number Popularity

import pandas as pd

def pre19import(yearstring):
    df = pd.read_excel('LOTTO6aus49.xls', sheet_name=yearstring, skiprows=6, parse_dates=[0])
    df = df.drop(range(52,60),axis=0).copy()
    df.rename(columns = {'Unnamed: 1':'Date'}, inplace = True) 
    pd.to_datetime(df['Date']).apply(lambda x: x.date())
    df.set_index(["Date"], inplace = True)
    df.drop(df.columns[[0,1]], axis=1,inplace=True)
    df['Quote Kl. 9'] = 5
    df.columns = ['G1', 'G2', 'G3', 'G4', 'G5', 'G6', 'S', 'Spieleinsatz', 'Anz. Kl. 1',
           'Quote Kl. 1', 'Anz. Kl. 2', 'Quote Kl. 2', 'Anz. Kl. 3', 'Quote Kl. 3',
           'Anz. Kl. 4', 'Quote Kl. 4', 'Anz. Kl. 5', 'Quote Kl. 5', 'Anz. Kl. 6',
           'Quote Kl. 6', 'Anz. Kl. 7', 'Quote Kl. 7', 'Anz. Kl. 8', 'Quote Kl. 8',
           'Anz. Kl. 9', 'Quote Kl. 9']
    df = df.replace({'Jackpot': 0, 'JP': 0, 'unbesetzt': 0}, regex=False)
    for c in df.columns:
        df[c] =  df[c].astype(float)
    for c in ['G1', 'G2', 'G3', 'G4', 'G5', 'G6', 'S']:
        df[c] =  df[c].astype(int)
    return df

df2017 = pre19import('2017')
df2018 = pre19import('2018')

df2017.head(5)

df19 = pd.read_excel('LOTTO_ab_2019.xls',  index_col=0, parse_dates=[0], decimal=",", thousands='.')
df19.drop(df19.columns[[0,7,8,10,11]], axis=1,inplace=True)
df19 = df19.iloc[0:-4,:].copy()
df19.columns = ['G1', 'G2', 'G3', 'G4', 'G5', 'G6', 'S', 'Spieleinsatz', 'Anz. Kl. 1',
           'Quote Kl. 1', 'Anz. Kl. 2', 'Quote Kl. 2', 'Anz. Kl. 3', 'Quote Kl. 3',
           'Anz. Kl. 4', 'Quote Kl. 4', 'Anz. Kl. 5', 'Quote Kl. 5', 'Anz. Kl. 6',
           'Quote Kl. 6', 'Anz. Kl. 7', 'Quote Kl. 7', 'Anz. Kl. 8', 'Quote Kl. 8',
           'Anz. Kl. 9', 'Quote Kl. 9']
df19 = df19.replace({'Jackpot': 0, 'JP': 0, 'unbesetzt': 0}, regex=False)
for c in ['G1', 'G2', 'G3', 'G4', 'G5', 'G6', 'S']:
    df19[c] =  df19[c].astype(int)

df19.shape

df = pd.concat([df2017, df2018, df19])
print(df.shape)

df.describe()

import numpy as np
#dfclass8 = df.iloc[:,np.r_[0:6,7,23:25]]
dfclass8 = pd.DataFrame()
dfclass8 = df.iloc[:,np.r_[0:6,23:24]].copy()
dfclass8.rename(columns = {'Quote Kl. 8':'quota_EUR_8'}, inplace=True)

# Winning class 6 + 8
dfclass6_8 = df.iloc[:,np.r_[0:6,19:20,23:24]].copy()
dfclass6_8.rename(columns = {'Quote Kl. 6':'quota_EUR_6'}, inplace=True)
dfclass6_8.rename(columns = {'Quote Kl. 8':'quota_EUR_8'}, inplace=True)

# Winning class 7 for 'Superzahl' analysis
dfclass7 = df.iloc[:,np.r_[6,21:22]].copy()
dfclass7.rename(columns = {'Quote Kl. 7':'quota_EUR_7'}, inplace=True)
dfclass7.rename(columns = {'S':'Superzahl'}, inplace=True)

# Winning class 5 for 'Superzahl' analysis
dfclass5 = df.iloc[:,np.r_[6,17:18]].copy()
dfclass5.rename(columns = {'Quote Kl. 5':'quota_EUR_5'}, inplace=True)
dfclass5.rename(columns = {'S':'Superzahl'}, inplace=True)

def sparsedf(df, numbers=6, pop=49): # custom hot encoding
    for i in range(1,pop+1):
        df.loc[:,'N' + str(i)] = pd.Series(np.zeros((len(df)))) #sparse cols
    df = df.fillna(0) # zero prefilled
    num = range(1,numbers+1)
 
    for index, row in df.iterrows(): # each drawing
        for h in num: # each drawn number
            dn = df.iloc[index,h].copy()
            df.loc[index,'N' + str(dn)] = 1 # sparsing
    return df

dfclass8.reset_index(inplace=True)
dfclass8 = sparsedf(dfclass8)

dfclass8

# expecting a uniform distribution of drawn numbers across population 1 -49
dfclass8.sum()[7:58]

# plotting
import numpy as np
import pandas as pd
from pandas import Series, DataFrame
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
plt.rcParams['figure.figsize'] = [12, 8] # for square canvas
plt.rcParams.update({'font.size': 16})

data = [23, 45, 56, 78, 213]
plt.bar(range(len(dfclass8.sum()[7:58])), dfclass8.sum()[7:58], color='blue')
nn = np.linspace(1,50)
plt.plot(range(len(dfclass8.sum()[7:58])),[37.8367 for w in dfclass8.sum()[7:58]],'r')
plt.title('Frequency of each number drawn 2017-2019')
figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.show()
print(dfclass8.sum()[7:58].mean())

# goodness of fit - test

# frequency of numbers drawn
number_draws_observed = dfclass8.sum()[7:58].to_numpy()

# 1854 numbers drawn in 309 drawing sessions across the 49 number population, 
# with an expected mean of each number to show up of 37.83673469387755 times (see also appendix)
# create a strictly uniform expected dist.
number_draws_expected = [38.44897959183673 for w in dfclass8.sum()[7:58]]
#print('number_draws_expected' ,number_draws_expected)
#print('   ')

# compare both distributions and compare
from scipy.stats import chisquare, chi2_contingency
test = chisquare(f_obs=number_draws_observed, f_exp=number_draws_expected)
print(test)
print(round(test[1],3))

dfclass8_X = pd.DataFrame(dfclass8.iloc[:,8:58])
dfclass8_y = dfclass8.iloc[:,7].to_numpy()

dfclass7_X = pd.DataFrame(dfclass7.iloc[:,2:12])
dfclass7_y = dfclass7.iloc[:,0].to_numpy()

from sklearn.model_selection import train_test_split
SEED = 1
# As usually we split the data into a training and test data set (X_train, y_train, and X_test, y_test)
# winning tier 8 ( three correct numbers )
X_train, X_test, y_train, y_test = train_test_split(dfclass8_X, dfclass8_y, test_size=0.11, random_state=SEED)
# winning tier 7 ( three correct numbers and correct Superzahl )
X_train_7, X_test_7, y_train_7, y_test_7 = train_test_split(dfclass7_X, dfclass7_y, test_size=0.15, random_state=SEED)

print('Train data shape 8 / 7:', X_train.shape, X_train_7.shape)
print('Test data shape 8 / 7: ', X_test.shape, X_test_7.shape)

#imports
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV
from sklearn.svm import LinearSVR
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import ElasticNet
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.svm import SVR
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import BayesianRidge
from sklearn.ensemble import AdaBoostRegressor
from sklearn.linear_model import OrthogonalMatchingPursuit
from sklearn.linear_model import ARDRegression
from sklearn.linear_model import SGDRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.svm import SVR
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.kernel_ridge import KernelRidge

def warn(*args, **kwargs):
    pass
import warnings
warnings.warn = warn

#model definitions
all_param = pd.Series()

queue = [
       
         ( 
            'linreg',
            LinearRegression(copy_X=True),
            {},
            'coef_'
         ),
        (
            'ridge',
            Ridge(random_state=SEED,
             copy_X=True
            ),
            {
                'alpha': [1, 1.25, 1.5, 1.75,],
                'max_iter': [ 200,300, 400, 500,600, 1000, 2000],
                'solver': ['auto', 'svd', 'cholesky', 'lsqr', 'sparse_cg', 'sag', 'saga'],
             },
            'coef_'
        ),
    (
            'lasso',
            Lasso(random_state=SEED,
             #copy_X=True
            ),
                {
                'alpha': [0.001, 0.0015,0.002, 0.0025, 0.003, 0.004, 0.005, 0.0075,  ],
                'max_iter': [  300, 400,500,700,900],
             },
            'coef_'
    ),
    (
            'bayesridge',
            BayesianRidge(
                copy_X=True,
                lambda_1=0.2,
            ),

                {
                'alpha_1': [
                                    0.001, 0.0015,0.002, 0.0025, 0.003, 0.0004, 0.0005, 0.000001, 0.0000001,
                                    #0.25, 0.5, 0.75, 1, 1.5, 2, 3, 5,
                                    #0.6
                                ],
                'lambda_1': [
                                    0.001, 0.0015,0.002, 0.0025, 0.003, 0.0004, 0.0005, 0.000001,
                                ],
                'n_iter': [100,200,300,400,500,1000],
             },
            'coef_'
    ),
    (
            'sgd',
            SGDRegressor(
                random_state=SEED
            ),

                {
                'alpha': [
                                    0.001, 0.0015,0.002, 0.0025, 0.003, 0.001, 0.0004, 0.0005, 0.000001,
                                    #0.25, 0.5, 0.75, 1, 1.5, 2, 3, 5,
                                    #0.6
                                ],
                'l1_ratio': [
                                    0.005,0.05, 0.1, 0.15, 0.25, 0.5, 0.75, 1,
                                ],
                'max_iter': [300,400,500,1000, 2000],
             },
            'coef_'
    ),
    (
            'linsvr',
            LinearSVR(
               
            ),

                {
                 
                'epsilon': [
                                    0.1, 0.125, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6,0.075,  0.05
                                    #0.25, 0.5, 0.75, 1, 1.5, 2, 3, 5,
                                    #0.6
                                ],
                'tol': [
                                    0.025,0.04,0.05, 0.06,0.07,0.875, 0.1,0.15
                                ],
               
                'C': [ 0.5, 0.7 , 0.9, 1,1.1,1.2,1.5, 2],
                 'max_iter': [ 400,700,1000,1200,1500,2000,3000],    
             },
            'coef_'
    )
]

%%time
#gs
from loky import get_reusable_executor
get_reusable_executor(max_workers=4, timeout=500)
                           
for name, model, param, attribute in queue[0:6]:
    
    gs = GridSearchCV(model,param,cv=10,n_jobs=-1,scoring='neg_mean_squared_error' )
    gs.fit(X_train, y_train)
    all_param[name] = gs.best_params_
    print(gs.best_score_, name, gs.best_params_)

print(' ')
print(all_param)

print(len(queue))

from sklearn import metrics

model_weights = pd.DataFrame({'number': range(1,50)})

for name, model, param, attribute in queue[0:6]:
    
    regr = model
    regr.set_params(**all_param[name]) # fit best scoring grid searched model
       
    regr.fit(X_train, y_train)
    
    y_pred = regr.predict(X_test)
    
    print(name, ' --------------')

    print('Mean Squared Error:', mean_squared_error(y_test, y_pred))  
    print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, y_pred)))
    
    model_weights[name] = getattr(regr,attribute)

model_weights =  model_weights.iloc[:,1:5]

# normalization
from sklearn.preprocessing import StandardScaler
from sklearn import preprocessing

weights = model_weights.iloc[:,:].copy()
scaler = StandardScaler()

for n in weights.columns:
    np_array = np.array(weights[n].values)
    arr = scaler.fit_transform(np_array.reshape(-1, 1))
    weights[n] = arr

# plotting popularity

import pandas as pd

model_weights = pd.DataFrame({'number': range(1,50)})
model_weights['number'] = model_weights['number'].astype(str) #convert
#model_weights.reindex([model_weights['number'])

model_weights = pd.concat([model_weights, weights], axis=1)
model_weights.set_index('number')
model_weights.index += 1 
print('Having a peak into the data...')
print(model_weights.head(7))

#prepare two df for df.plot   
dfplot_mw1 = model_weights[0:25][['linreg','ridge','lasso', 'bayesridge']].copy()
dfplot_mw2 = model_weights[25:50][['linreg','ridge','lasso', 'bayesridge']].copy()

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 16})

a = dfplot_mw1.plot(kind='bar', 
                    title ="Modeled Number Feature Coefficents (positive = higher quota / negative = lower quota)", figsize=(17, 9), legend=True, fontsize=14)
plt.xticks(rotation=0)

b = dfplot_mw2.plot(kind='bar', 
                    title ="Modeled Number Feature Coefficents (positive = higher quota / negative = lower quota)", figsize=(17, 9), legend=True, fontsize=14)
ax = [a,b]

for axis in ax: 
    axis.set_xlabel("number", fontsize=14)
    axis.set_ylabel("value", fontsize=14)
    
    major_ticks = np.arange(0, 24, 10)
    minor_ticks = np.arange(0, 24, 5)

    axis.grid(which='major', color='#CCCCCC', linestyle='--')
    axis.grid(which='minor', color='#CCCCCC', linestyle=':')

    axis.grid(True)


plt.xticks(rotation=0)
plt.show()

mean_weights = pd.DataFrame({'number': range(1,50)})
mean_weights.index += 1 
mean_weights['mean'] = model_weights.mean(axis=1)
sorted_mean_weights = mean_weights.sort_values('mean')

axx = sorted_mean_weights[0:50][['mean']].plot(kind='bar', title ="Sorted Modeled Number Feature Coefficents (positive = higher quota / negative = lower quota)", figsize=(16, 9), legend=True, fontsize=12)
plt.xticks(rotation=0, fontsize=13)
plt.xlabel("number", fontsize=16)
plt.ylabel("value", fontsize=16)
plt.show()

mean_weights.sum()

mean_weights.drop(mean_weights.columns[[0]], axis=1,inplace=True)

size = 7
data = mean_weights.to_numpy().reshape((size, size))
values = mean_weights.to_numpy().reshape((size, size))
idx = mean_weights.index.to_numpy().reshape((size, size))

# plot based on https://stackoverflow.com/a/45989759 Ramon Martinez
font = {'family' : ['Tahoma', 'DejaVu Sans', 'Lucida Grande', 'Verdana'],
        'weight' : 'normal',
        'size'   : 13}

plt.rc('font', **font)

# Limits for the extent
x_start = 1
x_end = 8
y_start = 8
y_end = 1

extent = [x_start, x_end, y_start, y_end]

# The normal figure
fig = plt.figure(figsize=(13, 9))
ax = fig.add_subplot(111)

im = ax.imshow(data, extent=extent, origin='upper', interpolation='None', cmap='hot')

# Add the text
jump_x = (x_end - x_start) / (2.0 * size)
jump_y = (y_end - y_start) / (2.0 * size)
x_positions = np.linspace(start=x_start, stop=x_end, num=size, endpoint=False)
y_positions = np.linspace(start=y_end, stop=y_start, num=size, endpoint=False)

for x_index, x in enumerate(x_positions):
    for y_index, y in enumerate(y_positions):
        label = np.round(values[y_index, x_index],2)
        label2 = idx[y_index, x_index]
        text_x = x + jump_x
        text_x2 = x + jump_x
        text_y = y - jump_y*1.35
        text_y2 = y - jump_y*0.75
        if label < 0:
            txtcol = 'white'
        else:
            txtcol = 'black'
        ax.text(text_x, text_y, label, color=txtcol, ha='center', va='center')
        ax.text(text_x2, text_y2, label2, color=txtcol, ha='center', va='center').set_fontsize(18)

# Labels for major ticks
ax.set_xticklabels(np.arange(1, 8, 1));
ax.set_yticklabels(np.arange(1, 8, 1));

# Gridlines based on minor ticks
gridcol= 'black'
ax.grid(which='major', color=gridcol, linestyle='-', linewidth=2)
ax.spines['bottom'].set_color(gridcol)
ax.spines['top'].set_color(gridcol) 
ax.spines['right'].set_color(gridcol)
ax.spines['left'].set_color(gridcol)
ax.tick_params(axis='x', colors=gridcol)
ax.tick_params(axis='y', colors=gridcol)
ax.set_title('\'Lotto 6 aus 49\' - mean quota / number correlations in winning tier 8 across all models\n')

fig.colorbar(im)
plt.setp(ax.get_xticklabels(), visible=False)
plt.setp(ax.get_yticklabels(), visible=False)
#plt.tight_layout(pad=0)
plt.show()

dfcl7mean = dfclass7.groupby('Superzahl').mean()
dfcl5mean = dfclass5.groupby('Superzahl').mean()
import matplotlib.cm as cm

ax = dfcl7mean[0:10].plot(kind='bar', title ="Mean quotas winning tier 7", figsize=(14, 8), 
                             legend=True, fontsize=12)

for p in ax.patches: 
    ax.annotate(np.round(p.get_height(),decimals=2), (p.get_x()+p.get_width()/2., p.get_height()), 
                                 ha='center', va='center', color='w', xytext=(0, -15), textcoords='offset points')
ax.legend(["winning tier 7"]);  
plt.xticks(rotation=0)
plt.show()
dfcl7mean.mean()

from sklearn.preprocessing import StandardScaler, RobustScaler, QuantileTransformer
from sklearn import preprocessing
from decimal import Decimal

dfcl7mean['quota_EUR_5'] = dfcl5mean['quota_EUR_5']
scaler = StandardScaler(with_std=True)
scaler.fit(dfcl7mean)
dfcl7mean_scaled = pd.DataFrame(scaler.transform(dfcl7mean), columns = dfcl7mean.columns)
ax = dfcl7mean_scaled[0:10].plot(kind='bar', title ="Studentized expected quotas for each Superzahl", figsize=(14, 8), legend=True, fontsize=12)

# from https://stackoverflow.com/a/48372659 thank you 'justfortherec'
def add_value_labels(ax, spacing=5):
    """Add labels to the end of each bar in a bar chart.

    Arguments:
        ax (matplotlib.axes.Axes): The matplotlib object containing the axes
            of the plot to annotate.
        spacing (int): The distance between the labels and the bars.
    """

    # For each bar: Place a label
    for rect in ax.patches:
        # Get X and Y placement of label from rect.
        y_value = rect.get_height()
        x_value = rect.get_x() + rect.get_width() / 2

        # Number of points between bar and label. Change to your liking.
        space = spacing
        # Vertical alignment for positive values
        va = 'bottom'

        # If value of bar is negative: Place label below bar
        if y_value < 0:
            # Invert space to place label below
            space *= -1
            # Vertically align label at top
            va = 'top'

        # Use Y value as label and format number with one decimal place
        label = "{:.1f}".format(y_value)

        # Create annotation
        ax.annotate(
            label,                      # Use `label` as label
            (x_value, y_value),         # Place label at end of the bar
            xytext=(0, space),          # Vertically shift label by `space`
            textcoords="offset points", # Interpret `xytext` as offset in points
            ha='center',                # Horizontally center label
            va=va)                      # Vertically align label differently for
                                        # positive and negative values.


# Call the function above. All the magic happens there.
add_value_labels(ax)
ax.legend(["winning tier 7", "winning tier 5"]);
plt.xticks(rotation=0)
plt.show()

%%time
dfsuperzahl = pd.DataFrame(df[['S','Anz. Kl. 7', 'Anz. Kl. 8','Anz. Kl. 5', 'Anz. Kl. 6','Anz. Kl. 3', 'Anz. Kl. 4']])
#print(dfsuperzahl['Anz. Kl. 7'])

def threeright(row):
    return row['Anz. Kl. 7'] / row['Anz. Kl. 8'] * 10

def fourright(row):
    return row['Anz. Kl. 5'] / row['Anz. Kl. 6'] * 10

def fiveright(row):
    return row['Anz. Kl. 3'] / row['Anz. Kl. 4'] * 10

dfsuperzahl['quotient_threeright'] = df.apply(threeright, axis=1)
dfsuperzahl['quotient_fourright'] = df.apply(fourright, axis=1)
dfsuperzahl['quotient_fiveright'] = df.apply(fiveright, axis=1)
dfsuperzahl.drop(dfsuperzahl.columns[[1,2,3,4,5,6]], axis=1,inplace=True)
#print(dfsuperzahl)
dfsz = pd.DataFrame()
dfsz = dfsuperzahl.groupby('S').mean()
dfsz

# standarize values with a mean of 0

from sklearn.preprocessing import StandardScaler, RobustScaler, QuantileTransformer
scaler = StandardScaler(with_std=True)
scaler.fit(dfsz)
dfsz_scaled = pd.DataFrame(scaler.transform(dfsz), columns = dfsz.columns)

ax = dfsz_scaled.plot(kind='bar', title ="Superzahl popularities", figsize=(16, 8), legend=True, fontsize=12)

# from https://stackoverflow.com/a/48372659 thank you 'justfortherec'
def add_value_labels(ax, spacing=5):
    """Add labels to the end of each bar in a bar chart.

    Arguments:
        ax (matplotlib.axes.Axes): The matplotlib object containing the axes
            of the plot to annotate.
        spacing (int): The distance between the labels and the bars.
    """

    # For each bar: Place a label
    for rect in ax.patches:
        # Get X and Y placement of label from rect.
        y_value = rect.get_height()
        x_value = rect.get_x() + rect.get_width() / 2

        # Number of points between bar and label. Change to your liking.
        if (x_value % 1) == 0.0:
            space = spacing * 5
        else:
            space = spacing
        # Vertical alignment for positive values
        va = 'bottom'

        # If value of bar is negative: Place label below bar
        if y_value < 0:
            # Invert space to place label below
            space *= -1
            # Vertically align label at top
            va = 'top'

        # Use Y value as label and format number with one decimal place
        label = "{:.2f}".format(y_value)

        # Create annotation
        ax.annotate(
            label,                      # Use `label` as label
            (x_value, y_value),         # Place label at end of the bar
            xytext=(0, space),          # Vertically shift label by `space`
            textcoords="offset points", # Interpret `xytext` as offset in points
            ha='center',                # Horizontally center label
            va=va)                      # Vertically align label differently for
                                        # positive and negative values.


# Call the function above. All the magic happens there.
add_value_labels(ax)
plt.xticks(rotation=0)
#plt.legend(loc=4)
plt.show()

# Split Data with numbers and superzahl combined into Train and Test Data
dfclass8sz = pd.DataFrame()
dfclass8sz = df.iloc[:,np.r_[0:7,23:24]].copy()
dfclass8sz.reset_index(inplace=True)
dfclass8sz.rename(columns = {'Quote Kl. 8':'quota_EUR_8'}, inplace=True)
dfclass8sz = sparsedf(dfclass8sz)
del dfclass8sz['index']
dfclass8sz['S'] = pd.Categorical(dfclass8sz.S)
dfclass8sz = pd.get_dummies(dfclass8sz, prefix_sep='')
dfclass8sz_X = pd.DataFrame(dfclass8sz.iloc[:,7:67]).copy()
dfclass8sz_y = dfclass8sz.iloc[:,6].to_numpy()
X_train, X_test, y_train, y_test = train_test_split(dfclass8sz_X, dfclass8sz_y, test_size=0.15, random_state=SEED)

%%time
all_param_sz = pd.Series()

#gs
from loky import get_reusable_executor
get_reusable_executor(max_workers=4, timeout=500)
                           
for name, model, param, attribute in queue[0:6]:
    
    gs = GridSearchCV(model,param,cv=5,n_jobs=1,scoring='neg_mean_squared_error' )
    gs.fit(X_train, y_train)
    all_param_sz[name] = gs.best_params_
    print(gs.best_score_, name, gs.best_params_)

print(' ')
print(all_param_sz)

# add non-linear high performing algos
svr_param = {'C': 0.075, 'coef0': 5, 'degree': 3, 'epsilon': 0.075, 'kernel': 'poly'}
kr_param = {'alpha': 0.6, 'coef0': 2, 'degree': 7, 'kernel': 'polynomial'}

models = [

    ('linreg', LinearRegression(**all_param_sz['linreg'])),
    
    ('lasso', Lasso(random_state=SEED,
          **all_param_sz['lasso']
            )),
    
    ('ridge', Ridge(random_state=SEED,
             copy_X=True,
            **all_param_sz['ridge'])),
   
    ('bridge', BayesianRidge(**all_param_sz['bayesridge'])),
    
    # add non-linear high performing algos
    ('svr', SVR(**svr_param)),
    
    ('kridge',KernelRidge(**kr_param))

]

from sklearn import metrics

for name, model in models:
    
    regr = model
    regr.fit(X_train, y_train)
    
    y_pred = regr.predict(X_test)
    
    print(name, '-------')
    print('Mean Squared Error:', mean_squared_error(y_test, y_pred))  
    print('Root Mean Squared Error:', np.sqrt(mean_squared_error(y_test, y_pred)))
    print('--------------')

%%time
from vecstack import StackingTransformer
stack = StackingTransformer(models,                     # list of models
                           regression=True,            # regression task (if you need 
                                                       #     classification - set to False)
                           metric=mean_squared_error, # metric: callable
                           n_folds=10,                  # number of folds
                           shuffle=True,               # shuffle the data
                           random_state=SEED,             # ensure reproducibility
                           verbose=0)                  # print all info

# Fit
stack = stack.fit(X_train, y_train)

# Get your stacked features
S_train = stack.transform(X_train)
S_test = stack.transform(X_test)

S_train[0:3]

S_test[0:3]

def rmspe(y_true, y_pred):
    '''
    Compute Root Mean Square Percentage Error between two arrays.
    '''
    return np.sqrt(np.mean(np.square(((y_true - y_pred) / y_true)), axis=0))*100

# Initialize 2nd level model
model_l2 = LinearRegression()
    
# Fit 2nd level model
model_l2 = model_l2.fit(S_train, y_train)

# Predict
y_pred = model_l2.predict(S_test)

# Final prediction score
print('Mean Squared Error: [%.8f]' % mean_squared_error(y_test, y_pred))
print('Root Mean Squared Error: [%.8f]' % np.sqrt(mean_squared_error(y_test, y_pred)))
print('Root Mean Squared Percentage Error: ', round(rmspe(y_test, y_pred),2), '%')
print('R2 score: [%.8f]' % metrics.r2_score(y_test, y_pred))

rmspe_var = round(rmspe(y_test, y_pred),2)
rmse_var = round(np.sqrt(mean_squared_error(y_test, y_pred)),2)

compare = pd.DataFrame({'Actual Testdata Quota': y_test.flatten(), 'Predicted Quota': np.round(y_pred.flatten(),2)})
compare

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
fig = plt.figure(figsize=(16, 9))
import numpy as np

rmse = np.sqrt(mean_squared_error(y_test, y_pred))
#y = np.sin(x) + dy * np.random.randn(50)
yp = compare['Predicted Quota'].to_numpy().flatten()
y = compare['Actual Testdata Quota'].to_numpy().flatten()
x = compare.index

plt.errorbar(x, yp, yerr=rmse, fmt='.k', marker="H", capsize=5, ms=8)
plt.scatter(x, y, c='blue', marker='x', s=100)
plt.legend(["Actual Testdata Quota", "Predicted Quota"], fontsize=17)
plt.title("Predicted and Actual Testdata Quotas for winning tier 8", fontsize=20)
plt.xlabel('Testdata observations', fontsize=16)
plt.ylabel('Quota in EUR', fontsize=16)
plt.show()

svrmodel = SVR(**svr_param)
svrmodel = svrmodel.fit(X_train, y_train)
y_pred = svrmodel.predict(X_test)
print('Mean Squared Error: [%.8f]' % mean_squared_error(y_test, y_pred))
print('Root Mean Squared Error: [%.8f]' % np.sqrt(mean_squared_error(y_test, y_pred)))
print('Root Mean Squared Percentage Error: ', round(rmspe(y_test, y_pred),2), '%')
print('R2 score: [%.8f]' % metrics.r2_score(y_test, y_pred))

###################################################################
#predict the quota for your numbers / or an observed number drawing
my_numbers = [3, 16, 23, 33, 41, 47]
my_superzahl = [4]
date = ('22.01.20')
#predict the quota for your numbers / or an observed number drawing
###################################################################

def sparse_superzahl_df(df, numbers=1, pop=10, start=0): # custom hot encoding
    for i in range(start,pop+start):
        df.loc[:,'S' + str(i)] = pd.Series(np.zeros((len(df)))) #sparse cols
    df = df.fillna(0) # zero prefilled
    num = range(start,numbers+start)
    for index, row in df.iterrows(): # each observation
            dn = df.loc[index,'S'].copy()
            df.loc[index,'S' + str(dn)] = 1 # sparsing
    return df

def DFprediction(my_numbers, my_superzahl=0, date=['dummy']):
    import itertools
    data = [date,my_numbers, my_superzahl]
    data = list(itertools.chain.from_iterable(data))
    cols = ['Date', 'G1', 'G2','G3', 'G4','G5', 'G6', 'S']
    df = pd.DataFrame([data], columns = cols, index=[0])
    df = sparsedf(df) # reg numbers
    df = sparse_superzahl_df(df) # superzahl
    df =  df.iloc[:,8:len(df.columns)]
    X_array = df.to_numpy()
    S_observation = stack.transform(X_array)
    return np.round(model_l2.predict(S_observation),2)

print('The estimated quota for the given numbers and winning tier is: %.2f EUR' % DFprediction(my_numbers, my_superzahl), 
      '+/- RMSE %.2f EUR' % np.sqrt(mean_squared_error(y_test, y_pred)))

# k draws out of n without replacement --> use the random.sample function

# simulate how many drawings?
sampling_cycles = 314

# 6 out of 49 Lotto draw function
def lottodraw():
    import random
    draw = np.array(random.sample(range(1,50), 6))
    return np.sort(draw)

# call cycles
yy = [lottodraw() for w in range(sampling_cycles)]
# show last three lotto draw sessions
yy[306:310]

yy = np.array(yy).flatten()
#frequency of numbers drawn altogether
len(yy)

# break it down to each number in the population
unique_elements, counts_elements = np.unique(yy, return_counts=True)

count = pd.DataFrame(counts_elements.reshape(1,-1), columns=[unique_elements])

# Expected frequency for each uniformly distributed number
counts_elements.mean()