From 7f59b0091625eeb4ce6728c8de4ac58dd8a7b730 Mon Sep 17 00:00:00 2001 From: Francesco Fiusco Date: Tue, 23 Jun 2026 16:15:05 +0200 Subject: [PATCH 1/5] Fixed schedule: separated linear from non-linear regression, removed redundant environment part --- content/data-science.rst | 178 +----- content/dataformats-dataframes.rst | 19 + content/guide.rst | 12 +- content/index.rst | 3 +- content/linear-regression.rst | 578 +++++++++++++++++ ...gression.rst => non-linear-regression.rst} | 583 +----------------- 6 files changed, 626 insertions(+), 747 deletions(-) create mode 100644 content/linear-regression.rst rename content/{regression.rst => non-linear-regression.rst} (62%) diff --git a/content/data-science.rst b/content/data-science.rst index 065ec97..32907b9 100644 --- a/content/data-science.rst +++ b/content/data-science.rst @@ -20,166 +20,6 @@ Data science and machine learning - 50 min exercises -Working with data ------------------ - -In the Data Formats and Dataframes lesson, we explored a Julian approach -to manipulation and visualisation of data. - - -Here we will learn and clustering, classification, machine learning and deep learning with some toy examples. - - -Download a dataset -^^^^^^^^^^^^^^^^^^ - -We start by downloading a dataset containing measurements -of characteristic features of different penguin species. - - -.. figure:: img/lter_penguins.png - :align: center - - Artwork by @allison_horst - -.. exercise:: - - To obtain the data we simply add the PalmerPenguins package. - - .. code-block:: julia - - using Pkg - Pkg.add("PalmerPenguins") - using PalmerPenguins - - - As it was done in the Data Formats and Dataframes lesson, we can - - .. code-block:: julia - - dropmissing!(df) - -The main features we are interested in for each penguin observation are -`bill_length_mm`, `bill_depth_mm`, `flipper_length_mm` and `body_mass_g`. -What the first three features mean is illustrated in the picture below. - -.. figure:: img/culmen_depth.png - :align: center - - Artwork by @allison_horst - - -Saving the Current Setup ------------------------- - -There are several ways to save the current setup in Julia. -This section will cover three parts: saving the environment to -have reproducible code and saving data using CSV files or ``JLD``. - -1. Saving the Environment -^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. exercise:: - To check the current status of your Julia environment, you can use the status command in the package manager. - - .. code-block:: julia - - using Pkg - Pkg.status() - - .. code-block:: text - - Status `~/.julia/environments/v1.9/Project.toml` - [336ed68f] CSV v0.10.11 - [aaaa29a8] Clustering v0.15.4 - [a93c6f00] DataFrames v1.6.1 - [682c06a0] JSON v0.21.4 - [8b842266] PalmerPenguins v0.1.4 - - This will display the list of packages in the current environment along with their versions. - - To save the state of your environment, Julia uses two files: ``Project.toml`` and ``Manifest.toml``. - The ``Project.toml`` file specifies the packages that you explicitly added to your environment, - while the ``Manifest.toml`` file records the exact versions of these packages and all their dependencies1. - - When you add packages using ``Pkg.add()``, Julia automatically updates these files. - Therefore, your environment’s state (i.e., the set of loaded packages) is automatically saved. - ``Project.toml`` and ``Manifest.toml`` are located in the directory of your current Julia environment; in our case, ``~/.julia/environments/v1.9/``. - - If you want to replicate this environment on another machine or in another folder, you can do the following: - - 1. Copy both ``Project.toml`` and ``Manifest.toml`` to the new location. - 2. In Julia, navigate to that folder and activate the environment using ``Pkg.activate(".")``. - 3. Use ``Pkg.instantiate()`` to download all the necessary packages. - - More information in section `Environments` at https://enccs.github.io/julia-intro/development/ - -2. Saving Data as a CSV File -^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -As shown in the Data Formats and DataFrames lesson, a DataFrame can easily dumped into a CSV file using -the ``CSV.jl`` package, which also allows for reading tabular data. - -.. exercise:: - - You can use the CSV.jl package to save a DataFrame as a CSV file, which can be re-read later. - - .. code-block:: julia - - # using Pkg - # Pkg.add("CSV") - using CSV - CSV.write("penguins.csv", df) - - And you can load it back with: - - .. code-block:: julia - - df = CSV.read("penguins.csv", DataFrame) - -3. Saving Data Using JLD/JLD2 -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - - Another option is to use `JLD.jl `_ - The ``JLD.jl`` package provides a way to save and load Julia variables while preserving native types. - It is based on HDF5, a cross-platform, multi-language data storage format most frequently used for scientific data. - However, it is written in pure Julia and does not require any of the original C HDF5 implementation. - - The ``JLD`` package can be imported in the usual way: - - .. code-block:: julia - - using Pkg - Pkg.add("JLD") - - A DataFrame can be saved to file in the following way: - - .. code-block:: julia - - using JLD - save("penguins.jld", "df", df) - - Here we're saving ``df`` as "df" within ``penguins.jld``. You can load this DataFrame back in with: - - .. code-block:: julia - - df = load("penguins.jld", "df") - - This will return the DataFrame ``df`` from the file and assign it back to ``df``. - In the past years, the ``JLD2.jl`` package came forward as an alternative to ``JLD``. It - is also based on HDF5 and can read h5 files saved by other HDF5 implementations. It exposes an interface - similar to ``JLD`` with ``save()`` and ``load()`` functions, but the more user-friendly function ``jldsave()`` - is also available: - - .. code-block:: julia - - using JLD2 - jldsave("penguins.jld2"; df) # This is equivalent to the save command above - df = load("penguins.jld2", "df") - - Moreover, a ``jldopen()`` function provides a file-like interface. More information can be found - `here `__. - Machine learning ---------------- @@ -698,6 +538,15 @@ Whereas for Lux: true_species = OneHotArrays.onecold(ytest, ["Adelie", "Gentoo", "Chinstrap"]) ConfusionMatrix()(predicted_species, true_species) + The model can be serialised for later inference using JLD2: + + .. code-block:: julia + + @save "model_trained.jld2" model_state=Flux.state(model) + @load "model_trained.jld2" model_state + model = model(); # So the definition has to be available + Flux.loadmodel!(model, model_state); + .. group-tab:: Lux .. code-block:: julia @@ -765,6 +614,15 @@ Whereas for Lux: true_species = OneHotArrays.onecold(ytest, ["Adelie", "Gentoo", "Chinstrap"]) ConfusionMatrix()(predicted_species, true_species) + The model parameters can be saved using JLD2: + + .. code-block:: julia + + @save "trained_model.jld2" ps st + @load "trained_model.jld2" ps st + + y, st = model(x, ps, st) + Exercises --------- diff --git a/content/dataformats-dataframes.rst b/content/dataformats-dataframes.rst index 34db559..59f56c8 100644 --- a/content/dataformats-dataframes.rst +++ b/content/dataformats-dataframes.rst @@ -172,6 +172,25 @@ Here's how you can create a new dataframe: df = DataFrame(Parquet2.Dataset("penguins.parquet")) +Finally, we would also like to introduce a Julia-native format called `JLD2 +`_. This format is based on HDF5, which is +a standard format for scientific data, but does not bring external dependencies +(like the HDF5 headers). It can be used to serialise native Julia *data*, +including custom structs. While not as interoperable as Arrow, JSON or CSV, it +*can*, in principle, be opened by any HDF5 compatible client. + +Our penguins dataset can be saved and loaded as follows: + +.. code-block:: julia + + using JLD2 + @save "penguins.jld2" df + @load "penguins.jld2" df + +It is possible to save and load several variables with a dictionary-like +syntax, as well as using a file-like syntax with a do-block. For more +information, the `documentation `_ is a +good place to start. Inspect dataset ^^^^^^^^^^^^^^^ diff --git a/content/guide.rst b/content/guide.rst index 126dfca..be838d4 100644 --- a/content/guide.rst +++ b/content/guide.rst @@ -86,15 +86,15 @@ The next three days were scheduled as follows: +-------------+--------------------------------------------+ | Time | Section | +=============+============================================+ -| 9:00-9:45 | Machine learning | +| 9:00-9:45 | Clustering and classification | +-------------+--------------------------------------------+ | 9:45-10:00 | Break | +-------------+--------------------------------------------+ -| 10:00-10:45 | Clustering and Classification | +| 10:00-10:45 | Deep learning | +-------------+--------------------------------------------+ | 10:45-11:00 | Break | +-------------+--------------------------------------------+ -| 11:00-12:00 | Deep learning | +| 11:00-12:00 | Non-linear regression | +-------------+--------------------------------------------+ @@ -103,15 +103,15 @@ The next three days were scheduled as follows: +-------------+--------------------------------------------+ | Time | Section | +=============+============================================+ -| 9:00-9:45 | Linear regression | +| 9:00-9:45 | Non-linear regression (cont.) | +-------------+--------------------------------------------+ | 9:45-10:00 | Break | +-------------+--------------------------------------------+ -| 10:00-10:45 | Non-linear regression | +| 10:00-10:45 | Scientific machine learning | +-------------+--------------------------------------------+ | 10:45-11:00 | Break | +-------------+--------------------------------------------+ -| 11:00-11:45 | Scientific machine learning | +| 11:00-11:45 | Scientific machine learning (cont.) | +-------------+--------------------------------------------+ | 11:45-12:00 | Conclusions and outlook | +-------------+--------------------------------------------+ diff --git a/content/index.rst b/content/index.rst index 304f3ec..ab398bd 100644 --- a/content/index.rst +++ b/content/index.rst @@ -53,7 +53,8 @@ please visit the lesson `Julia for high-performance scientific computing |t|) Lower 95% Upper 95% + ─────────────────────────────────────────────────────────────────────── + Intercept) 3.46467 0.448322 7.73 <1e-06 2.52278 4.40656 + cX 5.05127 0.0766497 65.90 <1e-22 4.89024 5.21231 + ─────────────────────────────────────────────────────────────────────── + +.. code-block:: julia + + # note the order in the formula argument + fit(LinearModel, @formula(cX ~ cy), df) # this would model line with slope 1/5 and intercept -3.4/5 + +Now let's plot the resulting prediction (green) together with the underlying line (blue) and data points. + +.. code-block:: julia + + using Plots, GLM, DataFrames + + X = Vector(range(0, 10, length=20)) + y = 5*X .+ 3.4 + y_noisy = @. 5*X + 3.4 + randn() + + plt = plot(X, y, label="linear") + plot!(X, y_noisy, seriestype=:scatter, label="data") + + df = DataFrame(cX=X, cy=y_noisy) + lm1 = fit(LinearModel, @formula(cy ~ cX), df) + + y_pred = GLM.predict(lm1) + + # alternative: do it explicitly + # coeffs = coeftable(lm1).cols[1] # intercept and slope + # y_pred = coeffs[1] .+ coeffs[2]*X + + plot!(X, y_pred, label="predicted") + + display(plt) + + lm1 + +.. figure:: img/linear_synth_2.png + :align: center + + Image of linear model prediction. The example shown has intercept 2.9 and slope 5.1 (the result depends on random added noise). + +Multivariate linear models are done in a similar way. Now we are fitting a multivariate linear function that minimizes the sum of +squares error. In the following example we generate a linear function of 4 variables with random coefficients (normally distributed). +On top of that we add normally distributed noise. + +.. code-block:: julia + + using Plots, GLM, DataFrames + + n = 4 + C = randn(n+1,1) + X = rand(100,n) + + y = X*C[2:end] .+ C[1] + y_noisy = y .+ 0.01*randn(100,1) + + df = DataFrame(cX1=X[:,1], cX2=X[:,2], cX3=X[:,3], cX4=X[:,4], cy=y_noisy[:,1]) + + lm2 = lm(@formula(cy ~ cX1+cX2+cX3+cX4), df) + + display(lm2) + println("Coefficient vector:") + print(C) + +.. code-block:: text + + cy ~ 1 + cX1 + cX2 + cX3 + cX4 + + Coefficients: + ─────────────────────────────────────────────────────────────────────────── + Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% + ─────────────────────────────────────────────────────────────────────────── + (Intercept) -1.21114 0.00350522 -345.52 <1e-99 -1.2181 -1.20418 + cX1 2.42963 0.00375007 647.89 <1e-99 2.42218 2.43707 + cX2 -0.399002 0.00354803 -112.46 <1e-99 -0.406046 -0.391959 + cX3 -0.500017 0.00358613 -139.43 <1e-99 -0.507136 -0.492897 + cX4 1.46202 0.00365527 399.98 <1e-99 1.45476 1.46928 + ─────────────────────────────────────────────────────────────────────────── + Coefficient vector: + [-1.2045802862085417; 2.423632187920813; -0.4006938351986558; -0.5016991252146699; 1.4622712737941417;;] + +Linear models with basis functions +---------------------------------- + +Using the package GLM, we can incorporate linear models with basis functions in a convenient way, +that is to model a function as a linear combination of given non-linear functions such polynomials +or trigonometric functions. + +.. code-block:: julia + + using Plots, GLM, DataFrames + + # try this polynomial + X = range(-6, 6, length=40) + y = X.^5 .- 34*X.^3 .+ 225*X + y_noisy = y .+ randn(40,) + + plt = plot(X, y, label="polynomial") + plot!(X, y_noisy, seriestype=:scatter, label="data") + + display(plt) + +.. figure:: img/linear_basis_1.png + :align: center + + A polynomial function with noisy data. + +Fitting a polynomial to data +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Fitting a linear model with basis functions means that we try to approximate our function with for example a polynomial +:math:`p(x)=ax^5+bx^4+cx^3+dx^2+ex+f`. We fit this model to the data in a least squares sense, which works since the model +is linear in the coefficients :math:`a,b,c,d,e,f`, even though non-linear in the data :math:`x`. The degree of the polynomial needed +to get a good fit is not known in advance but for this illustration we pick the same degree (5) as when generating the data. + +.. code-block:: julia + + using Plots, GLM, DataFrames + + # try this polynomial + X = range(-6, 6, length=40) + y = X.^5 .- 34*X.^3 .+ 225*X + y_noisy = y .+ randn(40,) + + plt = plot(X, y, label="polynomial") + plot!(X, y_noisy, seriestype=:scatter, label="data") + + df = DataFrame(cX=X, cy=y_noisy) + + lm3 = lm(@formula(cy ~ cX^5 + cX^4 + cX^3 + cX^2 + cX + 1), df) + + y_pred = GLM.predict(lm3) + + plot!(X, y_pred, label="predicted") + + display(plt) + + lm3 + +.. code-block:: text + + StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}} + + cy ~ 1 + :(cX ^ 5) + :(cX ^ 4) + :(cX ^ 3) + :(cX ^ 2) + cX + + Coefficients: + ─────────────────────────────────────────────────────────────────────────────────────── + Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% + ─────────────────────────────────────────────────────────────────────────────────────── + (Intercept) -0.0354375 0.343821 -0.10 0.9185 -0.734166 0.663291 + cX ^ 5 1.00118 0.000551333 1815.92 <1e-85 1.00006 1.0023 + cX ^ 4 -0.000992084 0.00169158 -0.59 0.5614 -0.00442979 0.00244563 + cX ^ 3 -34.054 0.0236797 -1438.11 <1e-82 -34.1021 -34.0058 + cX ^ 2 0.0230557 0.0571179 0.40 0.6890 -0.0930219 0.139133 + cX 225.511 0.226822 994.22 <1e-76 225.05 225.972 + ─────────────────────────────────────────────────────────────────────────────────────── + +.. figure:: img/linear_basis_1_pred.png + :align: center + + Fitting a polynomial to data. + +Exercises +--------- + +Let us illustrate linear regression on real data sets. The first dataset comes from the RDatasets package +and are data from chemical experiments for the production of formaldehyde. +The data columns are amount of Carbohydrate (ml) and Optical Density of a purple color on a spectrophotometer. + +Sources: + +- Bennett, N. A. and N. L. Franklin (1954), Statistical Analysis in Chemistry and the Chemical Industry, New York: Wiley. +- McNeil, D. R. (1977), Interactive Data Analysis, New York: Wiley. + +.. exercise:: + + In the exerises below we use the packages GLM, RDatasets, Plots and DataFrames: + + .. code-block:: julia + + using Pkg + Pkg.add("GLM") + Pkg.add("RDatasets") + Pkg.add("Plots") + Pkg.add("DataFrames") + +.. exercise:: Formaldehyde example + + To load the dataset, you can do: + + .. code-block:: julia + + using GLM, RDatasets, Plots + df = dataset("datasets", "Formaldehyde") + + The columns of the dataframe are called `Carb` and `OptDen` for the amount of Carbohydrate and Optical Density. + You can plot the data as follows: + + .. code-block:: julia + + plt = plot(df.Carb, df.OptDen, seriestype=:scatter, label="formaldehyde data") + display(plt) + + To model Density as a linear function of Carbohydrate you can do as follows. + The `predict` method is used to make model predictions. + + .. code-block:: julia + + model = fit(LinearModel, @formula(OptDen ~ Carb), df) + y_pred = GLM.predict(model) + + To add the prediction to the plot and print the model results you can do: + + .. code-block:: julia + + plot!(df.Carb, y_pred, label="model") + display(plt) + model + + .. solution:: A suggestion + + .. code-block:: julia + + using GLM, RDatasets, Plots + + df = dataset("datasets", "Formaldehyde") + + plt = plot(df.Carb, df.OptDen, seriestype=:scatter, label="formaldehyde data") + + display(plt) + + model = fit(LinearModel, @formula(OptDen ~ Carb), df) + + y_pred = GLM.predict(model) + + plot!(df.Carb, y_pred, label="model") + + display(plt) + + model + + .. figure:: img/linear_formaldehyde.png + :align: center + +.. exercise:: Changing hyperparameters + + Take a look at the code in the example `Fitting a polynomial to data`_. + This fit is pretty tight. + + - What happens if you increase the noise by say 100 times? + - What happens if if you use a degree 6 or 7 polynomial to fit the data instead? + + You can try the second experiment with the original noise level. + + .. solution:: + + You can change the following rows: + + .. code-block:: julia + + # y_noisy = y .+ randn(40,) + y_noisy = y .+ 100*randn(40,) + + # lm3 = lm(@formula(cy ~ cX^5 + cX^4 + cX^3 + cX^2 + cX + 1), df) + lm3 = lm(@formula(cy ~ cX^7 + cX^6 + cX^5 + cX^4 + cX^3 + cX^2 + cX + 1), df) + +Let us have a look at linear regression on real multidimensional data. For this we will use the Rdatasets +package and the "trees" dataset, which consists of measurements on +black cherry trees: girth, height and volume +(see Atkinson, A. C. (1985) Plots, Transformations and Regression. Oxford University Press). + +.. exercise:: Black cherry trees + + In this exercise we use also the package StatsBase: + + .. code-block:: julia + + using Pkg + Pkg.add("StatsBase") + + Load the trees data set as follows: + + .. code-block:: julia + + using GLM, RDatasets, StatsBase, Plots + # Girth Height and Volume of Black Cherry Trees + trees = dataset("datasets", "trees") + df = trees + + Randomly split the data set into a training and testing data set. + + .. code-block:: julia + + n_rows = size(df)[1] + rows_train = sample(1:n_rows, Int(round(n_rows*0.8)), replace=false) + rows_test = [x for x in 1:n_rows if ~(x in rows_train)] + + L_train = df[rows_train,:] + L_test = df[rows_test,:] + + It is reasonable to try to fit the logarithm of volume as a linear function of + the logarithm of the height and logarithm of the girth. This is because the + volume is presumably roughly proportional to the height times the girth squared. + + .. code-block:: julia + + # reasonable to look at logarithms since we can expect something like V~h*g^2 and + # log V = constant + log h + 2log g + model = fit(LinearModel, @formula(log(Volume) ~ log(Girth) + log(Height)), L_train) + + Lastly, make predictions on the training set according to the model and compute the + root mean squared error of the prediction (for instance on the training set). + + .. code-block:: julia + + Z = L_train + # Z = L_test + y_pred = GLM.predict(model, Z) + + # Root Mean Squared Error + rmse = sqrt(sum((exp.(y_pred) - Z.Volume).^2)/size(Z)[1]) + + .. solution:: The whole script + + .. code-block:: julia + + using GLM, RDatasets, StatsBase, Plots + # Girth Height and Volume of Black Cherry Trees + trees = dataset("datasets", "trees") + df = trees + + n_rows = size(df)[1] + rows_train = sample(1:n_rows, Int(round(n_rows*0.8)), replace=false) + rows_test = [x for x in 1:n_rows if ~(x in rows_train)] + + L_train = df[rows_train,:] + L_test = df[rows_test,:] + + # reasonable to look at logarithms since can expect something like V~h*r^2 and + # log V = constant + log h + 2log r + model = fit(LinearModel, @formula(log(Volume) ~ log(Girth) + log(Height)), L_train) + + Z = L_train + # Z = L_test + y_pred = GLM.predict(model, Z) + + # Root Mean Squared Error + rmse = sqrt(sum((exp.(y_pred) - Z.Volume).^2)/size(Z)[1]) + + println(rmse) + df + + .. code-block:: julia-repl + + 2.2631848027992776 # rmse + + 31×3 DataFrame + Row │ Girth Height Volume + │ Float64 Int64 Float64 + ─────┼────────────────────────── + 1 │ 8.3 70 10.3 + 2 │ 8.6 65 10.3 + 3 │ 8.8 63 10.2 + 4 │ 10.5 72 16.4 + 5 │ 10.7 81 18.8 + 6 │ 10.8 83 19.7 + 7 │ 11.0 66 15.6 + 8 │ 11.0 75 18.2 + 9 │ 11.1 80 22.6 + 10 │ 11.2 75 19.9 + 11 │ 11.3 79 24.2 + + And so on (31 data points). + + +.. exercise:: Trigonometric basis functions + + Try a similar example as the polynomial above but with trigonometric functions :math:`y(x)=\cos(x)+\cos(2x)`. + Here is a snippet that generates data for this example: + + .. code-block:: julia + + using Plots, GLM, DataFrames + + X = range(-6, 6, length=100) + y = cos.(X) .+ cos.(2*X) + y_noisy = y .+ 0.1*randn(100,) + + To make a dataframe out of the data and fit a linear model to it, you can do: + + .. code-block:: julia + + df = DataFrame(X=X, y=y_noisy) + lm1 = lm(@formula(y ~ 1 + cos(X) + cos(2*X) + cos(3*X) + cos(4*X)), df) + + .. solution:: A suggestion. + + .. code-block:: julia + + using Plots, GLM, DataFrames + + # try a cosine combination + X = range(-6, 6, length=100) + y = cos.(X) .+ cos.(2*X) + y_noisy = y .+ 0.1*randn(100,) + + plt = plot(X, y, label="waveform") + plot!(X, y_noisy, seriestype=:scatter, label="data") + + display(plt) + + df = DataFrame(X=X, y=y_noisy) + + lm1 = lm(@formula(y ~ 1 + cos(X) + cos(2*X) + cos(3*X) + cos(4*X)), df) + + .. code-block:: text + + StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}} + + y ~ 1 + :(cos(X)) + :(cos(2X)) + :(cos(3X)) + :(cos(4X)) + + Coefficients: + ──────────────────────────────────────────────────────────────────────────── + Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% + ──────────────────────────────────────────────────────────────────────────── + (Intercept) 0.0130408 0.0108222 1.21 0.2312 -0.00844393 0.0345256 + cos(X) 0.981561 0.015653 62.71 <1e-78 0.950486 1.01264 + cos(2X) 0.984984 0.0156219 63.05 <1e-78 0.953971 1.016 + cos(3X) -0.0135547 0.015573 -0.87 0.3863 -0.044471 0.0173616 + cos(4X) 0.0148532 0.0155105 0.96 0.3407 -0.015939 0.0456454 + ──────────────────────────────────────────────────────────────────────────── + + .. figure:: img/linear_basis_2.png + :align: center + + Fitting trigonometric functions to data. + +Loading data +------------ + +We will now have a look at a climate data set containing daily mean +temperature, humidity, wind speed and mean pressure at a location in +Delhi India over a period of several years. The data set is available +`here `__. +In the context of the Delhi dataset we have borrowed some elements of Sebastian Callh's personal +blog post *Forecasting the weather with neural ODEs* found `here +`__. + +.. code-block:: julia + + using DataFrames, CSV, Plots, Statistics + + # data_path = + # a string, full path to data file DailyDelhiClimateTrain.csv + # uploaded in julia-for-hpda/content/data + df_train = CSV.read(data_path, DataFrame) + df_train + + M = [df_train.meantemp df_train.humidity df_train.wind_speed df_train.meanpressure] + plottitles = ["meantemp" "humidity" "wind_speed" "meanpressure"] + plotylabels = ["C°" "g/m^3?" "km/h?" "hPa"] + # color=[1 2 3 4] gives default colors + plot(M, layout=(4,1), color=[1 2 3 4], legend=false, title=plottitles, + xlabel="time (days)", ylabel=plotylabels, size=(800,800)) + +.. figure:: img/climate_plots_first.png + :align: center + + Plots of measurements. + +The mean pressure data field seems to contain some unreasonably large values. Let us filter those out and consider these missing data. + +.. code-block:: julia + + using DataFrames, CSV, Plots, Statistics + + # data_path = + # a string, full path to data file DailyDelhiClimateTrain.csv + # uploaded in julia-for-hpda/content/data + df_train = CSV.read(data_path, DataFrame) + + M = [df_train.meantemp df_train.humidity df_train.wind_speed df_train.meanpressure] + + plottitles = ["meantemp" "humidity" "wind_speed" "meanpressure"] + plotylabels = ["C°" "g/m^3?" "km/h?" "hPa"] + + df_train[df_train.meanpressure .< 950,:meanpressure] .= NaN + df_train[1050 .< df_train.meanpressure,:meanpressure] .= NaN + + M = [df_train.meantemp df_train.humidity df_train.wind_speed df_train.meanpressure] + + # color=[1 2 3 4] gives default colors + plt = plot(M, layout=(4,1), color=[1 2 3 4], legend=false, title=plottitles, + xlabel="time (days)", ylabel=plotylabels, size=(800,800)) + + display(plt) + +.. figure:: img/climate_plots_second.png + :align: center + + Plots of cleaned up data. + diff --git a/content/regression.rst b/content/non-linear-regression.rst similarity index 62% rename from content/regression.rst rename to content/non-linear-regression.rst index 5572f2d..986fc15 100644 --- a/content/regression.rst +++ b/content/non-linear-regression.rst @@ -1,583 +1,5 @@ -.. _regression: - -Regression and time-series prediction -===================================== - -.. questions:: - - - How can I perform simple linear regression in Julia? - - How to do linear regression with non-linear basis functions? - - How do to basic Fourier based regression? - - How to perform non-linear regression and time-series prediction? - -.. instructor-note:: - - - 90 min teaching - - 60 min exercises - -.. callout:: - - The code in this lesson is written for Julia v1.12.6. - -Linear regression with synthetic data -------------------------------------- - -We begin with some simple examples of linear regression on generated data. -For the models we will use the package GLM (Generalized Linear Models), -which among other things contains linear regression models. - -Let's start by generating some data along a line and add normally distributed noise. - -.. code-block:: julia - - using Plots, GLM, DataFrames - - X = Vector(range(0, 10, length=20)) - y = 5*X .+ 3.4 - y_noisy = @. 5*X + 3.4 + randn() - - plt = plot(X, y, label="linear") - plot!(X, y_noisy, seriestype=:scatter, label="data") - - display(plt) - -.. figure:: img/linear_synth_1.png - :align: center - -Given data :math:`x_1,x_2,\ldots,x_k` and responses :math:`y_1,y_2,\ldots,y_k`, the ordinary least squares method -finds the linear function :math:`l(x) = ax+b` minimizing the sum of squares error :math:`\sum_i (l(x_i)-y_i)^2`. - -.. code-block:: julia - - using Plots, GLM, DataFrames - - X = Vector(range(0, 10, length=20)) - y = 5*X .+ 3.4 - y_noisy = @. 5*X + 3.4 + randn() - - df = DataFrame(cX=X, cy=y_noisy) - lm1 = fit(LinearModel, @formula(cy ~ cX), df) - - # the above is the same as @formula(cy ~ cX + 1), which also works - - # alternative syntax - # lm(@formula(cy ~ cX), df) - -.. code-block:: text - - StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}} - - cy ~ 1 + cX # the constant term (intercept) is there, same as if we do @formula(cy ~ cX + 1) - - Coefficients: - ─────────────────────────────────────────────────────────────────────── - Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% - ─────────────────────────────────────────────────────────────────────── - Intercept) 3.46467 0.448322 7.73 <1e-06 2.52278 4.40656 - cX 5.05127 0.0766497 65.90 <1e-22 4.89024 5.21231 - ─────────────────────────────────────────────────────────────────────── - -.. code-block:: julia - - # note the order in the formula argument - fit(LinearModel, @formula(cX ~ cy), df) # this would model line with slope 1/5 and intercept -3.4/5 - -Now let's plot the resulting prediction (green) together with the underlying line (blue) and data points. - -.. code-block:: julia - - using Plots, GLM, DataFrames - - X = Vector(range(0, 10, length=20)) - y = 5*X .+ 3.4 - y_noisy = @. 5*X + 3.4 + randn() - - plt = plot(X, y, label="linear") - plot!(X, y_noisy, seriestype=:scatter, label="data") - - df = DataFrame(cX=X, cy=y_noisy) - lm1 = fit(LinearModel, @formula(cy ~ cX), df) - - y_pred = GLM.predict(lm1) - - # alternative: do it explicitly - # coeffs = coeftable(lm1).cols[1] # intercept and slope - # y_pred = coeffs[1] .+ coeffs[2]*X - - plot!(X, y_pred, label="predicted") - - display(plt) - - lm1 - -.. figure:: img/linear_synth_2.png - :align: center - - Image of linear model prediction. The example shown has intercept 2.9 and slope 5.1 (the result depends on random added noise). - -Multivariate linear models are done in a similar way. Now we are fitting a multivariate linear function that minimizes the sum of -squares error. In the following example we generate a linear function of 4 variables with random coefficients (normally distributed). -On top of that we add normally distributed noise. - -.. code-block:: julia - - using Plots, GLM, DataFrames - - n = 4 - C = randn(n+1,1) - X = rand(100,n) - - y = X*C[2:end] .+ C[1] - y_noisy = y .+ 0.01*randn(100,1) - - df = DataFrame(cX1=X[:,1], cX2=X[:,2], cX3=X[:,3], cX4=X[:,4], cy=y_noisy[:,1]) - - lm2 = lm(@formula(cy ~ cX1+cX2+cX3+cX4), df) - - display(lm2) - println("Coefficient vector:") - print(C) - -.. code-block:: text - - cy ~ 1 + cX1 + cX2 + cX3 + cX4 - - Coefficients: - ─────────────────────────────────────────────────────────────────────────── - Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% - ─────────────────────────────────────────────────────────────────────────── - (Intercept) -1.21114 0.00350522 -345.52 <1e-99 -1.2181 -1.20418 - cX1 2.42963 0.00375007 647.89 <1e-99 2.42218 2.43707 - cX2 -0.399002 0.00354803 -112.46 <1e-99 -0.406046 -0.391959 - cX3 -0.500017 0.00358613 -139.43 <1e-99 -0.507136 -0.492897 - cX4 1.46202 0.00365527 399.98 <1e-99 1.45476 1.46928 - ─────────────────────────────────────────────────────────────────────────── - Coefficient vector: - [-1.2045802862085417; 2.423632187920813; -0.4006938351986558; -0.5016991252146699; 1.4622712737941417;;] - -Linear models with basis functions ----------------------------------- - -Using the package GLM, we can incorporate linear models with basis functions in a convenient way, -that is to model a function as a linear combination of given non-linear functions such polynomials -or trigonometric functions. - -.. code-block:: julia - - using Plots, GLM, DataFrames - - # try this polynomial - X = range(-6, 6, length=40) - y = X.^5 .- 34*X.^3 .+ 225*X - y_noisy = y .+ randn(40,) - - plt = plot(X, y, label="polynomial") - plot!(X, y_noisy, seriestype=:scatter, label="data") - - display(plt) - -.. figure:: img/linear_basis_1.png - :align: center - - A polynomial function with noisy data. - -Fitting a polynomial to data -^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Fitting a linear model with basis functions means that we try to approximate our function with for example a polynomial -:math:`p(x)=ax^5+bx^4+cx^3+dx^2+ex+f`. We fit this model to the data in a least squares sense, which works since the model -is linear in the coefficients :math:`a,b,c,d,e,f`, even though non-linear in the data :math:`x`. The degree of the polynomial needed -to get a good fit is not known in advance but for this illustration we pick the same degree (5) as when generating the data. - -.. code-block:: julia - - using Plots, GLM, DataFrames - - # try this polynomial - X = range(-6, 6, length=40) - y = X.^5 .- 34*X.^3 .+ 225*X - y_noisy = y .+ randn(40,) - - plt = plot(X, y, label="polynomial") - plot!(X, y_noisy, seriestype=:scatter, label="data") - - df = DataFrame(cX=X, cy=y_noisy) - - lm3 = lm(@formula(cy ~ cX^5 + cX^4 + cX^3 + cX^2 + cX + 1), df) - - y_pred = GLM.predict(lm3) - - plot!(X, y_pred, label="predicted") - - display(plt) - - lm3 - -.. code-block:: text - - StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}} - - cy ~ 1 + :(cX ^ 5) + :(cX ^ 4) + :(cX ^ 3) + :(cX ^ 2) + cX - - Coefficients: - ─────────────────────────────────────────────────────────────────────────────────────── - Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% - ─────────────────────────────────────────────────────────────────────────────────────── - (Intercept) -0.0354375 0.343821 -0.10 0.9185 -0.734166 0.663291 - cX ^ 5 1.00118 0.000551333 1815.92 <1e-85 1.00006 1.0023 - cX ^ 4 -0.000992084 0.00169158 -0.59 0.5614 -0.00442979 0.00244563 - cX ^ 3 -34.054 0.0236797 -1438.11 <1e-82 -34.1021 -34.0058 - cX ^ 2 0.0230557 0.0571179 0.40 0.6890 -0.0930219 0.139133 - cX 225.511 0.226822 994.22 <1e-76 225.05 225.972 - ─────────────────────────────────────────────────────────────────────────────────────── - -.. figure:: img/linear_basis_1_pred.png - :align: center - - Fitting a polynomial to data. - -Exercises ---------- - -Let us illustrate linear regression on real data sets. The first dataset comes from the RDatasets package -and are data from chemical experiments for the production of formaldehyde. -The data columns are amount of Carbohydrate (ml) and Optical Density of a purple color on a spectrophotometer. - -Sources: - -- Bennett, N. A. and N. L. Franklin (1954), Statistical Analysis in Chemistry and the Chemical Industry, New York: Wiley. -- McNeil, D. R. (1977), Interactive Data Analysis, New York: Wiley. - -.. exercise:: - - In the exerises below we use the packages GLM, RDatasets, Plots and DataFrames: - - .. code-block:: julia - - using Pkg - Pkg.add("GLM") - Pkg.add("RDatasets") - Pkg.add("Plots") - Pkg.add("DataFrames") - -.. exercise:: Formaldehyde example - - To load the dataset, you can do: - - .. code-block:: julia - - using GLM, RDatasets, Plots - df = dataset("datasets", "Formaldehyde") - - The columns of the dataframe are called `Carb` and `OptDen` for the amount of Carbohydrate and Optical Density. - You can plot the data as follows: - - .. code-block:: julia - - plt = plot(df.Carb, df.OptDen, seriestype=:scatter, label="formaldehyde data") - display(plt) - - To model Density as a linear function of Carbohydrate you can do as follows. - The `predict` method is used to make model predictions. - - .. code-block:: julia - - model = fit(LinearModel, @formula(OptDen ~ Carb), df) - y_pred = GLM.predict(model) - - To add the prediction to the plot and print the model results you can do: - - .. code-block:: julia - - plot!(df.Carb, y_pred, label="model") - display(plt) - model - - .. solution:: A suggestion - - .. code-block:: julia - - using GLM, RDatasets, Plots - - df = dataset("datasets", "Formaldehyde") - - plt = plot(df.Carb, df.OptDen, seriestype=:scatter, label="formaldehyde data") - - display(plt) - - model = fit(LinearModel, @formula(OptDen ~ Carb), df) - - y_pred = GLM.predict(model) - - plot!(df.Carb, y_pred, label="model") - - display(plt) - - model - - .. figure:: img/linear_formaldehyde.png - :align: center - -.. exercise:: Changing hyperparameters - - Take a look at the code in the example `Fitting a polynomial to data`_. - This fit is pretty tight. - - - What happens if you increase the noise by say 100 times? - - What happens if if you use a degree 6 or 7 polynomial to fit the data instead? - - You can try the second experiment with the original noise level. - - .. solution:: - - You can change the following rows: - - .. code-block:: julia - - # y_noisy = y .+ randn(40,) - y_noisy = y .+ 100*randn(40,) - - # lm3 = lm(@formula(cy ~ cX^5 + cX^4 + cX^3 + cX^2 + cX + 1), df) - lm3 = lm(@formula(cy ~ cX^7 + cX^6 + cX^5 + cX^4 + cX^3 + cX^2 + cX + 1), df) - -Let us have a look at linear regression on real multidimensional data. For this we will use the Rdatasets -package and the "trees" dataset, which consists of measurements on -black cherry trees: girth, height and volume -(see Atkinson, A. C. (1985) Plots, Transformations and Regression. Oxford University Press). - -.. exercise:: Black cherry trees - - In this exercise we use also the package StatsBase: - - .. code-block:: julia - - using Pkg - Pkg.add("StatsBase") - - Load the trees data set as follows: - - .. code-block:: julia - - using GLM, RDatasets, StatsBase, Plots - # Girth Height and Volume of Black Cherry Trees - trees = dataset("datasets", "trees") - df = trees - - Randomly split the data set into a training and testing data set. - - .. code-block:: julia - - n_rows = size(df)[1] - rows_train = sample(1:n_rows, Int(round(n_rows*0.8)), replace=false) - rows_test = [x for x in 1:n_rows if ~(x in rows_train)] - - L_train = df[rows_train,:] - L_test = df[rows_test,:] - - It is reasonable to try to fit the logarithm of volume as a linear function of - the logarithm of the height and logarithm of the girth. This is because the - volume is presumably roughly proportional to the height times the girth squared. - - .. code-block:: julia - - # reasonable to look at logarithms since we can expect something like V~h*g^2 and - # log V = constant + log h + 2log g - model = fit(LinearModel, @formula(log(Volume) ~ log(Girth) + log(Height)), L_train) - - Lastly, make predictions on the training set according to the model and compute the - root mean squared error of the prediction (for instance on the training set). - - .. code-block:: julia - - Z = L_train - # Z = L_test - y_pred = GLM.predict(model, Z) - - # Root Mean Squared Error - rmse = sqrt(sum((exp.(y_pred) - Z.Volume).^2)/size(Z)[1]) - - .. solution:: The whole script - - .. code-block:: julia - - using GLM, RDatasets, StatsBase, Plots - # Girth Height and Volume of Black Cherry Trees - trees = dataset("datasets", "trees") - df = trees - - n_rows = size(df)[1] - rows_train = sample(1:n_rows, Int(round(n_rows*0.8)), replace=false) - rows_test = [x for x in 1:n_rows if ~(x in rows_train)] - - L_train = df[rows_train,:] - L_test = df[rows_test,:] - - # reasonable to look at logarithms since can expect something like V~h*r^2 and - # log V = constant + log h + 2log r - model = fit(LinearModel, @formula(log(Volume) ~ log(Girth) + log(Height)), L_train) - - Z = L_train - # Z = L_test - y_pred = GLM.predict(model, Z) - - # Root Mean Squared Error - rmse = sqrt(sum((exp.(y_pred) - Z.Volume).^2)/size(Z)[1]) - - println(rmse) - df - - .. code-block:: text - - 2.2631848027992776 # rmse - - 31×3 DataFrame - Row │ Girth Height Volume - │ Float64 Int64 Float64 - ─────┼────────────────────────── - 1 │ 8.3 70 10.3 - 2 │ 8.6 65 10.3 - 3 │ 8.8 63 10.2 - 4 │ 10.5 72 16.4 - 5 │ 10.7 81 18.8 - 6 │ 10.8 83 19.7 - 7 │ 11.0 66 15.6 - 8 │ 11.0 75 18.2 - 9 │ 11.1 80 22.6 - 10 │ 11.2 75 19.9 - 11 │ 11.3 79 24.2 - - And so on (31 data points). - - -.. exercise:: Trigonometric basis functions - - Try a similar example as the polynomial above but with trigonometric functions :math:`y(x)=\cos(x)+\cos(2x)`. - Here is a snippet that generates data for this example: - - .. code-block:: julia - - using Plots, GLM, DataFrames - - X = range(-6, 6, length=100) - y = cos.(X) .+ cos.(2*X) - y_noisy = y .+ 0.1*randn(100,) - - To make a dataframe out of the data and fit a linear model to it, you can do: - - .. code-block:: julia - - df = DataFrame(X=X, y=y_noisy) - lm1 = lm(@formula(y ~ 1 + cos(X) + cos(2*X) + cos(3*X) + cos(4*X)), df) - - .. solution:: A suggestion. - - .. code-block:: julia - - using Plots, GLM, DataFrames - - # try a cosine combination - X = range(-6, 6, length=100) - y = cos.(X) .+ cos.(2*X) - y_noisy = y .+ 0.1*randn(100,) - - plt = plot(X, y, label="waveform") - plot!(X, y_noisy, seriestype=:scatter, label="data") - - display(plt) - - df = DataFrame(X=X, y=y_noisy) - - lm1 = lm(@formula(y ~ 1 + cos(X) + cos(2*X) + cos(3*X) + cos(4*X)), df) - - .. code-block:: text - - StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}} - - y ~ 1 + :(cos(X)) + :(cos(2X)) + :(cos(3X)) + :(cos(4X)) - - Coefficients: - ──────────────────────────────────────────────────────────────────────────── - Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95% - ──────────────────────────────────────────────────────────────────────────── - (Intercept) 0.0130408 0.0108222 1.21 0.2312 -0.00844393 0.0345256 - cos(X) 0.981561 0.015653 62.71 <1e-78 0.950486 1.01264 - cos(2X) 0.984984 0.0156219 63.05 <1e-78 0.953971 1.016 - cos(3X) -0.0135547 0.015573 -0.87 0.3863 -0.044471 0.0173616 - cos(4X) 0.0148532 0.0155105 0.96 0.3407 -0.015939 0.0456454 - ──────────────────────────────────────────────────────────────────────────── - - .. figure:: img/linear_basis_2.png - :align: center - - Fitting trigonometric functions to data. - -Loading data ------------- - -We will now have a look at a climate data set containing daily mean -temperature, humidity, wind speed and mean pressure at a location in -Delhi India over a period of several years. The data set is available -`here `__. -In the context of the Delhi dataset we have borrowed some elements of Sebastian Callh's personal -blog post *Forecasting the weather with neural ODEs* found `here -`__. - -.. code-block:: julia - - using DataFrames, CSV, Plots, Statistics - - # data_path = - # a string, full path to data file DailyDelhiClimateTrain.csv - # uploaded in julia-for-hpda/content/data - df_train = CSV.read(data_path, DataFrame) - df_train - - M = [df_train.meantemp df_train.humidity df_train.wind_speed df_train.meanpressure] - plottitles = ["meantemp" "humidity" "wind_speed" "meanpressure"] - plotylabels = ["C°" "g/m^3?" "km/h?" "hPa"] - # color=[1 2 3 4] gives default colors - plot(M, layout=(4,1), color=[1 2 3 4], legend=false, title=plottitles, - xlabel="time (days)", ylabel=plotylabels, size=(800,800)) - -.. figure:: img/climate_plots_first.png - :align: center - - Plots of measurements. - -The mean pressure data field seems to contain some unreasonably large values. Let us filter those out and consider these missing data. - -.. code-block:: julia - - using DataFrames, CSV, Plots, Statistics - - # data_path = - # a string, full path to data file DailyDelhiClimateTrain.csv - # uploaded in julia-for-hpda/content/data - df_train = CSV.read(data_path, DataFrame) - - M = [df_train.meantemp df_train.humidity df_train.wind_speed df_train.meanpressure] - - plottitles = ["meantemp" "humidity" "wind_speed" "meanpressure"] - plotylabels = ["C°" "g/m^3?" "km/h?" "hPa"] - - df_train[df_train.meanpressure .< 950,:meanpressure] .= NaN - df_train[1050 .< df_train.meanpressure,:meanpressure] .= NaN - - M = [df_train.meantemp df_train.humidity df_train.wind_speed df_train.meanpressure] - - # color=[1 2 3 4] gives default colors - plt = plot(M, layout=(4,1), color=[1 2 3 4], legend=false, title=plottitles, - xlabel="time (days)", ylabel=plotylabels, size=(800,800)) - - display(plt) - -.. figure:: img/climate_plots_second.png - :align: center - - Plots of cleaned up data. - -Non-linear regression ---------------------- +Non-linear regression and time-series prediction +------------------------------------------------ In this section we will have a look at non-linear regression methods. @@ -1489,3 +911,4 @@ To decrease overfitting, we may project to a lower dimensional subspace of basis :align: center Three models of varying crudeness and overfit. + From 7c473547d0906715a17fe84b9bec8f3cc5418159 Mon Sep 17 00:00:00 2001 From: Francesco Fiusco Date: Tue, 23 Jun 2026 16:37:55 +0200 Subject: [PATCH 2/5] Removed overlapping linear algebra stuff with julia-intro --- content/linear-algebra.rst | 143 +------------------------------------ 1 file changed, 1 insertion(+), 142 deletions(-) diff --git a/content/linear-algebra.rst b/content/linear-algebra.rst index 886cba6..04b8996 100644 --- a/content/linear-algebra.rst +++ b/content/linear-algebra.rst @@ -5,10 +5,8 @@ Linear algebra .. questions:: - - How can I create vectors and matrices in Julia? - - How can I perform vector and matrix operations in Julia? - - How to generate random matrices and perform sparse matrix computations? - How does Principle Component Analysis work? + - How to generate random matrices and perform sparse matrix computations? .. instructor-note:: @@ -19,145 +17,6 @@ Linear algebra The code in this lesson is written for Julia v1.12.6. -Vectors and matrices in Julia ------------------------------ - -We will start with a brief look at how we can create arrays -and vectors in Julia and how to perform vector and matrix operations. - -.. code-block:: julia - - # lazy range notation, list from 1 to 10 - 1:10 - - # make into vector - Vector(1:10) - - # another way to make ranges - range(1, 10) - -.. code-block:: julia-repl - - julia> Vector(1:10) - 10-element Vector{Int64}: - 1 - 2 - 3 - 4 - ... - 8 - 9 - 10 - -Indexing elements or parts of vectors and matrices can be done with slicing as in Python or Matlab. - -.. code-block:: julia - - # form vector and matrix - u = [2,3,5,7] - A = [1 2 3;4 5 6;7 8 9] - - # extract elements from vector - u[1] # first element: 2 - u[2] # second element: 3 - u[2:4] # range second to fourth: 3,5,7 - - # slicing - A[2,3] # second row third column: 6 - A[:,1] # first column: 1,4,7 - A[2,:] # second row: 4,5,6 - - # zeros - zeros(5) # [0,0,0,0,0] - zeros(5,5) # 5x5-matrix of zeros - - # ones - ones(5) # [1,1,1,1,1] - ones(5,5) # 5x5-matrix of ones - -.. code-block:: julia-repl - - julia> u - 4-element Vector{Int64}: - 2 - 3 - 5 - 7 - - julia> A - 3×3 Matrix{Int64}: - 1 2 3 - 4 5 6 - 7 8 9 - - julia> zeros(5,5) - 5×5 Matrix{Float64}: - 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 - - julia> ones(5,5) - 5×5 Matrix{Float64}: - 1.0 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 1.0 - -To perform vector and matrix operations we can use a syntax similar to Matlab or Python. - -.. code-block:: julia - - # forming vectors - a = [1,2,3,4] - b = [2,3,4,5] - - # scaling - 0.5*a - - # vector addition - a + b - a - b - - # powers - a^2 # MethodError - a.^2 # 1,4,9,16 - - # same as vector addition - a .+ b - - # element wise product - a.*b - - # applying functions - sin(a) # MethodError - sin.(a) # element wise computations - - # alternative way - @. a+a^2-sin(a)*sin(b) - - # forming matrix and vector - A = [1 2 3;4 5 6;7 8 9] - v = [1,2,3] - - # vector matrix multiplication - A*v - - # matrix multiplication - B = A*A - - # Matrix multiplication - A*B - - # matrix powers - A^3 - - # transpose - transpose(A) - A' - Eigenvectors and eigenvalues ---------------------------- From 34a5f5be9f106e109efa45fe6f19b94f7964cf11 Mon Sep 17 00:00:00 2001 From: Francesco Fiusco Date: Tue, 23 Jun 2026 16:48:10 +0200 Subject: [PATCH 3/5] Used OneHotArrays for both Flux and Lux for consistency --- content/data-science.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/content/data-science.rst b/content/data-science.rst index 32907b9..acc8d03 100644 --- a/content/data-science.rst +++ b/content/data-science.rst @@ -483,7 +483,7 @@ Whereas for Lux: # onecold (opposite to onehot) gives back the original representation function accuracy(x, y) - return sum(Flux.onecold(model(x)) .== Flux.onecold(y)) / size(y, 2) + return sum(OneHotArrays.onecold(model(x)) .== OneHotArrays.onecold(y)) / size(y, 2) end ``model`` will be our neural network, so we go ahead and define it: @@ -724,8 +724,8 @@ Exercises accuracy(xtest, ytest) # 0.9850746268656716 - predicted_species = Flux.onecold(model(xtest), ["Adelie", "Gentoo", "Chinstrap"]) - true_species = Flux.onecold(ytest, ["Adelie", "Gentoo", "Chinstrap"]) + predicted_species = OneHotArrays.onecold(model(xtest), ["Adelie", "Gentoo", "Chinstrap"]) + true_species = OneHotArrays.onecold(ytest, ["Adelie", "Gentoo", "Chinstrap"]) ConfusionMatrix()(predicted_species, true_species) .. group-tab:: Lux From 36076bcc1a31bc34ef9cbf3f8c1276456f12b601 Mon Sep 17 00:00:00 2001 From: ffrancesco94 Date: Wed, 24 Jun 2026 10:19:28 +0200 Subject: [PATCH 4/5] Implemented feedback for data science part --- content/data-science.rst | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/content/data-science.rst b/content/data-science.rst index acc8d03..4774ad7 100644 --- a/content/data-science.rst +++ b/content/data-science.rst @@ -512,7 +512,6 @@ Whereas for Lux: ytrain[:,1:5] # accuracy before training accuracy(xtrain, ytrain) - accuracy(xtest, ytest) Finally we are ready to train the model. Let's run 100 epochs: @@ -527,9 +526,6 @@ Whereas for Lux: accuracy(xtrain, ytrain) accuracy(xtest, ytest) - The performance of the model is probably somewhat underwhelming, but you will - fix that in an exercise below! - We finally create a confusion matrix to quantify the performance of the model: .. code-block:: julia @@ -586,7 +582,6 @@ Whereas for Lux: .. code-block:: julia accuracy(model, ps, st, xtrain, ytrain) - accuracy(model, ps, st, xtest, ytest) Now we can start the training loop! @@ -623,6 +618,16 @@ Whereas for Lux: y, st = model(x, ps, st) +As you can see, the performance of the model is less than optimal. This is due +to a combination of two effects: +- The features have very different ranges (from about 15 to 6000) +- We are using a sigmoid activation function + +The sigmoid function is capped between 0 and 1 and tends to "squish" the very +large and very small values. When the features have different ranges, sigmoid +leads to a loss of information, which prevents the network from learning +effectively. We will try to fix this, as well as apply other optimisations, in +the exercises below. Exercises --------- @@ -630,12 +635,9 @@ Exercises .. exercise:: Improve the deep learning model Improve the performance of the neural network we trained above! - The network is not improving much because of the large numerical - range of the input features (from around 15 to around 6000) combined - with the fact that we use a ``sigmoid`` activation function. A standard - method in machine learning is to normalize features by "batch - normalization". Replace the network definition with the following and - see if the performance improves: + A standard method in machine learning is to normalize features by "batch + normalization". Replace the network definition with the following and see if + the performance improves: .. tabs:: From b7ada1d01e2640ab519b93f6a9b76e6ae589f559 Mon Sep 17 00:00:00 2001 From: ffrancesco94 Date: Wed, 24 Jun 2026 10:28:11 +0200 Subject: [PATCH 5/5] Fixed imports for mean, std in data science --- content/data-science.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/content/data-science.rst b/content/data-science.rst index 4774ad7..82d0041 100644 --- a/content/data-science.rst +++ b/content/data-science.rst @@ -435,7 +435,7 @@ Whereas for Lux: .. code-block:: julia - using MLJ: partition, ConfusionMatrix + using MLJ: partition, ConfusionMatrix, mean, std using DataFrames using PalmerPenguins using OneHotArrays