Gaussian processes

Project.toml

@@ -7,12 +7,14 @@ authors = [
 version = "0.13.0"
 
 [deps]
+AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 Bootstrap = "e28b5b4c-05e8-5b66-bc03-6f0c0a0a06e0"
 CovarianceEstimation = "587fd27a-f159-11e8-2dae-1979310e6154"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
@@ -23,6 +25,7 @@ MeshAdaptiveDirectSearch = "f4d74008-4565-11e9-04bd-4fe404e6a92a"
 Monomials = "272bfe72-f66c-432f-a94d-600f29493792"
 Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+ParameterHandling = "2412ca09-6db7-441c-8e3a-88d5709968c5"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
@@ -31,13 +34,16 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+AbstractGPs = "0.5.24"
 Bootstrap = "2.2"
 CovarianceEstimation = "0.2"
 DataFrames = "0.22, 1.0"
 DelimitedFiles = "1"
 Dierckx = "0.5"
+DifferentiationInterface = "0.7.7"
 Distributions = "0.24, 0.25"
 FastGaussQuadrature = "0.4, 0.5, 1"
 FiniteDifferences = "0.12"
@@ -46,11 +52,13 @@ MeshAdaptiveDirectSearch = "0.1.0"
 Monomials = "1.0"
 Mustache = "1.0"
 Optim = "1.9.4"
+ParameterHandling = "0.5.0"
 Primes = "0.5"
 QuadGK = "2.11.1"
 QuasiMonteCarlo = "0.3"
 Reexport = "0.2, 1.0"
 Roots = "2.2.2"
 Statistics = "1"
 StatsBase = "0.33, 0.34"
+Zygote = "0.7.10"
 julia = "1.10"

demo/metamodels/gaussianprocess.jl

@@ -0,0 +1,43 @@
+using UncertaintyQuantification
+
+x = RandomVariable.(Uniform(-5, 5), [:x1, :x2])
+
+himmelblau = Model(
+    df -> (df.x1 .^ 2 .+ df.x2 .- 11) .^ 2 .+ (df.x1 .+ df.x2 .^ 2 .- 7) .^ 2, :y
+)
+
+design = LatinHypercubeSampling(80)
+
+mean_f = ConstMean(0.0)
+kernel = SqExponentialKernel() ∘ ARDTransform([1.0, 1.0])
+σ² = 1e-5
+
+gp_prior = with_gaussian_noise(GP(mean_f, kernel), σ²)
+
+using Optim
+
+optimizer = MaximumLikelihoodEstimation(Optim.Adam(alpha=0.005), Optim.Options(; iterations=10, show_trace=false))
+
+input_transform = ZScoreTransform()
+
+gp_model = GaussianProcess(
+    gp_prior,
+    x,
+    himmelblau,
+    :y,
+    design;
+    input_transform=input_transform,
+    optimization=optimizer
+)
+
+test_data = sample(x, 1000)
+evaluate!(gp_model, test_data; mode=:mean_and_var)
+
+test_data = sample(x, 1000)
+evaluate!(gp_model, test_data)
+evaluate!(himmelblau, test_data)
+
+mse = mean((test_data.y .- test_data.y_mean) .^ 2)
+println("MSE is:  $mse")
+
+# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

docs/Project.toml

@@ -5,10 +5,12 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 UncertaintyQuantification = "7183a548-a887-11e9-15ce-a56ab60bad7a"
 
 [compat]
 Documenter = "1.14.1"
 DocumenterCitations = "1.4.1"
 DocumenterVitepress = "0.2.6"
+Optim = "1.13.2"

docs/literate/metamodels/gaussianprocess.jl

@@ -0,0 +1,128 @@
+#===
+# Gaussian Process Regression 
+
+## Himmelblau's Function
+
+In this example, we will model the following test function (known as Himmelblau's function) in the range ``x1, x2 ∈ [-5, 5]`` with a Gaussian process (GP) regression model.
+
+It is defined as:
+
+ ```math
+f(x1, x2) = (x1^2 + x2 - 11)^2 + (x1 + x2^2 - 7)^2.
+```
+===#
+# ![](himmelblau.svg)
+#===
+Aanalogue to the response surface example, we create an array of random variables, that will be used when evaluating the points that our experimental design produces.
+===#
+
+using UncertaintyQuantification
+
+x = RandomVariable.(Uniform(-5, 5), [:x1, :x2])
+
+himmelblau = Model(
+    df -> (df.x1 .^ 2 .+ df.x2 .- 11) .^ 2 .+ (df.x1 .+ df.x2 .^ 2 .- 7) .^ 2, :y
+)
+
+#===
+Next, we chose a experimental design. In this example, we are using a `LatinHyperCube` design from which we draw 80 samples to train our model:
+===#
+
+design = LatinHypercubeSampling(80)
+
+#===
+After that, we construct a prior GP model. Here we assume a constant mean of 0.0 and a squared exponential kernel with automatic relevance determination (ARD). 
+We also assume a small Gaussian noise term in the observations for numerical stability:
+===#
+
+mean_f = ConstMean(0.0)
+kernel = SqExponentialKernel() ∘ ARDTransform([1.0, 1.0])
+σ² = 1e-5
+
+gp_prior = with_gaussian_noise(GP(mean_f, kernel), σ²)
+
+#===
+Next, we set up an optimizer used in the log marginal likelihood maximization to find the optimal hyperparameters of our GP model. Here we use the Adam optimizer from the `Optim.jl` package with a learning rate of 0.005 and run it for 10 iterations.:
+===#
+using Optim
+
+optimizer = MaximumLikelihoodEstimation(Optim.Adam(alpha=0.005), Optim.Options(; iterations=10, show_trace=false))
+
+#===
+Finally, we define an input standardization (here a z-score transform). While not strictly necessary for this example, standardization can help finding good hyperparameters. 
+Note that we can also define an output transform to scale the output for training the GP. When evaluating the GP model, the input will be automatically transformed with the fitted standardization.
+The output will be transformed back to the original scale automatically as well.
+===#
+
+input_transform = ZScoreTransform()
+
+#===
+The GP regression model is now constructed by calling the `GaussianProcess` constructor with the prior GP, the input random variables, the model, the output symbol, the experimental design, and the optional input and output transforms as well as the hyperparameter optimization method.
+The construction then samples the experimental design, evaluates the model at the sampled points, standardizes the input and output data, optimizes the hyperparameters of the GP, and constructs the posterior GP.
+===#
+#md using Random #hide
+#md Random.seed!(42) #hide
+
+gp_model = GaussianProcess(
+    gp_prior, 
+    x, 
+    himmelblau, 
+    :y, 
+    design; 
+    input_transform=input_transform, 
+    optimization=optimizer
+)
+
+#===
+To evaluate the `GaussianProcess`, use `evaluate!(gp::GaussianProcess, data::DataFrame)` with the `DataFrame` containing the points you want to evaluate. 
+The evaluation of a GP is not unique, and we can choose to evaluate the mean prediction, the prediction variance, a combination of both, or draw samples from the posterior distribution.
+The default is to evaluate the mean prediction.
+We can specify the evaluation mode via the `mode` keyword argument. Supported options are:
+- `:mean` - predictive mean (default)
+- `:var` - predictive variance
+- `:mean_and_var` - both mean and variance
+- `:sample` - random samples from the predictive distribution
+===#
+
+test_data = sample(x, 1000)
+evaluate!(gp_model, test_data; mode=:mean_and_var)
+
+#===
+The mean prediction of our model in this case has an mse of about 65 and looks like this in comparison to the original:
+===#
+
+#md using Plots #hide
+#md using DataFrames #hide
+#md a = range(-5, 5; length=200) #hide
+#md b = range(-5, 5; length=200) #hide
+#md A = repeat(collect(a)', length(b), 1) #hide
+#md B = repeat(collect(b), 1, length(a)) #hide
+#md df = DataFrame(x1 = vec(A), x2 = vec(B)) #hide
+#md evaluate!(gp_model, df; mode=:mean_and_var) #hide
+#md evaluate!(himmelblau, df) #hide
+#md gp_mean = reshape(df[:, :y_mean], length(b), length(a)) #hide
+#md gp_var = reshape(df[:, :y_var], length(b), length(a)) #hide
+#md himmelblau_values = reshape(df[:, :y], length(b), length(a)) #hide
+#md s1 = surface(a, b, himmelblau_values; plot_title="Himmelblau's function")
+#md s2 = surface(a, b, gp_mean; plot_title="GP posterior mean")
+#md plot(s1, s2, layout = (1, 2), legend = false)
+#md savefig("gp-mean-comparison.svg") # hide
+#md s3 = surface(a, b, gp_var; plot_title="GP posterior variance") # hide
+#md plot(s3, legend = false) #hide
+#md savefig("gp-variance.svg"); nothing # hide
+
+# ![](gp-mean-comparison.svg)
+
+#===
+Note that the mse in comparison to the response surface model (with an mse of about 1e-26) is significantly higher.
+However, the GP model also provides a measure of uncertainty in its predictions via the predictive variance.
+===#
+
+# ![](gp-variance.svg)
+
+#jl test_data = sample(x, 1000)
+#jl evaluate!(gp_model, test_data)
+#jl evaluate!(himmelblau, test_data)
+
+#jl mse = mean((test_data.y .- test_data.y_mean) .^ 2)
+#jl println("MSE is:  $mse")

docs/make.jl

@@ -59,6 +59,7 @@ makedocs(;
 			"Reliability" => "api/reliability.md",
 			"ResponseSurface" => "api/responsesurface.md",
 			"PolyharmonicSpline" => "api/polyharmonicspline.md",
+			"Gaussian Processes" => "api/gaussianprocesses.md",
 			"Simulations" => "api/simulations.md",
 			"Bayesian Updating" => "api/bayesianupdating.md",
 			"Power Spectral Density Functions" => "api/psd.md",

docs/references.bib

@@ -400,6 +400,20 @@ @book{raiffaAppliedStatisticalDecision1961
   pagetotal = {356}
 }
 
+@book{rasmussen2005gaussian,
+	title = {Gaussian {Processes} for {Machine} {Learning}},
+	copyright = {http://creativecommons.org/licenses/by-nc-nd/4.0/},
+	isbn = {978-0-262-25683-4},
+	url = {https://direct.mit.edu/books/book/2320/Gaussian-Processes-for-Machine-Learning},
+	language = {en},
+	urldate = {2025-09-04},
+	publisher = {The MIT Press},
+	author = {Rasmussen, Carl Edward and Williams, Christopher K. I.},
+	month = nov,
+	year = {2005},
+	doi = {10.7551/mitpress/3206.001.0001},
+}
+
 @article{schmelzer2023random,
   title     = {Random sets, copulas and related sets of probability measures},
   author    = {Schmelzer, Bernhard},

docs/src/api/gaussianprocesses.md

@@ -0,0 +1,28 @@
+# Gaussian Process Regression
+
+Methods for Gaussian process regression.
+
+## Index
+
+```@index
+Pages = ["gaussianprocesses.md"]
+```
+
+## Types
+
+```@docs
+GaussianProcess
+NoHyperparameterOptimization
+MaximumLikelihoodEstimation
+IdentityTransform
+ZScoreTransform
+UnitRangeTransform
+StandardNormalTransform
+```
+
+## Functions
+
+```@docs
+evaluate!(gp::GaussianProcess, data::DataFrame; mode::Symbol = :mean, n_samples::Int = 1)
+with_gaussian_noise(gp::AbstractGPs.GP, σ²::Real)
+```

docs/src/index.md

@@ -24,7 +24,7 @@ authors:
     link: https://github.com/mlsuh
   - name: Felix Mett
     platform: github
-    link: https://github.com/Cr0gan
+    link: https://github.com/felixmett
   - name: Andrea Perin
     platform: github
     link: https://github.com/andreaperin