You are seeing the HTML output generated by Documenter.jl and Literate.jl from the Julia source file . The corresponding notebook can be viewed in nbviewer .
In this notebook, we apply Gaussian process regression to the Mauna Loa CO₂ dataset. It is adapted from the corresponding AbstractGPs.jl tutorial . It is therefore instructive to read that first and then see how EasyGPs.jl simplifies the steps.
We make use of the following packages:
using CSV, DataFrames # data loading
using EasyGPs # handles all things related to GPs
using Plots # visualisationLet's load and visualize the dataset.
(xtrain, ytrain), (xtest, ytest) = let
- data = CSV.read(joinpath("/home/runner/work/EasyGPs.jl/EasyGPs.jl/examples/0-mauna-loa", "CO2_data.csv"), Tables.matrix; header=0)
+ data = CSV.read(joinpath("/home/runner/work/EasyGPs.jl/EasyGPs.jl/examples/0-mauna-loa", "CO2_data.csv"), Tables.matrix; header = 0)
year = data[:, 1]
co2 = data[:, 2]
@@ -14,1487 +14,1491 @@
end
function plotdata()
- plot(; xlabel="year", ylabel="CO₂ [ppm]", legend=:bottomright)
- scatter!(xtrain, ytrain; label="training data", ms=2, markerstrokewidth=0)
- return scatter!(xtest, ytest; label="test data", ms=2, markerstrokewidth=0)
+ plot(; xlabel = "year", ylabel = "CO₂ [ppm]", legend = :bottomright)
+ scatter!(xtrain, ytrain; label = "training data", ms = 2, markerstrokewidth = 0)
+ return scatter!(xtest, ytest; label = "test data", ms = 2, markerstrokewidth = 0)
end
plotdata()
-
+
-
-
+
-
-
+
-
+We construct the GP prior using the same kernels and initial parameters as in the original tutorial.
k_smooth_trend = exp(8.0) * with_lengthscale(SEKernel(), exp(4.0))#with_lengthscale(SEKernel(), exp(4.0))
-k_seasonality = exp(2.0) * PeriodicKernel(; r=[0.5]) *
- with_lengthscale(SEKernel(), exp(4.0))
-k_medium_term_irregularities = 1.0 * with_lengthscale(RationalQuadraticKernel(; α=exp(-1.0)), 1.0)
-k_noise_terms = exp(-4.0) * with_lengthscale(SEKernel(), exp(-2.0)) + exp(-4.0) * WhiteKernel()
-kernel = k_smooth_trend + k_seasonality + k_medium_term_irregularities + k_noise_termsWe construct the GP object with the kernel above:
gp = GP(kernel)To construct a posterior, we can call the GP object with the usual AbstractGPs.jl API:
fpost_init = posterior(gp(xtrain), ytrain)Let's visualize what the GP fitted to the data looks like, for the initial choice of kernel hyperparameters.
We use the following function to plot a GP f on a specific range, using the AbstractGPs plotting recipes . By setting ribbon_scale=2 we visualize the uncertainty band with $\pm 2$ (instead of the default $\pm 1$ ) standard deviations.
plot_gp!(f; label) = plot!(f(1920:0.2:2030); ribbon_scale=2, linewidth=1, label)
+k_seasonality =
+ exp(2.0) * PeriodicKernel(; r = [0.5]) * with_lengthscale(SEKernel(), exp(4.0))
+k_medium_term_irregularities =
+ 1.0 * with_lengthscale(RationalQuadraticKernel(; α = exp(-1.0)), 1.0)
+k_noise_terms =
+ exp(-4.0) * with_lengthscale(SEKernel(), exp(-2.0)) + exp(-4.0) * WhiteKernel()
+kernel = k_smooth_trend + k_seasonality + k_medium_term_irregularities + k_noise_termsWe construct the GP object with the kernel above:
gp = GP(kernel)To construct a posterior, we can call the GP object with the usual AbstractGPs.jl API:
fpost_init = posterior(gp(xtrain), ytrain)Let's visualize what the GP fitted to the data looks like, for the initial choice of kernel hyperparameters.
We use the following function to plot a GP f on a specific range, using the AbstractGPs plotting recipes . By setting ribbon_scale=2 we visualize the uncertainty band with $\pm 2$ (instead of the default $\pm 1$ ) standard deviations.
plot_gp!(f; label) = plot!(f(1920:0.2:2030); ribbon_scale = 2, linewidth = 1, label)
plotdata()
-plot_gp!(fpost_init; label="posterior f(⋅)")
-
+
-
-
+
-
-
+
-
+A reasonable fit to the data, but poor extrapolation away from the observations!
We can now call EasyGPs.fit in order to optimize the hyperparameters. This takes care of all the parameterizations, automatic differentiation, and runs the optimizer for us. We pass an option to choose the exact same optimizer as in the original tutorial.
@time fitted_gp = EasyGPs.fit(
- gp, xtrain, ytrain;
- optimizer=Optim.LBFGS(;
- alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
- linesearch=Optim.LineSearches.BackTracking(),
- )
-)109.345328 seconds (233.97 M allocations: 32.125 GiB, 12.48% gc time, 75.27% compilation time: <1% of which was recompilation)
+ gp,
+ xtrain,
+ ytrain;
+ optimizer = Optim.LBFGS(;
+ alphaguess = Optim.LineSearches.InitialStatic(; scaled = true),
+ linesearch = Optim.LineSearches.BackTracking(),
+ ),
+)105.932129 seconds (233.94 M allocations: 32.124 GiB, 13.50% gc time, 74.74% compilation time: <1% of which was recompilation)
Let's now construct the posterior GP with the optimized hyperparameters:
fpost_opt = posterior(fitted_gp(xtrain), ytrain)This is the kernel with the point-estimated hyperparameters:
fpost_opt.prior.kernelSum of 5 kernels:
Squared Exponential Kernel (metric = Distances.Euclidean(0.0))
- Scale Transform (s = 0.006409359405350664)
@@ -1512,737 +1516,737 @@
-
+
-
-
+
-
-
+
-
+Package and system information
@@ -2251,7 +2255,7 @@ Package and system information
Status `~/work/EasyGPs.jl/EasyGPs.jl/examples/0-mauna-loa/Project.toml`
[336ed68f] CSV v0.10.13
[a93c6f00] DataFrames v1.6.1
- [dcfb08e9] EasyGPs v0.1.0 `/home/runner/work/EasyGPs.jl/EasyGPs.jl#0fce35e`
+ [dcfb08e9] EasyGPs v0.1.0 `/home/runner/work/EasyGPs.jl/EasyGPs.jl#2596518`
[98b081ad] Literate v2.16.1
[91a5bcdd] Plots v1.40.1
@@ -2277,4 +2281,4 @@ Package and system information
JULIA_DEBUG = Documenter
JULIA_LOAD_PATH = :/home/runner/.julia/packages/JuliaGPsDocs/7M86H/src
- This page was generated using Literate.jl .
