lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update of p1

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2025-09-01 11:15:58 +02:00
parent d30110f659
commit da97fdfc3f
2 changed files with 988 additions and 1514 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/pub/week35/ipynb/week35.ipynb
+890 -1287
View File
File diff suppressed because one or more lines are too long
doc/src/Projects/2025/Project1/Project1.do.txt
+98 -227
View File
@@ -1,4 +1,4 @@
TITLE: Project 1 on Machine Learning, deadline October 7 (midnight), 2024
TITLE: Project 1 on Machine Learning, deadline October 6 (midnight), 2025
AUTHOR: "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html" at University of Oslo, Norway
DATE: September 2
@@ -21,7 +21,7 @@ also be discussed during the various lab sessions. Please do ask us if you are i
When using codes and material from other sources, you should refer to these in the bibliography of your report, indicating wherefrom you for example
got the code, whether this is from the lecture notes, softwares like
Scikit-Learn, TensorFlow, PyTorch or other sources. These should
Scikit-Learn, TensorFlow, PyTorch or other sources such AI software. These should
always be cited correctly. How to cite some of the libraries is often
indicated from their corresponding GitHub sites or websites, see for example how to cite Scikit-Learn at URL:"https://scikit-learn.org/dev/about.html".
@@ -44,114 +44,30 @@ regression methods, including the Ordinary Least Squares (OLS) method.
In addition to the scientific part, in this course we want also to
give you an experience in writing scientific reports.
_A small recommendation when developing the codes here_. Instead of
jumping on to the two-dimensional function described below, we
recommend to do the code development and testing with a simpler
one-dimensional function, similar to those discussed in the exercises
of weeks 35 and 36. A simple test, as discussed during the lectures the first
three weeks is to set the design matrix equal to the identity
matrix. Then your model should give a mean square error which is exactly equal to zero.
When you are sure that your codes function well, you can then replace
the one-dimensional test function with the two-dimensional _Franke_ function
discussed here.
The Franke function serves as a stepping stone towards the analysis of
real topographic data. The latter is the last part of this project.
=== Description of two-dimensional function ===
We will first study how to fit polynomials to a specific
two-dimensional function called "Franke's
function":"http://www.dtic.mil/dtic/tr/fulltext/u2/a081688.pdf". This
is a function which has been widely used when testing various
interpolation and fitting algorithms. Furthermore, after having
established the model and the method, we will employ resamling
techniques such as cross-validation and/or bootstrap in order to perform a
proper assessment of our models. We will also study in detail the
so-called Bias-Variance trade off.
The Franke function, which is a weighted sum of four exponentials reads as follows
We will first study how to fit polynomials to specific
one-dimensional functions. We will start with a function given by Runge's function (see URL:"https://en.wikipedia.org/wiki/Runge%27s_phenomenon" for a discussion). The one-dimensional function we will study first is
!bt
\begin{align*}
f(x,y) &= \frac{3}{4}\exp{\left(-\frac{(9x-2)^2}{4} - \frac{(9y-2)^2}{4}\right)}+\frac{3}{4}\exp{\left(-\frac{(9x+1)^2}{49}- \frac{(9y+1)}{10}\right)} \\
&+\frac{1}{2}\exp{\left(-\frac{(9x-7)^2}{4} - \frac{(9y-3)^2}{4}\right)} -\frac{1}{5}\exp{\left(-(9x-4)^2 - (9y-7)^2\right) }.
\end{align*}
\[
f(x) = \frac{1}{1+25x^2}.
\]
!et
The function will be defined for $x,y\in [0,1]$. In a sense, our data are thus scaled to a particular domain for the input values.
Our first step will
be to perform an OLS regression analysis of this function, trying out
a polynomial fit with an $x$ and a $y$ dependence of the form $[x, y,
x^2, y^2, xy, \dots]$. We will also include bootstrap first as a
resampling technique. After that we will include the cross-validation
technique.
We can
use a uniform distribution to set up the arrays of values for $x$ and
$y$, or as in the example below just a set of fixed values for $x$ and
$y$ with a given step size. We will fit a function (for example a
polynomial) of $x$ and $y$. Thereafter we will repeat much of the
Our first step will be to perform an OLS regression analysis of this
function, trying out a polynomial fit with an $x$ dependence of the
form $[x,x^2,\dots]$. We can use a uniform distribution to set up the
arrays of values for $x \in [-5,5]$, or alternatively use a fixed step size.
Thereafter we will repeat much of the
same procedure using the Ridge and Lasso regression methods,
introducing thus a dependence on the bias (penalty) $\lambda$.
Finally we are going to use (real) digital terrain data and try to
reproduce these data using the same methods. We will also try to go
beyond the second-order polynomials metioned above and explore
which polynomial fits the data best.
We will also include bootstrap first as a
resampling technique. After that we will include the cross-validation
technique.
The Python code for the Franke function is included here (it performs also a three-dimensional plot of it)
!bc pycod
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
from random import random, seed
fig = plt.figure()
ax = fig.add_subplot(projection = '3d')
# Make data.
x = np.arange(0, 1, 0.05)
y = np.arange(0, 1, 0.05)
x, y = np.meshgrid(x,y)
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
z = FrankeFunction(x, y)
# Plot the surface.
surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
# Customize the z axis.
ax.set_zlim(-0.10, 1.40)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
!ec
If you wish to compare your results with other on the Franke function or other popular functions tested with linear regression, see the list in Figure 1 of the article by Cook et al at URL:"https://arxiv.org/abs/2401.11694".
=== Part a) : Ordinary Least Square (OLS) on the Franke function ===
We will generate our own dataset for a function
@@ -191,8 +107,8 @@ where we have defined the mean value of $\bm{y}$ as
\]
!et
Plot the resulting scores (MSE and R$^2$) as functions of the polynomial degree (here up to polymial degree five).
Plot also the parameters $\beta$ as you increase the order of the polynomial. Comment your results.
Plot the resulting scores (MSE and R$^2$) as functions of the polynomial degree (here up to polymial degree 20).
Plot also the parameters $\theta$ as you increase the order of the polynomial. Comment your results.
Your code has to include a scaling/centering of the data (for example by
subtracting the mean value), and
@@ -211,23 +127,27 @@ data and say test data. An accepted rule of thumb is to use
approximately $2/3$ to $4/5$ of the data as training data.
You can easily reuse the solutions to your exercises from week 35 and week 36.
You can easily reuse the solutions to your exercises from week 35.
See also the lecture slides from week 35 and week 36.
On scaling, we recommend reading the following section from the scikit-learn software description, see URL:"https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html#plot-all-scaling-standard-scaler-section".
=== Part b): Adding Ridge regression for the Franke function ===
=== Part b): Adding Ridge regression for the Runge function ===
Write your own code for the Ridge method, either using matrix
inversion or the singular value decomposition as done in the previous
exercise. The lecture notes from week 35 and 36 contain more information. Furthermore, the numerical exercise from week 36 is something you can reuse here.
exercise. The lecture notes from week 35 and 36 contain more information. Furthermore, the exercise from week 36 is something you can reuse here.
Perform the same analysis as you did in the previous exercise but now for different values of $\lambda$. Compare and
analyze your results with those obtained in part a) with the ordinary least squares method. Study the
dependence on $\lambda$.
Add to the
=== Part c): Adding Lasso for the Franke function ===
This exercise is essentially a repeat of the previous two ones, but now
@@ -239,65 +159,6 @@ model fits the data best.
=== Part d): Paper and pencil part ===
This exercise deals with various mean values and variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of "Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer":"https://www.springer.com/gp/book/9780387848570"). The exercise is also part of the weekly exercises for week 37.
The assumption we have made is
that there exists a continuous function $f(\bm{x})$ and a normal distributed error $\bm{\varepsilon}\sim N(0, \sigma^2)$
which describes our data
!bt
\[
\bm{y} = f(\bm{x})+\bm{\varepsilon}
\]
!et
We then approximate this function $f(\bm{x})$ with our model $\bm{\tilde{y}}$ from the solution of the linear regression equations (ordinary least squares OLS), that is our
function $f$ is approximated by $\bm{\tilde{y}}$ where we minimized $(\bm{y}-\bm{\tilde{y}})^2$, with
!bt
\[
\bm{\tilde{y}} = \bm{X}\bm{\beta}.
\]
!et
The matrix $\bm{X}$ is the so-called design or feature matrix.
Show that the expectation value of $\bm{y}$ for a given element $i$
!bt
\[
\mathbb{E}(y_i) =\sum_{j}x_{ij} \beta_j=\mathbf{X}_{i, \ast} \, \bm{\beta},
\]
!et
and that
its variance is
!bt
\[
\mbox{Var}(y_i) = \sigma^2.
\]
!et
Hence, $y_i \sim N( \mathbf{X}_{i, \ast} \, \bm{\beta}, \sigma^2)$, that is $\bm{y}$ follows a normal distribution with
mean value $\bm{X}\bm{\beta}$ and variance $\sigma^2$.
With the OLS expressions for the optimal parameters $\bm{\hat{\beta}}$ show that
!bt
\[
\mathbb{E}(\bm{\hat{\beta}}) = \bm{\beta}.
\]
!et
Show finally that the variance of $\bm{\beta}$ is
!bt
\[
\mbox{Var}(\bm{\hat{\beta}}) = \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1}.
\]
!et
We can use the last expression when we define a so-called confidence interval for the parameters $\beta$.
A given parameter $\beta_j$ is given by the diagonal matrix element of the above matrix.
=== Part e): Bias-variance trade-off and resampling techniques ===
@@ -388,80 +249,18 @@ You can follow the code example in the jupyter-book at URL:"https://compphysics.
The aim here is to implement another widely popular
resampling technique, the so-called cross-validation method.
Implement the $k$-fold cross-validation algorithm (write your own
code or use the functionality of _Scikit-Learn_) and evaluate again the MSE function resulting
Implement the $k$-fold cross-validation algorithm (feel free to use the functionality of _Scikit-Learn_ or write your own code) and evaluate again the MSE function resulting
from the test folds.
Compare the MSE you get from your cross-validation code with the one
you got from your _bootstrap_ code. Comment your results. Try $5-10$
folds.
In addition to using the ordinary least squares method, you should include both Ridge and Lasso regression.
In addition to using the ordinary least squares method, you should include both Ridge and Lasso regression in the analysis.
=== Part g): Analysis of real data ===
With our codes functioning and having been tested properly on a
simpler function we are now ready to look at real data. We will
essentially repeat in this exercise what was done in exercises a-f. However, we
need first to download the data and prepare properly the inputs to our
codes. We are going to download digital terrain data from the website
URL:"https://earthexplorer.usgs.gov/",
Or, if you prefer, we have placed selected datafiles at URL:"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Projects/2023/Project1/DataFiles"
In order to obtain data for a specific region, you need to register as
a user (free) at this website and then decide upon which area you want
to fetch the digital terrain data from. In order to be able to read
the data properly, you need to specify that the format should be _SRTM
Arc-Second Global_ and download the data as a _GeoTIF_ file. The
files are then stored in *tif* format which can be imported into a
Python program using
!bc pycod
scipy.misc.imread
!ec
Here is a simple part of a Python code which reads and plots the data
from such files
!bc pycod
import numpy as np
from imageio import imread
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
# Load the terrain
terrain1 = imread('SRTM_data_Norway_1.tif')
# Show the terrain
plt.figure()
plt.title('Terrain over Norway 1')
plt.imshow(terrain1, cmap='gray')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
!ec
If you should have problems in downloading the digital terrain data,
we provide two examples under the data folder of project 1. One is
from a region close to Stavanger in Norway and the other Møsvatn
Austfjell, again in Norway.
Feel free to produce your own terrain data.
Alternatively, if you would like to use another data set, feel free to do so. This could be data close to your reseach area or simply a data set you found interesting. See for example "kaggle.com":"https://www.kaggle.com/datasets" for examples.
Our final part deals with the parameterization of your digital terrain
data (or your own data). We will apply all three methods for linear regression, the same type (or higher order) of polynomial
approximation and cross-validation as resampling technique to evaluate which
model fits the data best.
At the end, you should present a critical evaluation of your results
and discuss the applicability of these regression methods to the type
of data presented here (either the terrain data we propose or other data sets).
@@ -548,3 +347,75 @@ encompass communities of developers in the thousands or more. And the number
of code developers and contributors keeps increasing.
"Franke's
function":"http://www.dtic.mil/dtic/tr/fulltext/u2/a081688.pdf". This
is a function which has been widely used when testing various
interpolation and fitting algorithms. Furthermore, after having
established the model and the method, we will employ resamling
techniques such as cross-validation and/or bootstrap in order to perform a
proper assessment of our models. We will also study in detail the
so-called Bias-Variance trade off.
The Franke function, which is a weighted sum of four exponentials reads as follows
!bt
\begin{align*}
f(x,y) &= \frac{3}{4}\exp{\left(-\frac{(9x-2)^2}{4} - \frac{(9y-2)^2}{4}\right)}+\frac{3}{4}\exp{\left(-\frac{(9x+1)^2}{49}- \frac{(9y+1)}{10}\right)} \\
&+\frac{1}{2}\exp{\left(-\frac{(9x-7)^2}{4} - \frac{(9y-3)^2}{4}\right)} -\frac{1}{5}\exp{\left(-(9x-4)^2 - (9y-7)^2\right) }.
\end{align*}
!et
The function will be defined for $x,y\in [0,1]$. In a sense, our data are thus scaled to a particular domain for the input values.
Finally we are going to use (real) digital terrain data and try to
reproduce these data using the same methods. We will also try to go
beyond the second-order polynomials metioned above and explore
which polynomial fits the data best.
The Python code for the Franke function is included here (it performs also a three-dimensional plot of it)
!bc pycod
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
from random import random, seed
fig = plt.figure()
ax = fig.add_subplot(projection = '3d')
# Make data.
x = np.arange(0, 1, 0.05)
y = np.arange(0, 1, 0.05)
x, y = np.meshgrid(x,y)
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
z = FrankeFunction(x, y)
# Plot the surface.
surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
# Customize the z axis.
ax.set_zlim(-0.10, 1.40)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
!ec
If you wish to compare your results with other on the Franke function or other popular functions tested with linear regression, see the list in Figure 1 of the article by Cook et al at URL:"https://arxiv.org/abs/2401.11694".