update on today's lectures
This commit is contained in:
Morten Hjorth-Jensen
@@ -0,0 +1,44 @@
|
||||
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
from matplotlib import cm
|
||||
from matplotlib.ticker import LinearLocator, FormatStrFormatter
|
||||
import sys
|
||||
|
||||
# the number of datapoints
|
||||
n = 100
|
||||
x = 2*np.random.rand(n,1)
|
||||
y = 4+3*x*x+np.random.randn(n,1)
|
||||
|
||||
X = np.c_[np.ones((n,1)), x, x*x]
|
||||
XT_X = X.T @ X
|
||||
|
||||
#Ridge parameter lambda
|
||||
lmbda = 0.001
|
||||
Id = lmbda* np.eye(XT_X.shape[0])
|
||||
|
||||
beta_linreg = np.linalg.inv(XT_X+Id) @ X.T @ y
|
||||
print(beta_linreg)
|
||||
# Start plain gradient descent
|
||||
beta = np.random.randn(2,1)
|
||||
|
||||
eta = 0.1
|
||||
Niterations = 100
|
||||
|
||||
for iter in range(Niterations):
|
||||
gradients = 2.0/n*X.T @ (X @ (beta)-y)+2*lmbda*beta
|
||||
beta -= eta*gradients
|
||||
|
||||
print(beta)
|
||||
ypredict = X @ beta
|
||||
ypredict2 = X @ beta_linreg
|
||||
plt.plot(x, ypredict, "r-")
|
||||
plt.plot(x, ypredict2, "b-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Gradient descent example for Ridge')
|
||||
plt.show()
|
||||
@@ -0,0 +1,77 @@
|
||||
"""
|
||||
Code to test Ridge and NNs using Scikit-Learn only
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
from sklearn.neural_network import MLPRegressor
|
||||
from sklearn.metrics import accuracy_score
|
||||
import seaborn as sns
|
||||
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
# A seed just to ensure that the random numbers are the same for every run.
|
||||
# Useful for eventual debugging.
|
||||
np.random.seed(315)
|
||||
|
||||
n = 100
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
|
||||
|
||||
Maxpolydegree = 5
|
||||
X = np.zeros((n,Maxpolydegree-1))
|
||||
|
||||
for degree in range(1,Maxpolydegree): #No intercept column
|
||||
X[:,degree-1] = x**(degree)
|
||||
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
# Decide which values of lambda to use
|
||||
|
||||
nlambdas = 10
|
||||
lmbd_vals = np.logspace(-4, 0, nlambdas)
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lmbd_vals[i]
|
||||
RegRidge = linear_model.Ridge(lmb)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lmbd_vals), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
# Neural Network part
|
||||
|
||||
n_hidden_neurons = 50
|
||||
epochs = 100
|
||||
# store models for later use
|
||||
eta_vals = np.logspace(-4, 0, 10)
|
||||
# store the models for later use
|
||||
DNN_scikit = np.zeros((len(eta_vals), len(lmbd_vals)), dtype=object)
|
||||
test_accuracy = np.zeros((len(eta_vals), len(lmbd_vals)))
|
||||
sns.set()
|
||||
for i, eta in enumerate(eta_vals):
|
||||
for j, lmbd in enumerate(lmbd_vals):
|
||||
dnn = MLPRegressor(hidden_layer_sizes=(n_hidden_neurons), activation='logistic',
|
||||
alpha=lmbd, learning_rate_init=eta, max_iter=epochs)
|
||||
dnn.fit(X_train, y_train)
|
||||
ypredictMLP = dnn.predict(X_test)
|
||||
test_accuracy[i][j] = MSE(ypredictMLP, y_test)
|
||||
|
||||
fig, ax = plt.subplots(figsize = (10, 10))
|
||||
sns.heatmap(test_accuracy, annot=True, ax=ax, cmap="viridis")
|
||||
ax.set_title("Training Accuracy")
|
||||
ax.set_ylabel("$\eta$")
|
||||
ax.set_xlabel("$\lambda$")
|
||||
plt.show()
|
||||
@@ -0,0 +1,84 @@
|
||||
"""
|
||||
Code to test Ridge with own gradient descent and SGD
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
from sklearn.neural_network import MLPRegressor
|
||||
from sklearn.metrics import accuracy_score
|
||||
import seaborn as sns
|
||||
import autograd.numpy as np
|
||||
from autograd import grad
|
||||
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
# A seed just to ensure that the random numbers are the same for every run.
|
||||
# Useful for eventual debugging.
|
||||
np.random.seed(315)
|
||||
|
||||
n = 100
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
|
||||
|
||||
Maxpolydegree = 5
|
||||
X = np.zeros((n,Maxpolydegree-1))
|
||||
|
||||
for degree in range(1,Maxpolydegree): #No intercept column
|
||||
X[:,degree-1] = x**(degree)
|
||||
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
|
||||
nlambdas = 10
|
||||
lmbd_vals = np.logspace(-4, 0, nlambdas)
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lmbd_vals[i]
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
|
||||
beta = np.random.randn(X_train.shape[1],1)
|
||||
loss = np.mean((y_train.reshape(-1,1) - X_train@beta)**2)
|
||||
print(loss)
|
||||
get_grad = grad(loss,argnum=2)
|
||||
grad_beta = get_grad(X_train,y_train,beta)
|
||||
#print(grad_beta)
|
||||
|
||||
"""
|
||||
print(beta)
|
||||
print( (X_train.T @ y_train).T)
|
||||
# Make own gradient descent and define precalculated quantities, saves cycles
|
||||
XT_X = X_train.T @ X_train
|
||||
XTy = X_train.T @ y_train
|
||||
MSERidgeGDPredict = np.zeros(nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lmbd_vals[i]
|
||||
Id = lmb* np.eye(XT_X.shape[0])
|
||||
beta = np.random.randn(X_train.shape[1],1)
|
||||
eta = 0.01
|
||||
Niterations = 2
|
||||
# beta_linreg = np.linalg.pinv(XT_X+Id) @ X_train.T @ y_train
|
||||
for iter in range(Niterations):
|
||||
XX = XT_X @ beta-XTy
|
||||
gradients = (2.0/n)*XX *lmb*beta
|
||||
beta -= eta*gradients
|
||||
ypredictRidgeGD = X_test @ beta
|
||||
MSERidgeGDPredict[i] = MSE(y_test,ypredictRidgeGD)
|
||||
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lmbd_vals), MSERidgePredict, 'g--', label = 'MSE Sklearn Ridge Test')
|
||||
plt.plot(np.log10(lmbd_vals), MSERidgeGDPredict, 'r', label = 'MSE GD Ridge Test')
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
"""
|
||||
|
||||
@@ -2651,9 +2651,17 @@ Let us remind of this and recast it in terms of the mathematical operation of co
|
||||
!split
|
||||
===== Convolution Examples: Principle of Superposition and Periodic Forces (Fourier Transforms) =====
|
||||
|
||||
For problems with so-called harmonic oscillations, given by for example the following differential equation
|
||||
!bt
|
||||
\[
|
||||
m\frac{d^2x(t)}(dt^2}+\eta\frac{dx}{dt}+x(t)=F(t),
|
||||
\]
|
||||
!et
|
||||
where $F(t)$ is an applied external force acting on the system (often called a driving force), one can use the theory of Fourier transformations to find the solutions of this type of equations.
|
||||
|
||||
If one has several driving forces, $F(t)=\sum_n F_n(t)$, one can find
|
||||
the particular solution to each $F_n$, $x_{pn}(t)$, and the particular
|
||||
solution for the entire driving force is
|
||||
solution for the entire driving force is then given by a series like
|
||||
|
||||
!bt
|
||||
\begin{equation}
|
||||
@@ -2661,13 +2669,16 @@ x_p(t)=\sum_nx_{pn}(t).
|
||||
\end{equation}
|
||||
!et
|
||||
|
||||
This is known as the principal of superposition. It only applies when
|
||||
!split
|
||||
===== Principle of Superposition =====
|
||||
|
||||
This is known as the principle of superposition. It only applies when
|
||||
the homogenous equation is linear. If there were an anharmonic term
|
||||
such as $x^3$ in the homogenous equation, then when one summed various
|
||||
solutions, $x=(\sum_n x_n)^2$, one would get cross
|
||||
terms. Superposition is especially useful when $F(t)$ can be written
|
||||
as a sum of sinusoidal terms, because the solutions for each
|
||||
sinusoidal (sine or cosine) term is analytic, as we saw above.
|
||||
sinusoidal (sine or cosine) term is analytic.
|
||||
|
||||
Driving forces are often periodic, even when they are not
|
||||
sinusoidal. Periodicity implies that for some time $\tau$
|
||||
@@ -2683,6 +2694,9 @@ components in electric circuits are non-linear, e.g. diodes, which
|
||||
makes many wave forms non-sinusoidal even when the circuits are being
|
||||
driven by purely sinusoidal sources.
|
||||
|
||||
!split
|
||||
===== Simple Code Example =====
|
||||
|
||||
The code here shows a typical example of such a square wave generated using the functionality included in the _scipy_ Python package. We have used a period of $\tau=0.2$.
|
||||
|
||||
!bc pycod
|
||||
@@ -2718,6 +2732,9 @@ F(t)=\frac{f_0}{2}+\sum_{n>0} f_n\cos(2n\pi t/\tau)+g_n\sin(2n\pi t/\tau).
|
||||
\end{equation}
|
||||
!et
|
||||
|
||||
!split
|
||||
===== Wrapping up Fourier transforms =====
|
||||
|
||||
We can write down the answer for
|
||||
$x_{pn}(t)$, by substituting $f_n/m$ or $g_n/m$ for $F_0/m$. By
|
||||
writing each factor $2n\pi t/\tau$ as $n\omega t$, with $\omega\equiv
|
||||
@@ -2745,11 +2762,14 @@ x_p(t)&=&\frac{f_0}{2k}+\sum_{n>0} \alpha_n\cos(n\omega t-\delta_n)+\beta_n\sin(
|
||||
\end{eqnarray}
|
||||
!et
|
||||
|
||||
!split
|
||||
===== Finding the Coefficients =====
|
||||
|
||||
Because the forces have been applied for a long time, any non-zero
|
||||
damping eliminates the homogenous parts of the solution, so one need
|
||||
only consider the particular solution for each $n$.
|
||||
|
||||
The problem will considered solved if one can find expressions for the
|
||||
The problem is considered solved if one can find expressions for the
|
||||
coefficients $f_n$ and $g_n$, even though the solutions are expressed
|
||||
as an infinite sum. The coefficients can be extracted from the
|
||||
function $F(t)$ by
|
||||
|
||||
Reference in New Issue
Block a user