added codes
This commit is contained in:
mhjensen
@@ -0,0 +1,590 @@
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
from matplotlib.ticker import LinearLocator, FormatStrFormatter
|
||||
from matplotlib import cm
|
||||
import numpy as np
|
||||
from scipy import linalg
|
||||
import matplotlib.pyplot as plt
|
||||
import time
|
||||
|
||||
# Variance
|
||||
def var(f_model):
|
||||
n = np.size(f_model)
|
||||
f_model_mean = np.sum(f_model)/n
|
||||
#f_model_mean = np.mean(f_model)
|
||||
return np.sum((f_model-f_model_mean)**2)/n
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
# Bias
|
||||
def bias(f_true,f_model):
|
||||
n = np.size(f_model)
|
||||
#f_model_mean = np.sum(f_model)/n
|
||||
f_model_mean = np.mean(f_model)
|
||||
return np.sum((f_true-f_model_mean)**2)/n
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
# MSE
|
||||
def MSE(f_true,f_model):
|
||||
n = np.size(f_model)
|
||||
return np.sum((f_true-f_model)**2)/n
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
# Extra term
|
||||
def extra_term(f_true,f_model):
|
||||
n = np.size(f_model)
|
||||
f_model_mean = np.mean(f_model)
|
||||
return 2.0/n*np.sum((f_model_mean-f_true)*(f_model-f_model_mean))
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
# SVD invert
|
||||
def SVDinv(A):
|
||||
''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).
|
||||
SVD is numerically more stable (at least in our case) than the inversion algorithms provided by
|
||||
numpy and scipy.linalg at the cost of being slower.
|
||||
'''
|
||||
U, s, VT = linalg.svd(A)
|
||||
D = np.zeros((len(U),len(VT)))
|
||||
for i in range(0,len(VT)):
|
||||
D[i,i]=s[i]
|
||||
UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)
|
||||
return np.matmul(V,np.matmul(invD,UT))
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
# R2 score
|
||||
def R2(x_true,x_predict):
|
||||
n = np.size(x_true)
|
||||
x_avg = np.sum(x_true)/n
|
||||
enumerator = np.sum ((x_true-x_predict)**2)
|
||||
denominator = np.sum((x_true-x_avg)**2)
|
||||
return 1.0 - enumerator/denominator
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
## Mean
|
||||
#def mean(x):
|
||||
# n = np.size(x)
|
||||
# x_avg = np.sum(x)/n
|
||||
# return x_avg
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
# get sub-entries of matrix A
|
||||
def get_subset(A,indices):
|
||||
'''given an indexing set "indices", return the vector consisting of
|
||||
entries A[i,j] where (i,j) is an entry in indices.'''
|
||||
N = len(indices)
|
||||
B = np.zeros(N)
|
||||
for k in range(0,N):
|
||||
i = indices[k][0]
|
||||
j = indices[k][1]
|
||||
B[k] = A[j,i]
|
||||
return B
|
||||
|
||||
|
||||
#============================================================================================================================
|
||||
|
||||
class k_cross_validation:
|
||||
'''An k-cross validation object is initialized by passing to it data of the type linreg,
|
||||
and a paritition of the data. The class function R2 calculates the mean R2 scores
|
||||
of test and training data for the given model. The function MSE calculates the mean MSE, bias,
|
||||
variance and error terms of the test data for the given model. These quantities are stored
|
||||
as self variables.'''
|
||||
|
||||
def __init__(self, data, partition,*args):
|
||||
self.data = data; self.partition = partition; self.args = args;
|
||||
#f = data.f; X = data.X; z = data.z; correspondence = data.correspondence;
|
||||
self.k = len(partition)
|
||||
self.test_R2, self.test_var, self.test_bias, self.test_MSE, self.test_extra_terms = 0, 0, 0, 0, 0
|
||||
self.train_R2 = 0
|
||||
|
||||
#self.train_var, self.train_bias, self.train_MSE, self.train_extra_terms = 0, 0, 0, 0
|
||||
|
||||
def R2(self):
|
||||
data = self.data
|
||||
f = data.f; X = data.X; z = data.z; correspondence = data.correspondence; partition = self.partition
|
||||
k = self.k
|
||||
args = self.args
|
||||
|
||||
test_R2, train_R2 = 0, 0
|
||||
|
||||
for i, test_data in enumerate(partition):
|
||||
train_data = [x for j,x in enumerate(partition) if j!=i]
|
||||
train_data = sum(train_data, [])
|
||||
beta = data.get_beta(X[train_data],z[train_data],*args)
|
||||
freg = data.model(beta)
|
||||
test_data = [correspondence[j] for j in test_data]
|
||||
train_data = [correspondence[j] for j in train_data]
|
||||
|
||||
# test errors:
|
||||
ftest = get_subset(f,test_data); fregtest = get_subset(freg,test_data)
|
||||
test_R2 += R2(ftest,fregtest)
|
||||
|
||||
#training errors:
|
||||
ftrain = get_subset(f,train_data); fregtrain = get_subset(freg,train_data)
|
||||
train_R2 += R2(ftrain,fregtrain)
|
||||
|
||||
# self variables
|
||||
self.test_R2 = test_R2/k
|
||||
self.train_R2 = train_R2/k
|
||||
|
||||
def MSE(self):
|
||||
data = self.data
|
||||
f = data.f; X = data.X; z = data.z; correspondence = data.correspondence; partition = self.partition
|
||||
k = self.k
|
||||
args = self.args
|
||||
|
||||
test_var, test_bias, test_MSE, test_extra_terms = 0, 0, 0, 0
|
||||
#train_var, train_bias, train_MSE, train_extra_terms = 0, 0, 0, 0
|
||||
|
||||
for i, test_data in enumerate(partition):
|
||||
train_data = [x for j,x in enumerate(partition) if j!=i]
|
||||
train_data = sum(train_data, [])
|
||||
beta = data.get_beta(X[train_data],z[train_data],*args)
|
||||
freg = data.model(beta)
|
||||
test_data = [correspondence[j] for j in test_data]
|
||||
# train_data = [correspondence[j] for j in train_data]
|
||||
|
||||
# test errors:
|
||||
ftest = get_subset(f,test_data); fregtest = get_subset(freg,test_data)
|
||||
test_var += var(fregtest)
|
||||
test_bias += bias(ftest,fregtest)
|
||||
test_MSE += MSE(ftest,fregtest)
|
||||
test_extra_terms += extra_term(ftest,fregtest)
|
||||
|
||||
##training errors:
|
||||
#ftrain = get_subset(f,train_data); fregtrain = get_subset(freg,train_data)
|
||||
#train_var += var(fregtrain)
|
||||
#train_bias += bias(ftrain,fregtrain)
|
||||
#train_MSE += MSE(ftrain,fregtrain)
|
||||
#train_extra_terms += extra_term(ftrain,fregtrain)
|
||||
|
||||
# self variables
|
||||
self.test_var = test_var/k
|
||||
self.test_bias = test_bias/k
|
||||
self.test_MSE = test_MSE/k
|
||||
self.test_extra_terms = test_extra_terms/k
|
||||
|
||||
#self.train_var = train_var/k
|
||||
#self.train_bias = train_bias/k
|
||||
#self.train_MSE = train_MSE/k
|
||||
#self.train_extra_terms = train_extra_terms/k
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
class regdata:
|
||||
def __init__(self, f, degree):
|
||||
# initializing variables
|
||||
m = len(f[0,:]); n = len(f); mn = m*n;
|
||||
x = np.linspace(0, 1, m); y = np.linspace(0, 1, n); z = np.zeros(mn); xy = np.zeros((mn,2));
|
||||
|
||||
# initializing some self variables
|
||||
self.f = f; self.degree = degree; self.xm, self.ym = np.meshgrid(x,y); self.n=n;self.m=m; self.mn = mn; self.correspondence = []
|
||||
|
||||
# Making a sequence xy containing the pairs (x_i,y_j) for i,j=0,...,n, and a sequence z with matching pairs z_ij = f(x_i, y_j)
|
||||
counter = 0
|
||||
for i in range(0,m):
|
||||
for j in range(0,n):
|
||||
z[counter]=f[j,i] #wtf
|
||||
xy[counter,:] = [x[i],y[j]]
|
||||
self.correspondence.append([i,j]) #Saves the 1-1 correspondence: {counter} <-> {(i,j)} for later
|
||||
counter+=1
|
||||
self.z = z
|
||||
|
||||
# Make X
|
||||
number_basis_elts=int((degree+2)*(degree+1)/2) #(degree+1)th triangular number (number of basis elements for R[x,y] of degree <= degree)
|
||||
X = np.zeros((mn,number_basis_elts))
|
||||
powers = []
|
||||
for i in range(0,mn):
|
||||
counter = 0
|
||||
for j in range(0,degree+1):
|
||||
k = 0
|
||||
while j+k <= degree:
|
||||
xi = xy[i,0]
|
||||
yi = xy[i,1]
|
||||
X[i,counter]= (xi**j)*(yi**k)
|
||||
powers.append([j , k])
|
||||
k+=1
|
||||
counter+=1
|
||||
self.X = X
|
||||
self.powers = powers
|
||||
self.number_basis_elts = number_basis_elts
|
||||
self.invXTX = linalg.inv(np.matmul(np.transpose(X),X))
|
||||
|
||||
# Regression
|
||||
def get_reg(self, *args):
|
||||
'''Returns the polynomial fit as a numpy array. If *args is empty the fit is based on an ordinary least square.
|
||||
If *args contains a number LAMBDA, then the fit is found using Ridge for the given bias LAMBDA. If *args contains
|
||||
two numbers LAMBDA and epsilon, then the fit is found using lasso. See the function " __get_beta" for more details.'''
|
||||
|
||||
X=self.X; z=self.z #relabeling self variables
|
||||
beta = self.get_beta(X,z,*args) #obtaining beta
|
||||
reg = self.model(beta) #obtaining model from coefficients beta
|
||||
return reg
|
||||
|
||||
# Get beta (given X and z)
|
||||
def get_beta(self, X, z,*args):
|
||||
'''Returns coefficients for a given beta as a numpy array, found using either ordinary least square,
|
||||
Ridge or Lasso regression depending on the arguments. If *args is empty, then beta is found using
|
||||
ordinary least square. If *args contains a number it will be treated as a bias LAMBDA for a Ridge regression.
|
||||
If *args contains two numbers, then the first will count as a LAMBDA and the second as a tolerance epsilon.
|
||||
In this case beta is found using a shooting algorithm that runs until it converges up to the set tolerance.
|
||||
'''
|
||||
|
||||
XT = np.transpose(X)
|
||||
beta = np.matmul(XT,X)
|
||||
if len(args) >= 1: #Ridge parameter LAMBDA
|
||||
LAMBDA = args[0]
|
||||
beta[np.diag_indices_from(beta)]+=LAMBDA
|
||||
beta = SVDinv(beta)
|
||||
beta = np.matmul(beta,XT)
|
||||
beta = np.matmul(beta,z)
|
||||
|
||||
#Shooting algorithm for Lasso
|
||||
if len(args)>=2:
|
||||
epsilon = args[1]
|
||||
D = self.number_basis_elts
|
||||
ints = np.arange(0,D,1)
|
||||
beta_old = 0.0
|
||||
while np.linalg.norm(beta-beta_old)>=epsilon:
|
||||
beta_old = np.copy(beta)
|
||||
for j in range(0,D):
|
||||
aj = 2*np.sum(X[:,j]**2)
|
||||
no_j = ints[np.arange(D)!=j]
|
||||
cj = 2*np.sum(np.multiply(X[:,j],(z-np.matmul(X[:,no_j],beta[no_j]))))
|
||||
if cj<-LAMBDA:
|
||||
beta[j]=(cj+LAMBDA)/aj
|
||||
elif cj > LAMBDA:
|
||||
beta[j]=(cj-LAMBDA)/aj
|
||||
else:
|
||||
beta[j]=0.0
|
||||
return beta
|
||||
|
||||
# Get model given beta
|
||||
def model(self,beta):
|
||||
'''Returns heigh values based on the coefficients beta as a matrix
|
||||
that matches the grid xm, ym. The degree of the polynomial equals self.degree.
|
||||
'''
|
||||
xm = self.xm; ym = self.ym; degree = self.degree #relabeling self variables
|
||||
s=0
|
||||
counter = 0
|
||||
# loop that adds terms of the form beta*x^j*y^k such that j+k<=5
|
||||
for j in range(0,degree + 1):
|
||||
k = 0
|
||||
while j+k <= degree:
|
||||
s+= beta[counter]*(xm**j)*(ym**k)
|
||||
counter +=1
|
||||
k+=1
|
||||
return s
|
||||
|
||||
def get_data_partition(self,k):
|
||||
''' Creates a random partition of k (almost) equally sized parts of the array
|
||||
{1,2,...,mn}. This can be used to make training/testing data.
|
||||
'''
|
||||
mn = self.mn; correspondence = self.correspondence
|
||||
indices = np.arange(mn)
|
||||
indices_shuffle = np.arange(mn)
|
||||
np.random.shuffle(indices_shuffle)
|
||||
partition = []
|
||||
for step in range(0,k):
|
||||
part = list(indices_shuffle[step:mn:k])
|
||||
#part = [correspondence[i] for i in part]
|
||||
partition.append(part)
|
||||
return partition
|
||||
|
||||
def bootstrap_step(self, samplesize, *args):
|
||||
'''Finds and returns the coefficient that determines a model (ols, Ridge or Lasso),
|
||||
depending on args*.
|
||||
'''
|
||||
mn = self.mn; X = self.X; z = self.z; #relabeling self variables
|
||||
integers = np.random.randint(low=0, high=mn-1, size=samplesize)
|
||||
znew = z[integers]
|
||||
Xnew = X[integers,:]
|
||||
betanew = self.get_beta(Xnew,znew,*args)
|
||||
return betanew
|
||||
|
||||
# Variance/ covariance matrix
|
||||
def var_covar_matrix(self,reg):
|
||||
''' Returns the variance/covariance matrix for beta based on the given data.
|
||||
This matrix is derived from a statistical viewpoint, where one assumes beta to
|
||||
have a normal distribution.
|
||||
'''
|
||||
p = self.number_basis_elts; invXTX = self.invXTX; N = self.mn; f = self.f # Relabeling self variables
|
||||
sigma2=1.0/(N-p-1)*np.sum((f-reg)*(f-reg))
|
||||
return sigma2*invXTX # OBS! Based on matrix inversion. Inaccurate for N,p>>0.
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
def plot_3D(f,plottitle):
|
||||
''' Simple function to create 3d plot of the given data f,
|
||||
with plotitle.
|
||||
'''
|
||||
|
||||
m = len(f[0,:]); n = len(f);
|
||||
x = np.linspace(0, 1, m)
|
||||
y = np.linspace(0, 1, n);
|
||||
xm, ym = np.meshgrid(x,y)
|
||||
|
||||
# Plot f
|
||||
fig = plt.figure()
|
||||
ax = fig.gca(projection="3d")
|
||||
surf = ax.plot_surface(xm, ym, f, cmap=cm.coolwarm, linewidth=0, antialiased=False)
|
||||
|
||||
# Customize the z axis.
|
||||
ax.zaxis.set_major_locator(LinearLocator(10))
|
||||
ax.zaxis.set_major_formatter(FormatStrFormatter("%.02f"))
|
||||
ax.text2D(0.05, 0.95, plottitle, transform=ax.transAxes)
|
||||
ax.view_init(30, 60)
|
||||
|
||||
# Add a color bar which maps values to colors.
|
||||
fig.colorbar(surf, shrink=0.5, aspect=5)
|
||||
plt.show(block=False)
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
def numerical_error(data,LAMBDA):
|
||||
'''Rough numerical analysis of matrix inversions for this problem. Comparison of error and time usage
|
||||
of SVD (singular values decomposition) for matrix inversion against scipy.linalg inversion algorithm.
|
||||
Printing results to terminal.
|
||||
'''
|
||||
return_items = []
|
||||
|
||||
degree = data.degree; m = data.m; n = data.n
|
||||
# Study numerical error and time for SVD
|
||||
print("Polynomial fit of FrankeFunction in x, y of degree ", degree," with grid size ", (m,n)," analysis:")
|
||||
print("")
|
||||
X = data.X; XT = np.transpose(X); XTX = np.matmul(XT,X) #Obtaining XTX
|
||||
start_time = time.time() # start meassuring time
|
||||
inv_XTX = linalg.inv(XTX) # inversion using scipi.linalg
|
||||
end_time = time.time()
|
||||
print("Inverting XTX without SVD", "--- %s seconds ---" % (end_time - start_time)); return_items.append(end_time - start_time)
|
||||
inv_XTX_ = np.copy(inv_XTX) # storing inversion of XTX for later
|
||||
start_time = time.time()
|
||||
inv_XTX = SVDinv(XTX)
|
||||
end_time = time.time()
|
||||
print("Inverting XTX with SVD", "--- %s seconds ---" % (end_time - start_time)); return_items.append(end_time - start_time)
|
||||
print(' ')
|
||||
I_approx_ = np.matmul(inv_XTX_,XTX); # approximate I (no SVD)
|
||||
I = np.identity(len(I_approx_)); # obtaining analytical I
|
||||
output = np.linalg.norm(I_approx_-I)
|
||||
print("|(X^TX)^-1(X^TX)-I| = ",output, " (no SVD)"); return_items.append(output)
|
||||
I_approx = np.matmul(inv_XTX,XTX) # approximate I (SVD)
|
||||
output = np.linalg.norm(I_approx-I)
|
||||
print("|(X^TX)^-1(X^TX)-I| = ",np.linalg.norm(I_approx-I), " (SVD)"); return_items.append(output)
|
||||
XTX[np.diag_indices_from(XTX)]+=LAMBDA
|
||||
inv_XTX = linalg.inv(XTX)
|
||||
I_approx_ = np.matmul(inv_XTX,XTX) # approximate I (no SVD)
|
||||
output = np.linalg.norm(I_approx_-I)
|
||||
print("|(X^TX + I LAMBDA)^-1(X^TX + I LAMBDA)-I| = ",output , ", LAMBDA = ", LAMBDA, " (no SVD)"); return_items.append(output)
|
||||
inv_XTX = SVDinv(XTX)
|
||||
I_approx = np.matmul(inv_XTX,XTX)
|
||||
output = np.linalg.norm(I_approx-I)
|
||||
print("|(X^TX + I LAMBDA)^-1(X^TX + I LAMBDA)-I| = ",output, ", LAMBDA = ", LAMBDA, " (SVD)"); return_items.append(output)
|
||||
print(' ')
|
||||
|
||||
return return_items
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
def plot_R2_scores(data,Nstart,Nstop,name, epsilon = 0.001):
|
||||
''' This function makes a plot of the R2 scores vs Lambda of the different regression methods,
|
||||
for a given dataset.'''
|
||||
|
||||
degree = data.degree; f = data.f # obtaining class data
|
||||
N = Nstop-Nstart # number of lambdas
|
||||
lambdas = np.zeros(N)
|
||||
R2_ols = np.zeros(N)
|
||||
R2_Ridge = np.zeros(N)
|
||||
R2_Lasso = np.zeros(N)
|
||||
for i in range(0,N):
|
||||
LAMBDA = 10**(Nstart+i)
|
||||
lambdas[i]=LAMBDA
|
||||
R2_ols[i]=R2(f, data.get_reg())
|
||||
R2_Ridge[i]=R2(f, data.get_reg(LAMBDA))
|
||||
R2_Lasso[i]=R2(f, data.get_reg(LAMBDA,epsilon))
|
||||
print("Completed lambda: ", LAMBDA, " Completion: {:.1%}".format(float(i)/(N-1)))
|
||||
plotitle = '$R^2$ score of degree {} polynomial fit on {}'.format(degree,name)
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas),R2_ols)
|
||||
plt.plot(np.log10(lambdas),R2_Ridge)
|
||||
plt.plot(np.log10(lambdas),R2_Lasso,'--')
|
||||
plt.axis([Nstart, N+Nstart-1, 0, 1])
|
||||
plt.xlabel('log $\lambda$')
|
||||
plt.ylabel('$R^2$ score')
|
||||
plt.legend(('Ordinary least square','Ridge','Lasso'))
|
||||
plt.title(plotitle)
|
||||
plt.grid(True)
|
||||
plt.show(block=False)
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
def plot_R2_scores_k_cross_validation(data,Nstart,Nstop,k,name, epsilon = 0.001):
|
||||
''' This function makes a plot of the R2 scores vs LAMBDA of the best iteration from a k-fold cross validation on
|
||||
the data set from the given data. Best in the sense that the fit had the highest R2 score on testing data. The same
|
||||
partition of the data set is used for each lambda, and each time we select the best training data on which we base the model.
|
||||
See "k_cross_validation" for more details.'''
|
||||
|
||||
degree = data.degree; f = data.f # obtaining class data
|
||||
N = Nstop-Nstart # number of lambdas
|
||||
|
||||
# Comparing R2 scores, regression with fixed degree, variable LAMBDA
|
||||
lambdas = np.zeros(N)
|
||||
partition = data.get_data_partition(k)
|
||||
kval = k_cross_validation(data,partition)
|
||||
kval.R2()
|
||||
|
||||
R2_Lasso_test_data = np.zeros(N)
|
||||
R2_Lasso_training_data = np.zeros(N)
|
||||
R2_Ridge_test_data = np.zeros(N)
|
||||
R2_Ridge_training_data = np.zeros(N)
|
||||
|
||||
# OLS R2 score
|
||||
R2score_ols_test, R2score_ols_train = kval.test_R2, kval.train_R2
|
||||
R2_ols_test_data = np.ones(N)*R2score_ols_test
|
||||
R2_ols_training_data = np.ones(N)*R2score_ols_train
|
||||
|
||||
for i in range(0,N):
|
||||
LAMBDA = 10**(Nstart+i)
|
||||
lambdas[i]=LAMBDA
|
||||
kval = k_cross_validation(data,partition,LAMBDA)
|
||||
kval.R2()
|
||||
|
||||
# Ridge R2 score
|
||||
R2score_ridge_test, R2score_ridge_train = kval.test_R2, kval.train_R2
|
||||
R2_Ridge_test_data[i] = R2score_ridge_test
|
||||
R2_Ridge_training_data[i] = R2score_ridge_train
|
||||
|
||||
kval = k_cross_validation(data,partition,LAMBDA,epsilon)
|
||||
kval.R2()
|
||||
|
||||
# Lasso R2 score
|
||||
R2score_lasso_test, R2score_lasso_train = kval.test_R2, kval.train_R2
|
||||
R2_Lasso_test_data[i] = R2score_lasso_test
|
||||
R2_Lasso_training_data[i] = R2score_lasso_train
|
||||
|
||||
print("Completed lambda: ", LAMBDA, " Completion: {:.1%}".format(float(i)/(N-1)))
|
||||
|
||||
plotitle = '$R^2$ scores of degree {} polynomial fit on {}, $k=${}'.format(degree,name,k)
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas),R2_ols_test_data)
|
||||
plt.plot(np.log10(lambdas),R2_ols_training_data,'--')
|
||||
plt.plot(np.log10(lambdas),R2_Ridge_test_data)
|
||||
plt.plot(np.log10(lambdas),R2_Ridge_training_data,'--')
|
||||
plt.plot(np.log10(lambdas),R2_Lasso_test_data)
|
||||
plt.plot(np.log10(lambdas),R2_Lasso_training_data,'--')
|
||||
plt.axis([Nstart, Nstart+N-2, 0, 1])
|
||||
plt.xlabel('log $\lambda$')
|
||||
plt.ylabel('$R^2$ score')
|
||||
if (np.amax(R2_ols_test_data)> 0 and np.amax(R2_ols_training_data)> 0):
|
||||
plt.legend(('OLS: test data', 'OLS: training data','Ridge: test data', 'Ridge: training data','Lasso: test data', 'Lasso: training data'))
|
||||
elif (np.amax(R2_ols_test_data)<= 0 and np.amax(R2_ols_training_data)> 0):
|
||||
plt.legend(('OLS: test data (negative)', 'OLS: training data','Ridge: test data', 'Ridge: training data','Lasso: test data', 'Lasso: training data'))
|
||||
elif (np.amax(R2_ols_test_data)> 0 and np.amax(R2_ols_training_data)<= 0):
|
||||
plt.legend(('OLS: test data', 'OLS: training data (negative)','Ridge: test data', 'Ridge: training data','Lasso: test data', 'Lasso: training data'))
|
||||
elif (np.amax(R2_ols_test_data)<= 0 and np.amax(R2_ols_training_data)<= 0):
|
||||
plt.legend(('OLS: test data (negative)', 'OLS: training data (negative)','Ridge: test data', 'Ridge: training data','Lasso: test data', 'Lasso: training data'))
|
||||
plt.title(plotitle)
|
||||
plt.grid(True)
|
||||
plt.show(block=False)
|
||||
|
||||
#return ols_best, ridge_best, lasso_best
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
def plot_R2_complexity(degstart,degend,degstep,f,name, LAMBDA = 0.00001, epsilon = 0.001):
|
||||
''' Comparing R2 scores, regression with fixed LAMBDA, variable degree as well as variance and Bias
|
||||
Plotting the result.
|
||||
'''
|
||||
degrees = np.arange(degstart,degend+1,degstep)
|
||||
N = len(degrees)
|
||||
R2_ols, R2_Ridge, R2_Lasso = np.zeros(N), np.zeros(N), np.zeros(N)
|
||||
for i, degree in enumerate(degrees):
|
||||
data_f = regdata(f,degree)
|
||||
R2_ols[i]=R2(f, data_f.get_reg())
|
||||
R2_Ridge[i]=R2(f, data_f.get_reg(LAMBDA))
|
||||
R2_Lasso[i]=R2(f, data_f.get_reg(LAMBDA,epsilon))
|
||||
print("Completed degree: ", degree, " Completion: {:.1%}".format(float(i)/(N-1)))
|
||||
plotitle = '$R^2$ score of polynomial fit on {} with $\lambda=${}'.format(name,LAMBDA)
|
||||
plt.figure()
|
||||
plt.plot(degrees,R2_ols)
|
||||
plt.plot(degrees,R2_Ridge)
|
||||
plt.plot(degrees,R2_Lasso,'--')
|
||||
plt.xlabel('degree of fitting polynomial')
|
||||
plt.ylabel('$R^2$ score')
|
||||
plt.axis([degstart,degend, 0, 1])
|
||||
plt.legend(('Ordinary least square','Ridge','Lasso'))
|
||||
plt.title(plotitle)
|
||||
plt.grid(True)
|
||||
plt.show(block=False)
|
||||
|
||||
#================================================================================================================
|
||||
|
||||
def plot_MSE_variance(degstart, degend, degstep, f, LAMBDA = 0.01, epsilon = 0.001, k=10):
|
||||
# Comparing MSE, bias, variance and additional terms as function of complexity.
|
||||
degrees = np.arange(degstart,degend+1,degstep)
|
||||
N = len(degrees)
|
||||
data = regdata(f,5)
|
||||
fvar = np.zeros(N); fbias = np.zeros(N); fMSE = np.zeros(N); fextra_terms = np.zeros(N)
|
||||
|
||||
# function for plotting
|
||||
def makeplot(methodname, *args, partition = None):
|
||||
print(methodname)
|
||||
for i, degree in enumerate(degrees):
|
||||
data = regdata(f,degree)
|
||||
if partition == None:
|
||||
freg = data.get_reg(*args)
|
||||
fvar[i], fbias[i], fMSE[i], fextra_terms[i] = var(freg), bias(f,freg), MSE(f,freg), extra_term(f,freg)
|
||||
else:
|
||||
kval = k_cross_validation(data, partition, *args)
|
||||
kval.MSE()
|
||||
fvar[i] = kval.test_var
|
||||
fbias[i] = kval.test_bias
|
||||
fMSE[i] = kval.test_MSE
|
||||
fextra_terms[i] =kval.test_extra_terms
|
||||
|
||||
#fvar[i], fbias[i], fMSE[i], fextra_terms[i], train_var, train_bias, train_MSE, train_extra_terms
|
||||
print("Completed degree: ", degree, " Completion: {:.1%}".format(float(degree-degstart)/(degend-degstart)))
|
||||
plt.figure()
|
||||
plt.plot(degrees, fvar)
|
||||
plt.plot(degrees, fbias)
|
||||
plt.plot(degrees, fMSE,'--')
|
||||
plt.plot(degrees, fextra_terms)
|
||||
plt.xlabel('degree')
|
||||
plt.ylabel('Variance, bias, and MSE')
|
||||
plt.legend(('Variance','Bias','MSE','Additional term'))
|
||||
plt.grid(True)
|
||||
plt.show(block=False)
|
||||
|
||||
#It is a good idea to comment out the plots that you dont need
|
||||
|
||||
|
||||
## Ordinary least square plot
|
||||
#makeplot("Ordinary least squares")
|
||||
#plt.title("Error of ordinary least squares")
|
||||
|
||||
## Ridge plot
|
||||
#makeplot("Ridge regression",LAMBDA)
|
||||
#plt.title("Error of Ridge regression, $\lambda=${}".format(LAMBDA))
|
||||
|
||||
## Lasso plot
|
||||
#makeplot("Lasso regression",LAMBDA,epsilon)
|
||||
#plt.title("Error of lasso regression, $\lambda=${}".format(LAMBDA))
|
||||
|
||||
# k-cross validation
|
||||
partition_ = data.get_data_partition(k)
|
||||
|
||||
# Ordinary least square plot
|
||||
# makeplot("Ordinary least squares {}-fold cross validation".format(k), partition = partition_)
|
||||
# plt.title("Error OLS using {}-fold cross validation".format(k))
|
||||
|
||||
## Ridge plot
|
||||
#makeplot("Ridge regression {}-fold cross validation".format(k), LAMBDA, partition=partition_)
|
||||
#plt.title("Error Ridge using {}-fold cross validation, $\lambda=${}".format(k,LAMBDA))
|
||||
|
||||
# Lasso plot
|
||||
makeplot("Lasso regression {}-fold cross validation".format(k), LAMBDA, epsilon, partition_)
|
||||
plt.title("Error Lasso using {}-fold cross validation, $\lambda=${}".format(k,LAMBDA))
|
||||
|
||||
@@ -0,0 +1,51 @@
|
||||
# Importing functions from folder with common functions for project 1
|
||||
import sys
|
||||
sys.path.append('../functions')
|
||||
from functions import *
|
||||
from regression import OLS
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from sklearn.linear_model import Ridge as OLS_sklearn
|
||||
|
||||
|
||||
|
||||
# Making meshgrid of datapoints and compute Franke's function
|
||||
n = 3
|
||||
N = 1000
|
||||
x = np.sort(np.random.uniform(0, 1, N))
|
||||
y = np.sort(np.random.uniform(0, 1, N))
|
||||
x_mesh_, y_mesh_ = np.meshgrid(x,y)
|
||||
z = FrankeFunction(x_mesh_, y_mesh_)
|
||||
|
||||
# Add noise
|
||||
z_noise = z + np.random.normal(scale = 1, size = (N,N))
|
||||
|
||||
# Perform regression
|
||||
X = create_X(x_mesh_, y_mesh_, n=n)
|
||||
model = OLS()
|
||||
beta = model.fit(X, z_noise, ret=True)
|
||||
|
||||
# Perform regression with Scikit learn using ridge with alpha = 0
|
||||
# Because of inconsistencies in linear_regression
|
||||
model2 = OLS_sklearn(alpha = 0, fit_intercept = False)
|
||||
model2.fit(X, np.ravel(z_noise))
|
||||
|
||||
# Print beta-values of the two models
|
||||
print('============================')
|
||||
print('Calculated beta-values:', beta)
|
||||
print('Scikit-learn beta-values:', model2.coef_)
|
||||
|
||||
# Create best-fit matrix for plotting
|
||||
x_r = np.linspace(0,1,N)
|
||||
y_r = np.linspace(0,1,N)
|
||||
x_mesh, y_mesh = np.meshgrid(x,y)
|
||||
X_r = create_X(x_mesh, y_mesh, n=n)
|
||||
|
||||
# Predict
|
||||
z_reg = (model.predict(X_r)).reshape((N,N))
|
||||
plot_surface(x_mesh, y_mesh, z_reg, "OLS regression")
|
||||
print('============================ \n')
|
||||
print("MSE: %.5f" %MSE(z, z_reg))
|
||||
print("R2_Score: %.5f" %R2_Score(z, z_reg))
|
||||
|
||||
File diff suppressed because one or more lines are too long
@@ -0,0 +1,188 @@
|
||||
import sys
|
||||
import numpy as np
|
||||
from matplotlib import cm
|
||||
"""
|
||||
A file for all common functions used in project 1
|
||||
"""
|
||||
|
||||
|
||||
def FrankeFunction(x,y):
|
||||
term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
|
||||
term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
|
||||
term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
|
||||
term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
|
||||
return term1 + term2 + term3 + term4
|
||||
|
||||
def MSE(y, y_tilde):
|
||||
"""
|
||||
Function for computing mean squared error.
|
||||
Input is y: analytical solution, y_tilde: computed solution.
|
||||
"""
|
||||
return np.sum((y-y_tilde)**2)/y.size
|
||||
|
||||
def R2_Score(y, y_tilde):
|
||||
"""
|
||||
Function for computing the R2 score.
|
||||
Input is y: analytical solution, y_tilde: computed solution.
|
||||
"""
|
||||
|
||||
return 1 - np.sum((y[:-2]-y_tilde[:-2])**2)/np.sum((y[:-2]-np.average(y))**2)
|
||||
|
||||
def create_X(x, y, n = 5):
|
||||
"""
|
||||
Function for creating a X-matrix with rows [1, x, y, x^2, xy, xy^2 , etc.]
|
||||
Input is x and y mesh or raveled mesh, keyword agruments n is the degree of the polinomial you want to fit.
|
||||
"""
|
||||
if len(x.shape) > 1:
|
||||
x = np.ravel(x)
|
||||
y = np.ravel(y)
|
||||
|
||||
N = len(x)
|
||||
l = int((n+1)*(n+2)/2) # Number of elements in beta
|
||||
X = np.ones((N,l))
|
||||
|
||||
for i in range(1,n+1):
|
||||
q = int((i)*(i+1)/2)
|
||||
for k in range(i+1):
|
||||
X[:,q+k] = (x**(i-k))*(y**k)
|
||||
|
||||
return X
|
||||
|
||||
|
||||
def plot_surface(x, y, z, title = "", show = False, cmap=cm.coolwarm, figsize = None):
|
||||
"""
|
||||
Function to plot surfaces of z, given an x and y.
|
||||
Input: x, y, z (NxN'Modeler' matrices), and a title (string)
|
||||
"""
|
||||
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.ticker import LinearLocator, FormatStrFormatter
|
||||
|
||||
if figsize:
|
||||
fig = plt.figure(figsize = figsize)
|
||||
else:
|
||||
fig = plt.figure()
|
||||
|
||||
ax = fig.gca(projection='3d')
|
||||
|
||||
# Plot the surface.of the best fit
|
||||
|
||||
surf = ax.plot_surface(x, y, z, cmap=cmap,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
|
||||
# Customize the z axis automatically
|
||||
z_min = np.min(z)
|
||||
z_min = z_min*1.01 if z_min < 0 else z_min*.99
|
||||
z_max = np.max(z)
|
||||
z_max = z_max*1.01 if z_max > 0 else z_max*.99
|
||||
|
||||
ax.set_zlim(z_min, z_max)
|
||||
ax.zaxis.set_major_locator(LinearLocator(10))
|
||||
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
|
||||
ax.view_init(azim=20,elev=45)
|
||||
# Add a color bar which maps values to colors.
|
||||
fig.colorbar(surf, shrink=0.5, aspect=5)
|
||||
ax.set_title(title)
|
||||
|
||||
if show:
|
||||
plt.show()
|
||||
|
||||
return fig, ax ,surf
|
||||
|
||||
def train_test_data(x_,y_,z_,i):
|
||||
"""
|
||||
Takes in x,y and z arrays, and a array with random indesies iself.
|
||||
returns learning arrays for x, y and z with (N-len(i)) dimetions
|
||||
and test data with length (len(i))
|
||||
"""
|
||||
x_learn=np.delete(x_,i)
|
||||
y_learn=np.delete(y_,i)
|
||||
z_learn=np.delete(z_,i)
|
||||
x_test=np.take(x_,i)
|
||||
y_test=np.take(y_,i)
|
||||
z_test=np.take(z_,i)
|
||||
|
||||
return x_learn,y_learn,z_learn,x_test,y_test,z_test
|
||||
|
||||
|
||||
def K_fold(x,y,z,k,alpha,model,m=5, ret_std = False):
|
||||
"""Function to who calculate the average MSE and R2 using k-fold.
|
||||
Takes in x,y and z varibles for a dataset, k number of folds, alpha and which method beta shall use. (OLS,Ridge or Lasso)
|
||||
Returns average MSE and average R2"""
|
||||
print(m)
|
||||
if len(x.shape) > 1:
|
||||
x = np.ravel(x)
|
||||
y = np.ravel(y)
|
||||
z = np.ravel(z)
|
||||
n=len(x)
|
||||
n_k=int(n/k)
|
||||
if n_k*k!=n:
|
||||
print("k needs to be a multiple of ", n,k)
|
||||
i=np.arange(n)
|
||||
np.random.shuffle(i)
|
||||
|
||||
MSE_=0
|
||||
R2_=0
|
||||
Variance_=0
|
||||
Bias_=0
|
||||
betas = np.zeros((k,int((m+1)*(m+2)/2)))
|
||||
for t in range(k):
|
||||
x_,y_,z_,x_test,y_test,z_test=train_test_data(x,y,z,i[t*n_k:(t+1)*n_k])
|
||||
X= create_X(x_,y_,n=m)
|
||||
X_test= create_X(x_test,y_test,n=m)
|
||||
|
||||
|
||||
model.fit(X,z_)
|
||||
betas[t] = model.beta
|
||||
z_predict=model.predict(X_test)
|
||||
|
||||
MSE_+=MSE(z_test,z_predict)
|
||||
R2_+=R2_Score(z_test,z_predict)
|
||||
Bias_+=bias(z_test,z_predict)
|
||||
Variance_+=variance(z_predict)
|
||||
|
||||
return (MSE_/k, R2_/k, Bias_/k, Variance_/k, np.std(betas, axis = 0), np.mean(betas, axis = 0))
|
||||
|
||||
|
||||
def variance(y_tilde):
|
||||
"""
|
||||
Calculates the variance of the predicted values y_tilde.
|
||||
"""
|
||||
return np.sum((y_tilde - np.mean(y_tilde))**2)/np.size(y_tilde)
|
||||
|
||||
def bias(y, y_tilde):
|
||||
"""
|
||||
Calculates the bias of the predicted values y_tilde compared to
|
||||
the actual data y.
|
||||
"""
|
||||
return np.sum((y - np.mean(y_tilde))**2)/np.size(y_tilde)
|
||||
|
||||
|
||||
def update_progress(job_title, progress):
|
||||
"""
|
||||
Shows the progress of an for-loop.
|
||||
"""
|
||||
length = 20 # modify this to change the length
|
||||
block = int(round(length*progress))
|
||||
msg = "\r{0}: [{1}] {2}%".format(job_title, "#"*block + "-"*(length-block), round(progress*100, 2))
|
||||
if progress >= 1: msg += " DONE\r\n"
|
||||
sys.stdout.write(msg)
|
||||
sys.stdout.flush()
|
||||
|
||||
|
||||
def savefigure(name, figure = "gcf"):
|
||||
"""
|
||||
Function for saving figures as a .tex-file for easier integration with latex.
|
||||
"""
|
||||
try:
|
||||
from matplotlib2tikz import save as tikz_save
|
||||
tikz_save(name.replace(" ", "_") + ".tex", figure = figure, figureheight='\\figureheight', figurewidth='\\figurewidth')
|
||||
except ImportError:
|
||||
print("Please install matplotlib2tikz to save figure as a .tex-file.")
|
||||
import matplotlib.pyplot as plt
|
||||
if figure == "gcf":
|
||||
plt.savefig(name+".pdf")
|
||||
else:
|
||||
fig.savefig(name+".pdf")
|
||||
@@ -0,0 +1,81 @@
|
||||
import scipy as scipy
|
||||
import scipy.special as special
|
||||
import numpy as np
|
||||
import itertools as it
|
||||
from pandas import *
|
||||
import matplotlib.pylab as plt
|
||||
from scipy.integrate import ode
|
||||
import time
|
||||
|
||||
def unique_rows(a):
|
||||
a = np.ascontiguousarray(a)
|
||||
unique_a = np.unique(a.view([('', a.dtype)]*a.shape[1]))
|
||||
return unique_a.view(a.dtype).reshape((unique_a.shape[0], a.shape[1]))
|
||||
|
||||
def hamiltonian(n_pairs,n_basis,delta,g):
|
||||
"""
|
||||
n_pairs - Number of electron pairs
|
||||
n_basis - Number of spacial basis states
|
||||
"""
|
||||
n_SD = int(special.binom(n_basis,n_pairs))
|
||||
print("n = ", n_SD)
|
||||
H_mat = np.zeros((n_SD,n_SD))
|
||||
S = stateMatrix(n_pairs,n_basis)
|
||||
for row in range(n_SD):
|
||||
bra = S[row,:]
|
||||
for col in range(n_SD):
|
||||
ket = S[col,:]
|
||||
if np.sum(np.equal(bra,ket)) == bra.shape:
|
||||
H_mat[row,col] += 2*delta*np.sum(bra - 1) - 0.5*g*n_pairs
|
||||
if n_pairs - np.intersect1d(bra,ket).shape[0] == 1:
|
||||
H_mat[row,col] += -0.5*g
|
||||
return(H_mat)
|
||||
|
||||
|
||||
def stateMatrix(n_pairs,n_basis):
|
||||
L = []
|
||||
states = range(1,n_basis+1)
|
||||
for perm in it.permutations(states,n_pairs):
|
||||
L.append(perm)
|
||||
L = np.array(L)
|
||||
L.sort(axis=1)
|
||||
L = unique_rows(L)
|
||||
return(L)
|
||||
|
||||
|
||||
|
||||
g = 0.5
|
||||
|
||||
H = hamiltonian(4,8,1,g)
|
||||
print("Hamiltonian calculated")
|
||||
|
||||
A = H
|
||||
x = np.zeros(A.shape[0])
|
||||
x[0] = 1
|
||||
|
||||
def f(t,x):
|
||||
return(-(x.T@x)*A@x + (x.T@A@x)*x)
|
||||
|
||||
|
||||
|
||||
r = ode(f)
|
||||
r.set_initial_value(x,0) #langsos algorithm
|
||||
|
||||
t1=10
|
||||
dt=0.1
|
||||
start = time.time()
|
||||
while r.successful() and r.t < t1:
|
||||
r.integrate(r.t+dt)
|
||||
end = time.time()
|
||||
|
||||
|
||||
|
||||
print("RNN eig: ",r.y.T@A@r.y/(r.y.T@r.y))
|
||||
print("calculation time with RNN: ", end - start)
|
||||
|
||||
|
||||
start = time.time()
|
||||
eigvals, eigvecs = np.linalg.eig(H)
|
||||
end = time.time()
|
||||
print("numpy eig: ",np.sort(eigvals))
|
||||
print("calculation time with numpy eig: ", end - start)
|
||||
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user