<jupyter_start><jupyter_text>Importing Libraries<jupyter_code>%load_ext autoreload
%autoreload 2
%matplotlib inline

from regression_analysis.fit_model import linear_regression
from regression_analysis.utils import franke
from regression_analysis.utils.plots import triangulation_for_triheatmap as triheatmap
from regression_analysis.fit_model.apply_linear_regression import plot_stat, plot_stat_sgd, apply_regression, apply_regression_sgd
import numpy as np
import matplotlib.pyplot as plt

from matplotlib.tri import Triangulation
from matplotlib.ticker import AutoMinorLocator, MultipleLocator
from matplotlib import cm
import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn import linear_model

import ipywidgets as widget
from ipywidgets import interact, interactive, fixed, interact_manual

import glob as glob
from IPython.display import Image, display, HTML, Video

import multiprocessing as mp

import os
if not os.path.exists('data_linear_regression'):
    os.makedirs('data_linear_regression')
if not os.path.exists('data_linear_regression_sgd'):
    os.makedirs('data_linear_regression_sgd')
# Needed for parallelization later
print("Number of processors: ", mp.cpu_count())<jupyter_output>Number of processors:  8<jupyter_text>Demo for applying regression to Franke data   Visualisation of Franke Data<jupyter_code>n = 100 #number of data points along each axis. Total data points = n*n
x1 = np.linspace(0,1,n)
x2 = np.linspace(0,1,n)
xx1, xx2 = np.meshgrid(x1, x2)
xx1 = xx1.reshape((n*n),1)
xx2 = xx2.reshape((n*n),1)

noise_var = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5])
noise_var = np.arange(0,6)*0.0002

fig, ax = plt.subplots(2, 3,subplot_kw={"projection": "3d"}, figsize=(12,8))
for ind, var in enumerate(noise_var):
    y = franke.Franke(xx1, xx2, noise_var=np.sqrt(var)) #zero mean gaussian noise
    xx1 = xx1.reshape([n,n])
    xx2 = xx2.reshape([n,n])

    y = y.reshape(n,n)
    #rescaling y
    y = (y-np.amin(y))/(np.amax(y)-np.amin(y))
    
    if(ind<3): a = 0 
    else: a = 1

    #ax[a, ind%3].plot_surface(xx1, xx2, y, rstride=8, cstride=8, alpha=0.2, cmap=cm.coolwarm)
    ax[a, ind%3].title.set_text(str(var))
    cset = ax[a, ind%3].contourf(xx1, xx2, y, zdir='z', offset=0, cmap=cm.coolwarm)
    cset = ax[a, ind%3].contour(xx1, xx2, y, zdir='x', offset=0, cmap=cm.coolwarm_r)
    cset = ax[a, ind%3].contour(xx1, xx2, y, zdir='y', offset=1, cmap=cm.coolwarm_r)

n = 100 #number of data points along each axis. Total data points = n*n
x1 = np.linspace(0,1,n)
x2 = np.linspace(0,1,n)
xx1, xx2 = np.meshgrid(x1, x2)
xx1 = xx1.reshape((n*n),1)
xx2 = xx2.reshape((n*n),1)

y = franke.Franke(xx1, xx2, noise_var=0.0) #zero mean gaussian noise has variance = var
print(np.max(y))
#create linear regression object by passing input and output data
linear_reg = linear_regression.linear_regression2D(xx1, xx2, y)<jupyter_output><empty_output><jupyter_text>Apply OLS<jupyter_code>print("own OLS")
linear_reg.apply_leastsquares(order=5, test_ratio=0.1, reg_method="ols")
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)

print('\n'+ "Scikit OLS")
linear_reg.apply_leastsquares(order=5, test_ratio=0.1, reg_method="scikit_ols")
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)
linear_reg.apply_leastsquares_bootstrap(order=5, test_ratio=0.1, n_boots=30, reg_method="ols")
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)<jupyter_output>Train MSE 0.0013923439730719196
Test MSE 0.0014268410392499957
Train R2 0.9748131730949406
Test R2 0.9743388006719232
Train bias 0.05712902034615127
Test bias 0.05753024980678943
Train model variance 0.055690013258889205
Test model variance 0.05603923933975953<jupyter_text>Apply OLS with cross validation sampling<jupyter_code>linear_reg.apply_leastsquares_crossvalidation(order=5, kfolds=10, reg_method="ols")
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)<jupyter_output>Train MSE 0.001402214692532389
Test MSE 0.0014085433301206554
Train R2 0.9747316032316593
Test R2 0.9744605953791732
Train bias 0.005625327734845175
Test bias 0.004864263812675437
Train model variance 0.005485610179965326
Test model variance 0.004798420334117806<jupyter_text>Apply Ridge regression<jupyter_code>print("own ridge")
linear_reg.apply_leastsquares(order=5, test_ratio=0.1, reg_method="ridge", lmbda=0.9)
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)

print('\n'+"scikit ridge")
linear_reg.apply_leastsquares(order=5, test_ratio=0.1, reg_method="scikit_ridge", lmbda=0.9)
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)<jupyter_output>own ridge
Train MSE 0.005645842602311908
Test MSE 0.005858756280193477
Train R2 0.8974437046616717
Test R2 0.9013792344324465
Train bias 0.05505116320842965
Test bias 0.059408558205137255
Train model variance 0.04744486737106172
Test model variance 0.051949014647433005

scikit ridge
Train MSE 0.005640374249617838
Test MSE 0.006089762069647348
Train R2 0.8981461835011143
Test R2 0.8923215383141608
Train bias 0.05537715902987056
Test bias 0.05656646625061504
Train model variance 0.047812272295902486
Test model variance 0.047531258160739255<jupyter_text>Apply Ridge regression with bootstrap<jupyter_code>linear_reg.apply_leastsquares_bootstrap(order=5, test_ratio=0.1, n_boots=30, reg_method="ridge", lmbda=0.9)
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)<jupyter_output>Train MSE 0.005681618690677311
Test MSE 0.005704325899161004
Train R2 0.8977681072876926
Test R2 0.897362709450519
Train bias 0.057422852193160574
Test bias 0.05746212792023058
Train model variance 0.04953835987437372
Test model variance 0.049552357131788044<jupyter_text>Apply Ridge regression with cross validation<jupyter_code>linear_reg.apply_leastsquares_crossvalidation(order=5, kfolds=10, reg_method="ridge", lmbda=0.1)
print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)<jupyter_output>Train MSE 0.003587352408226471
Test MSE 0.003599872334739232
Train R2 0.9353575568181884
Test R2 0.9350153642432717
Train bias 0.005564117476790335
Test bias 0.005419981970152453
Train model variance 0.005011375572516438
Test model variance 0.004949380692700832<jupyter_text>Apply Lasso regression (utilises scikit learn)<jupyter_code>linear_reg.apply_leastsquares(order=5, test_ratio=0.1, reg_method="scikit_lasso", lmbda=0.001)

print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Train bias", linear_reg.trainbias)
print("Test bias", linear_reg.testbias)
print("Train model variance", linear_reg.trainvar)
print("Test model variance", linear_reg.testvar)<jupyter_output>Train MSE 0.009256854773192463
Test MSE 0.008999053959781882
Train R2 0.8341321073038298
Test R2 0.8291186523537881
Train bias 0.05580958212499437
Test bias 0.05266383511240413
Train model variance 0.04052579876935987
Test model variance 0.038035797684885395<jupyter_text>Widgets to display results for various parameter combinations polynomial order: The order of 2D polynomial used for fittingnoise_var = The variance of the gaussian noise added to the Franke dataratio: The testing ratio. For example, if ratio is 0.1 then the test set comprises of 10 percent of the dataset.num: The number of points taken along each direction of the input vector. If num = 20, then the total number of data points is 20*20 = 400stat: Displays the chosen statistic    * training MSE    * testing MSE    * training R2 error    * testing R2 error    * bias in test data    * variance in test predictionmethod: To choose the regression type    - Direct solution and stochastic gradient descent solution available            * OLS        * OLS with bootstrap resampling        * OLS with crossvalidation resampling        * Ridge        * Ridge with bootstrap resampling        * Ridge with crossvalidation resampling            - Direct Solution only (using scikit learn)            * Lasso        * Lasso with bootstrap resampling        * Lasso with crossvalidation resamplingn-boot = number of times bootstrap sampling is performed in the bootstrap method. Changing this for methods that doesn't involve bootstrap sampling will not have any effectk_fold = number of folds in the cross validation method. Changing this for methods that doesn't involve cross validation will not have any effectridge_lambda = Regularisation parameter for ridge regressionlasso_lambda = Regularisation parameter for LASSO regressionlearn_rates = Parameter that controls the descent jump size in gradient descent methodsepoch = Parameter that controls the number of times the algorithm works though the entire dataset in gradient descent methodsnumber of batches = Parameter that controls the batch size in gradient descent methods      Without Stochastic Gradient Descent<jupyter_code>#load data
order = np.load("data_linear_regression/order.npy")
num_points = np.load("data_linear_regression/num_points.npy")
noise_var = np.load("data_linear_regression/noise_var.npy")
test_ratio = np.load("data_linear_regression/test_ratio.npy")
k_folds = np.load("data_linear_regression/k_folds.npy")
n_boots = np.load("data_linear_regression/n_boots.npy")
ridge_lambda = np.load("data_linear_regression/ridge_lambda.npy")
lasso_lambda = np.load("data_linear_regression/lasso_lambda.npy")
methods = ["ols", "ols_bootstrap", "ols_crossvalidation", "ridge", 
           "ridge_bootstrap", "ridge_crossvalidation", "lasso", "lasso_bootstrap",
           "lasso_crossvalidation"]
stats = ["train MSE", "test MSE", "test R2", "train R2", "test bias", "test variance"]

widget.interact(plot_stat, ratio=test_ratio.tolist(), num=num_points.tolist(), stat=stats, 
                method=methods, k_fold=k_folds.tolist(), n_boot=n_boots.tolist(),
                ridge_lmb=ridge_lambda.tolist(), lasso_lmb=lasso_lambda.tolist())<jupyter_output><empty_output><jupyter_text>With Stochastic Gradient Descent<jupyter_code># load data
order = np.load("data_linear_regression_sgd/order.npy")
num_points = np.load("data_linear_regression_sgd/num_points.npy")
noise_var = np.load("data_linear_regression_sgd/noise_var.npy")
test_ratio = np.load("data_linear_regression_sgd/test_ratio.npy")
k_folds = np.load("data_linear_regression_sgd/k_folds.npy")
n_boots = np.load("data_linear_regression_sgd/n_boots.npy")
ridge_lambda = np.load("data_linear_regression_sgd/ridge_lambda.npy")
lasso_lambda = np.load("data_linear_regression_sgd/lasso_lambda.npy")
learn_rates= np.load("data_linear_regression_sgd/learn_rates.npy")
num_min_batches= np.load("data_linear_regression_sgd/num_min_batches.npy")
epochs = np.load("data_linear_regression_sgd/epochs.npy")
methods = [ "ols_sgd", "ols_bootstrap_sgd", "ols_crossvalidation_sgd", "ridge_sgd", "ridge_bootstrap_sgd",
            "ridge_crossvalidation_sgd"]
stats = ["train MSE", "test MSE", "test R2", "train R2", "test bias", "test variance"]

widget.interact(plot_stat_sgd, ratio=test_ratio.tolist(), num=num_points.tolist(), stat=stats, 
                method=methods, k_fold=k_folds.tolist(), n_boot=n_boots.tolist(),
                ridge_lmb=ridge_lambda.tolist(), lasso_lmb=lasso_lambda.tolist(), 
                learn_rate=learn_rates.tolist(), batch=num_min_batches.tolist(), 
                epoch=epochs.tolist())<jupyter_output><empty_output><jupyter_text>Edge cases study of OLS  Worst case studyWhat is the lowest possible MSE that we can get for data normalised to the range [0,1] such that the polynomial fitted is the best polynomial fit?This will happen when half of the data is 0 and the other half is 1, and if we fit a polynomial of order p such that p<<n where n is the number of data points. Then we will end up fitting a constant which is 0.5. What is MSE in this case?<jupyter_code>n = 1000 #number of data points along each axis. Total data points = n*n
x1 = np.linspace(0,1,n)
x2 = np.linspace(0,1,n)
xx1, xx2 = np.meshgrid(x1, x2)
xx1 = xx1.reshape((n*n),1)
xx2 = xx2.reshape((n*n),1)
y = np.zeros(xx1.shape)
y[::2, :]+=1.0
print("mean of the output is ", np.mean(y))

print("fitting a constant")
linear_reg = linear_regression.linear_regression2D(xx1, xx2, y) 
linear_reg.apply_leastsquares(order=0, test_ratio=0.1, reg_method="ols")

print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Test bias", linear_reg.testbias)
print("Test model variance", linear_reg.testvar)

print("fitting a 10th order polynomial")
linear_reg.apply_leastsquares(order=10, test_ratio=0.1, reg_method="ols")

print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Test bias", linear_reg.testbias)
print("Test model variance", linear_reg.testvar)

print("The worst possible MSE for data normalised to range [0,1] is ", 
      linear_reg.testMSE, "and is almost independent of the order of the polynomial used to fit such that the order of polynomial is much less than the number of data points")<jupyter_output>mean of the output is  0.5
fitting a constant
Train MSE 0.24999997945555558
Test MSE 0.25000039034444466
Train R2 -2.220446049250313e-16
Test R2 -8.217832479617826e-06
Test bias 0.25000039034444443
Test model variance 4.930380657631324e-32
fitting a 10th order polynomial
Train MSE 0.2499850967951986
Test MSE 0.25001597682331467
Train R2 5.946842779280814e-05
Test R2 -7.560457756028427e-05
Test bias 0.2500007034197615
Test model variance 1.5161495592950492e-05
The worst possible MSE for data normalised to range [0,1] is  0.25001597682331467 and is almost independent of the order of the polynomial used to fit such that the order of polynomial is much less than the number of data points<jupyter_text>Fitting uniform random numbersWhat will happen if we use OLS to fit data which is essentially a random uniform number between 0 and 1 irrespective of the value of the dependent variables. What is the MSE in this situation?<jupyter_code>n = 1000 #number of data points along each axis. Total data points = n*n
x1 = np.linspace(0,1,n)
x2 = np.linspace(0,1,n)
xx1, xx2 = np.meshgrid(x1, x2)
xx1 = xx1.reshape((n*n),1)
xx2 = xx2.reshape((n*n),1)
y = np.random.uniform(size=xx1.shape)
print("mean of the output is ", np.mean(y))

print("fitting a constant")
linear_reg = linear_regression.linear_regression2D(xx1, xx2, y) 
linear_reg.apply_leastsquares(order=0, test_ratio=0.1, reg_method="ols")

print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Test bias", linear_reg.testbias)
print("Test model variance", linear_reg.testvar)

print("fitting a 10th order polynomial")
linear_reg.apply_leastsquares(order=10, test_ratio=0.1, reg_method="ols")

print("Train MSE", linear_reg.trainMSE)
print("Test MSE", linear_reg.testMSE)
print("Train R2", linear_reg.trainR2)
print("Test R2", linear_reg.testR2)
print("Test bias", linear_reg.testbias)
print("Test model variance", linear_reg.testvar)<jupyter_output>mean of the output is  0.49978607727778634
fitting a constant
Train MSE 0.0832655097108108
Test MSE 0.08327085848269622
Train R2 0.0
Test R2 -7.364200145243416e-06
Test bias 0.08327085848269622
Test model variance 1.1093356479670479e-31
fitting a 10th order polynomial
Train MSE 0.08329262983002073
Test MSE 0.08298314606081367
Train R2 6.373448297614104e-05
Test R2 -5.1651273969000044e-05
Test bias 0.08297895083118705
Test model variance 5.3961492524665985e-06<jupyter_text>Fitting gaussian dataWhat will happen if we use OLS to fit data which is essentially a truncated gaussian number between 0 and 1 where the output is independent of the dependant variables. What is the MSE in this situation?<jupyter_code>n = 500 #number of data points along each axis. Total data points = n*n
x1 = np.linspace(0,1,n)
x2 = np.linspace(0,1,n)
xx1, xx2 = np.meshgrid(x1, x2)
xx1 = xx1.reshape((n*n),1)
xx2 = xx2.reshape((n*n),1)

var = np.linspace(0.025, 1.0, 21)
MSE = np.zeros(len(var))
MSE_full = np.zeros(len(var))
max_gauss = np.zeros(len(var))


for ind, variance in enumerate(var):
    print()
    print("truncated")
    y = np.random.normal(loc=0.5, scale = np.sqrt(variance), size=xx1.shape)
    y[np.where(y>1.0)] = 1.0
    y[np.where(y<0.0)] = 0.0

    print("mean of the output is ", np.mean(y))
    print("variance of the output is ", np.var(y))

    print("fitting a constant")
    linear_reg = linear_regression.linear_regression2D(xx1, xx2, y) 
    linear_reg.apply_leastsquares(order=0, test_ratio=0.1, reg_method="ols")

    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """

    print("fitting a 10th order polynomial")
    linear_reg.apply_leastsquares(order=10, test_ratio=0.1, reg_method="ols")
    
    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """
    
    MSE[ind] = linear_reg.testMSE
    print()
    print("non truncated")
    
    y = np.random.normal(loc=0.5, scale = np.sqrt(variance), size=xx1.shape)
    max_gauss[ind] = np.amax(y)
    print("mean of the output is ", np.mean(y))
    print("variance of the output is ", np.var(y))

    print("fitting a constant")
    linear_reg = linear_regression.linear_regression2D(xx1, xx2, y) 
    linear_reg.apply_leastsquares(order=0, test_ratio=0.1, reg_method="ols")


    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """
    print("fitting a 10th order polynomial")
    linear_reg.apply_leastsquares(order=10, test_ratio=0.1, reg_method="ols")
    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """
    MSE_full[ind] = linear_reg.testMSE

plt.plot(var, MSE_full, 'k', label="non truncated")
plt.plot(var, MSE, 'b', label = "truncated")
plt.xlabel("variance")
plt.ylabel("testing MSE")
plt.legend()

var = np.linspace(0.025, 1.0)
max_gauss = np.zeros(len(var))
true_var = np.zeros(len(var))
true_var_trunc = np.zeros(len(var))

for ind, variance in enumerate(var):   
    y = np.random.normal(loc=0.5, scale = np.sqrt(variance), size=10000)
    max_gauss[ind] = np.amax(y)
    true_var[ind] = np.var(y)
    #truncating
    y[np.where(y>1.0)] = 1.0
    y[np.where(y<0.0)] = 0.0
    true_var_trunc[ind] = np.var(y)
plt.plot(var, max_gauss)
plt.xlabel("variance")
plt.ylabel("Max sample")
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(var, true_var, label = "non truncated")
ax.plot(var, true_var_trunc, label = "truncated")
ax.set_xlabel("variance")
ax.set_ylabel("True variance")
ax.legend()<jupyter_output><empty_output><jupyter_text>Franke with and without truncated gaussian<jupyter_code>n = 500 #number of data points along each axis. Total data points = n*n
x1 = np.linspace(0,1,n)
x2 = np.linspace(0,1,n)
xx1, xx2 = np.meshgrid(x1, x2)
xx1 = xx1.reshape((n*n),1)
xx2 = xx2.reshape((n*n),1)

var = np.linspace(0.01, 0.1, 21)
MSE = np.zeros(len(var))
MSE_full = np.zeros(len(var))
max_gauss = np.zeros(len(var))


for ind, variance in enumerate(var):
    print()
    print("truncated")
    y = np.random.normal(loc=0.5, scale = np.sqrt(variance), size=xx1.shape)
    y[np.where(y>1.0)] = 1.0
    y[np.where(y<0.0)] = 0.0
    y+=franke.Franke(xx1, xx2, noise_var=0.0)
    
    print("mean of the output is ", np.mean(y))
    print("variance of the output is ", np.var(y))

    print("fitting a constant")
    linear_reg = linear_regression.linear_regression2D(xx1, xx2, y) 
    linear_reg.apply_leastsquares(order=0, test_ratio=0.1, reg_method="ols")

    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """

    print("fitting a 10th order polynomial")
    linear_reg.apply_leastsquares(order=10, test_ratio=0.1, reg_method="ols")
    
    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """
    
    MSE[ind] = linear_reg.testMSE
    print()
    print("non truncated")
    
    y = np.random.normal(loc=0.5, scale = np.sqrt(variance), size=xx1.shape)
    y+=franke.Franke(xx1, xx2, noise_var=0.0)
    max_gauss[ind] = np.amax(y)
    print("mean of the output is ", np.mean(y))
    print("variance of the output is ", np.var(y))

    print("fitting a constant")
    linear_reg = linear_regression.linear_regression2D(xx1, xx2, y) 
    linear_reg.apply_leastsquares(order=0, test_ratio=0.1, reg_method="ols")


    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """
    print("fitting a 10th order polynomial")
    linear_reg.apply_leastsquares(order=10, test_ratio=0.1, reg_method="ols")
    """
    print("Train MSE", linear_reg.trainMSE)
    print("Test MSE", linear_reg.testMSE)
    print("Train R2", linear_reg.trainR2)
    print("Test R2", linear_reg.testR2)
    print("Test bias", linear_reg.testbias)
    print("Test model variance", linear_reg.testvar)
    """
    MSE_full[ind] = linear_reg.testMSE

plt.plot(var, MSE_full, 'k', label="non truncated")
plt.plot(var, MSE, 'b', label = "truncated")
plt.xlabel("variance")
plt.ylabel("testing MSE")
plt.legend()<jupyter_output><empty_output><jupyter_text>generate data for widgets to display for your own set of parametersNote: This will overwrite the existent data and will take a long time to run. We do not recommend running it unless you are extremely curious about the parameters we didn't cover  Without SGD<jupyter_code># For all methods we apply the following parameters
order = np.arange(0,13)
num_points = np.array([10, 20, 30])
noise_var = np.arange(0, 6)*0.0005
test_ratio = np.round(np.arange(1,5)*0.1, 2)

n_boots = np.array([2, 5, 10], dtype=int)
k_folds = np.array([2, 5, 10], dtype=int)
ridge_lambda = np.array([0.0001, 0.001, 0.01, 0.1, 1.0])
lasso_lambda = np.array([0.000001, 0.00001, 0.0001, 0.001, 0.01])
# Save values in numpy file format
np.save("data_linear_regression/order.npy", order)
np.save("data_linear_regression/num_points.npy", num_points)
np.save("data_linear_regression/noise_var.npy", noise_var)
np.save("data_linear_regression/test_ratio.npy", test_ratio)
np.save("data_linear_regression/k_folds.npy", k_folds)
np.save("data_linear_regression/n_boots.npy", n_boots)
np.save("data_linear_regression/ridge_lambda.npy", ridge_lambda)
np.save("data_linear_regression/lasso_lambda.npy", lasso_lambda)<jupyter_output><empty_output><jupyter_text>Perform* Ordinary Least Squares(OLS) with scikit-learn test train split, own bootstrap, own cross-validation* Ridge Regression(RR) with scikit-learn test train split, own bootstrap, own cross-validation* Scikit Lasso Regression(LR) with scikit-learn test train split, own bootstrap, own cross-validation* Scikit OLS with scikit-learn test train split. Note: this is just for comparing OLS to Scikit-OLS* Scikit RR with scikit-learn test train split. Note: this is just for comparing OLS to Scikit-OLSand save output.<jupyter_code># Define for which regression methods and resampling techniques the statistical indicators should be calculated
methods = ["ols", "ols_bootstrap", "ols_crossvalidation", "ridge", 
           "ridge_bootstrap", "ridge_crossvalidation", "lasso", "lasso_bootstrap",
           "lasso_crossvalidation"]
methods = ["ols_bootstrap"]

# Loop over specified methods
for method in methods:
    if(method == "ols_crossvalidation" or method == "ridge_crossvalidation" or method == "lasso_crossvalidation"):
        test_ratio_array = np.ones(1)*0.1
    if(method != "ols_bootstrap" and method != "ridge_bootstrap" and method != "lasso_bootstrap"):
        n_boots = np.ones(1, dtype=int)
    if(method != "ols_crossvalidation" and method != "ridge_crossvalidation" and method != "lasso_crossvalidation"):
        k_folds = np.ones(1, dtype=int)
    if(method != "ridge" and method != "ridge_bootstrap" and method != "ridge_crossvalidation"):
        ridge_lambda = np.ones(1)
    if(method != "lasso" and method != "lasso_bootstrap" and method != "lasso_crossvalidation"):
        lasso_lambda = np.ones(1)
    # Calculate statistical indicators
    train_MSE, test_MSE, train_R2, test_R2, test_bias, test_var = apply_regression(order, num_points, noise_var, 
                                                                                   test_ratio, 
                                                                                   n_boots=n_boots, 
                                                                                   k_folds=k_folds, 
                                                                                   ridge_lambda=ridge_lambda, 
                                                                                   lasso_lambda=lasso_lambda, 
                                                                                   reg_type=method)
    # Save output
    np.save("data_linear_regression/train_MSE"+str(method)+".npy", train_MSE)
    np.save("data_linear_regression/test_MSE"+str(method)+".npy", test_MSE)
    np.save("data_linear_regression/train_R2"+str(method)+".npy", train_R2)
    np.save("data_linear_regression/test_R2"+str(method)+".npy", test_R2)
    np.save("data_linear_regression/test_bias"+str(method)+".npy", test_bias)
    np.save("data_linear_regression/test_variance"+str(method)+".npy", test_var)
    
    #reload variables
    test_ratio_array=np.load("data_linear_regression/test_ratio.npy")
    ridge_lambda=np.load("data_linear_regression/ridge_lambda.npy")
    lasso_lambda=np.load("data_linear_regression/lasso_lambda.npy")
    k_folds=np.load("data_linear_regression/k_folds.npy")
    n_boots=np.load("data_linear_regression/n_boots.npy")
    # To track loop progress print size of MSE output
    print(train_MSE.shape)<jupyter_output><empty_output><jupyter_text>with SGD<jupyter_code># For all methods we apply the following parameters
order = np.arange(0,13)
num_points = np.array([10, 20, 30])
noise_var = np.arange(0, 3)*0.0005
test_ratio = np.round(np.arange(1,4)*0.1, 2)

n_boots = np.array([10], dtype=int)
k_folds = np.array([10], dtype=int)
ridge_lambda = np.array([0.001, 0.01, 0.1])

# Parameters for stochastic gradient descent
learn_rates = np.array([0.001, 0.01, 0.1, 1.0])
num_min_batches = np.array([1, 2, 5, 10, 32])
epochs = np.array([10, 100, 500, 1000])
# Save values in numpy file format
np.save("data_linear_regression_sgd/order.npy", order)
np.save("data_linear_regression_sgd/num_points.npy", num_points)
np.save("data_linear_regression_sgd/noise_var.npy", noise_var)
np.save("data_linear_regression_sgd/test_ratio.npy", test_ratio)
np.save("data_linear_regression_sgd/k_folds.npy", k_folds)
np.save("data_linear_regression_sgd/n_boots.npy", n_boots)
np.save("data_linear_regression_sgd/ridge_lambda.npy", ridge_lambda)
np.save("data_linear_regression_sgd/learn_rates.npy", learn_rates)
np.save("data_linear_regression_sgd/num_min_batches.npy", num_min_batches)
np.save("data_linear_regression_sgd/epochs.npy", epochs)<jupyter_output><empty_output><jupyter_text>methods = [ "ols_sgd", "ols_bootstrap_sgd", "ols_crossvalidation_sgd", "ridge_sgd", "ridge_bootstrap_sgd"            "ridge_crossvalidation_sgd"] Loop over specified methodsfor method in methods:    if "crossvalidation" in method:        test_ratio_array = np.ones(1)*0.1    if "bootstrap" not in method:        n_boots = np.ones(1, dtype=int)    if "crossvalidation" not in method:        k_folds = np.ones(1, dtype=int)    if "ridge" not in method:        ridge_lambda = np.ones(1)    if "lasso" not in method:        lasso_lambda = np.ones(1)     Calculate statistical indicators    train_MSE, test_MSE, train_R2, test_R2, test_bias, test_var = apply_regression_sgd(order, num_points, noise_var,                                                                                    test_ratio,                                                                                    n_boots=n_boots,                                                                                    k_folds=k_folds,                                                                                    ridge_lambda=ridge_lambda,                                                                                    lasso_lambda=lasso_lambda,                                                                                    learn_rate=learn_rates,                                                                                   num_min_batch=num_min_batches,                                                                                   epochs=epochs,                                                                                   reg_type=method)     Save output    np.save("data_linear_regression_sgd/train_MSE"+str(method)+".npy", train_MSE)    np.save("data_linear_regression_sgd/test_MSE"+str(method)+".npy", test_MSE)    np.save("data_linear_regression_sgd/train_R2"+str(method)+".npy", train_R2)    np.save("data_linear_regression_sgd/test_R2"+str(method)+".npy", test_R2)    np.save("data_linear_regression_sgd/test_bias"+str(method)+".npy", test_bias)    np.save("data_linear_regression_sgd/test_variance"+str(method)+".npy", test_var)        reload variables    test_ratio=np.load("data_linear_regression_sgd/test_ratio.npy")    ridge_lambda=np.load("data_linear_regression_sgd/ridge_lambda.npy")    lasso_lambda=np.load("data_linear_regression_sgd/lasso_lambda.npy")    k_folds=np.load("data_linear_regression_sgd/k_folds.npy")    n_boots=np.load("data_linear_regression_sgd/n_boots.npy")         To track loop progress print size of MSE output    print(train_MSE.shape)<jupyter_code>pool = mp.Pool(processes=3)
methods = methods = ["e","e","d", "d", "s","d", "d", "s", "ols", "ridge_crossvalidation", "lasso", "lasso_bootstrap",
           "lasso_crossvalidation", "ols_bootstrap", "ols_crossvalidation", "ridge", 
           "ridge_bootstrap"]
#methods = methods = ["e","e","d", "d", "s","d", "d", "s", "lasso", "lasso_bootstrap",
#           "lasso_crossvalidation"]
def run_regression(method):
    print(method)
    order = np.load("data/order.npy")
    num_points = np.load("data/num_points.npy")
    noise_var = np.load("data/noise_var.npy")
    test_ratio_array = np.load("data/test_ratio_array.npy")
    k_folds = np.load("data/k_folds.npy")
    n_boots = np.load("data/n_boots.npy")
    ridge_lambda = np.load("data/ridge_lambda.npy")
    lasso_lambda = np.load("data/lasso_lambda.npy")
    
    if(method == "ols_crossvalidation" or method == "ridge_crossvalidation" or method == "lasso_crossvalidation"):
        test_ratio_array = np.ones(1)*0.1
    if(method != "ols_bootstrap" and method != "ridge_bootstrap" and method != "lasso_bootstrap"):
        n_boots = np.ones(1, dtype=int)
    if(method != "ols_crossvalidation" and method != "ridge_crossvalidation" and method != "lasso_crossvalidation"):
        k_folds = np.ones(1, dtype=int)
    if(method != "ridge" and method != "ridge_bootstrap" and method != "ridge_crossvalidation"):
        ridge_lambda = np.ones(1)
    if(method != "lasso" and method != "lasso_bootstrap" and method != "lasso_crossvalidation"):
        lasso_lambda = np.ones(1)
    # Calculate statistical indicators
    train_MSE, test_MSE, train_R2, test_R2, test_bias, test_var = apply_regression(order, num_points, noise_var, 
                                                                                   test_ratio_array, 
                                                                                   n_boots=n_boots, 
                                                                                   k_folds=k_folds, 
                                                                                   ridge_lambda=ridge_lambda, 
                                                                                   lasso_lambda=lasso_lambda, 
                                                                                   reg_type=method)
    # Save output
    test_ratio_array=np.load("data/test_ratio.npy")
    print(lasso_lambda)
    ridge_lambda=np.load("data/ridge_lambda.npy")
    lasso_lambda=np.load("data/lasso_lambda.npy")
    k_folds=np.load("data/k_folds.npy")
    n_boots=np.load("data/n_boots.npy")
    np.save("data/train_MSE"+str(method)+".npy", train_MSE)
    np.save("data/test_MSE"+str(method)+".npy", test_MSE)
    np.save("data/train_R2"+str(method)+".npy", train_R2)
    np.save("data/test_R2"+str(method)+".npy", test_R2)
    np.save("data/test_bias"+str(method)+".npy", test_bias)
    np.save("data/test_variance"+str(method)+".npy", test_var)
    # To track loop progress print size of MSE output
    print(train_MSE.shape)
    print("finished ", method)
    return 0
    
[pool.map(run_regression, methods)]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11701244813278008, 'fraction_numerical': 0.06011065006915629, 'mean_word_length': 4.14312135438896, 'pattern_counts': {'":': 1, '<': 70, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 63, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '259291', 'license': 'Apache-2.0', 'n_tokens_mistral': 13884, 'n_tokens_neox': 12479, 'n_words': 2749, 'path': 'regression_analysis/examples/linear_regression_analysis.ipynb', 'repo_name': 'am-kaiser/CompSci-Project-1', 'size': '36150'}