02582

# To embed plots in the notebooks %matplotlib inline import matplotlib.pyplot as plt import numpy as np # numpy library import scipy.linalg as lng # linear algebra from scipy library from sklearn import preprocessing as preproc # load preprocessing function # seaborn can be used to "prettify" default matplotlib plots by importing and setting as default import seaborn as sns sns.set() # Set searborn as default

prostatePath = 'Prostate.txt' T = np.loadtxt(prostatePath, delimiter = ' ', skiprows = 1, usecols=[1,2,3,4,5,6,7,8,9]) y = T[:, 0] X = T[:,1:] [n, p] = np.shape(X) #Our chosen normalization centers and normalize the variables of a data matrix to unit length. # We can use sklearn "Normalizer" to do this, but we must transpose the matrices to act on the variables instead of samples X_pre = X - np.mean(X,axis=0) y_pre = y - np.mean(y,axis=0) normalizer = preproc.Normalizer().fit(X_pre.T) X_pre = normalizer.transform(X_pre.T).T

def ridgeMulti(X, _lambda, p, y): inner_prod = np.linalg.inv(np.matmul(X.T, X) + _lambda * np.eye(p,p)) outer_prod = np.matmul(X.T, y) betas = np.matmul(inner_prod, outer_prod) return betas

k = 100; # try k values of lambda lambdas = np.logspace(-4, 4, k) betas = np.zeros((p,k)) for i in range(k): betas[:,i] = ridgeMulti(X,lambdas[i],p,y) plt.figure() plt.semilogx(lambdas, betas.T ) plt.xlabel("Lambdas") plt.ylabel("Betas") plt.suptitle("Regularized beta estimates", fontsize=20) plt.show()

def centerData(data): mu = np.mean(data,axis=0) data = data - mu return data, mu

K = 10 #lambdas = np.logspace(-4, 4, k) #Create a vector of length n that contains equal amounts of numbers from 1 to K x = list(range(1, K)) array = ((np.repeat(x, int(n/K), axis=0))) pp = np.random.permutation(array) p1 = pp == 1 print(pp[p1])

# Calculate the standard error for the best lambda, and find a the largest lambda with a MSE that is within # the range of the optimal lambda +- the standard error. Lambda_CV_1StdRule = print("CV lambda with 1-std-rule %0.2f" % Lambda_CV_1StdRule)

D = np.zeros(100) AIC = np.zeros(100) BIC = np.zeros(100) # calculate the AIC and BIC for ridge regression for different lambdas # pick the lambda based on themodels with the lowest AIC and BIC jAIC = jBIC = print("AIC at %0.2f" % lambdas[jAIC]) print("BIC at %0.2f" % lambdas[jBIC])

#plot different methods Information criteria: AIC BIC CV fig, axs = plt.subplots(1, 2, figsize=(15,5)) axs[0].set_title('Information criteria') _ = axs[0].semilogx(lambdas,AIC,'-r',label='AIC') axs[0].semilogx(lambdas,BIC/300,'-b',label='BIC') axs[0].semilogx(lambdas,meanMSE,'-k',label='CV') axs[0].semilogx(lambdas[jAIC],np.min(AIC),'*r',label='AIC') axs[0].semilogx(lambdas[jBIC],np.min(BIC)/300,'*b',label='BIC') axs[0].semilogx(lambdas[jOpt],np.min(meanMSE),'*k',label='CV') axs[0].legend() axs[0].set_xlabel('Lambda') axs[0].set_ylabel('Criteria') #plot the degree of freedom axs[1].semilogx(lambdas,D,'-k'); axs[1].set_title('Degrees of freedom'); axs[1].set_xlabel('Lambda'); axs[1].set_ylabel('d');

NBoot = 100 Beta = np.zeros((p, len(lambdas), NBoot)) # Run bootstrap on ridge for the same range of lambdas as before # similarly as with cross validation you run the bootstrap x times and each time you sample with replacement # ridge is then run on that sample for each lambda value # Plot the varriance of the betas for each lambda value plt.figure() for i in range(8): plt.semilogx(lambdas, stdBeta[i,:]) plt.title("Bootstrapped variance") plt.ylabel("beta") plt.xlabel("lambda") plt.show()