Sampling from a GP

Preliminary steps

Loading necessary packages

using AugmentedGaussianProcesses
using Distributions
using LinearAlgebra
using Plots
default(lw=3.0)

Generating some random binary data

kernel = with_lengthscale(SqExponentialKernel(), 1.0)
N_train = 50
x = range(-10, 10; length=50)
x_test = range(-10, 10; length=500)
K = kernelmatrix(kernel, vcat(x, x_test))
f_all = rand(MvNormal(K + 1e-8I)) # Sample a random GP
f = f_all[1:N_train]
y = rand.(Bernoulli.(AGP.logistic.(f)))
y_sign = Int.(sign.(y .- 0.5));

We create a function to visualize the data

function plot_data(x, y; size=(300, 500), kwargs...)
    return Plots.scatter(x, y; alpha=0.5, markerstrokewidth=0.0, lab="", size=size, kwargs...)
end
plot_data(x, y; size=(600, 500), xlabel="x", ylabel="y")

Model initialization and training

Run the variational gaussian process approximation

@info "Running full model"
mfull = VGP(x, y_sign, kernel, LogisticLikelihood(), AnalyticVI(); optimiser=false)
@time train!(mfull, 5)
(Variational Gaussian Process with a BernoulliLikelihood{GPLikelihoods.LogisticLink}(GPLikelihoods.LogisticLink(LogExpFunctions.logistic)) infered by Analytic Variational Inference , (local_vars = (c = [0.8142691922636827, 0.752563365282887, 0.7894929148703664, 0.9528403689976539, 1.1702435934033402, 1.358121221102231, 1.4782841218893823, 1.536593504500029, 1.5588578880365507, 1.563685819333194  …  1.3724685843205982, 1.4088625811952509, 1.2945272245162387, 1.0354804103255384, 0.7713378409056546, 0.7509189635077701, 0.9002201941065535, 0.9616125735281037, 0.9029063985647964, 0.8466512533378797], θ = [0.23704506843197112, 0.23883302749216445, 0.2377759643785148, 0.2326581109205659, 0.22490107618682584, 0.21754624892847121, 0.21262033657618967, 0.21018504523097856, 0.20924896923798728, 0.20904557935501003  …  0.21696576601384443, 0.2154833870209969, 0.22008998939590535, 0.22982280004040592, 0.2383005406478847, 0.2388791766634075, 0.2343812678952892, 0.23236442530673088, 0.23429494173304496, 0.23606437027515187]), opt_state = (NamedTuple(),), hyperopt_state = (NamedTuple(),), kernel_matrices = ((K = LinearAlgebra.Cholesky{Float64, Matrix{Float64}}([1.0000499987500624 0.9200303474974239 … 4.46855233855175e-84 1.3838273370995804e-87; 0.9200763478648182 0.3919746926572805 … 3.1154230234065755e-80 1.1397426112260664e-83; … ; 4.468775760583266e-84 1.2215781024732928e-80 … 0.09632063884559038 0.29129437964209737; 1.3838965267367376e-87 4.468775760583266e-84 … 0.9200763478648182 0.09632063884675626], 'U', 0),),)))

We can also create a sampling based model

@info "Sampling from model"
mmcmc = MCGP(x, y, kernel, LogisticLikelihood(), GibbsSampling(); optimiser=false)
m = mmcmc
@time samples = sample(mmcmc, 1000);
[ Info: Sampling from model
Sampling with Gibbs Sampler   0%|                       |  ETA: N/A
Sampling with Gibbs Sampler   0%|▏                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   1%|▎                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   1%|▍                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   2%|▍                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   2%|▌                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   3%|▋                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   3%|▊                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   4%|▉                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   4%|█                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   5%|█                      |  ETA: 0:00:00
Sampling with Gibbs Sampler   5%|█▏                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   5%|█▎                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   6%|█▍                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   6%|█▌                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   7%|█▋                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   7%|█▋                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   8%|█▊                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   8%|█▉                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   9%|██                     |  ETA: 0:00:00
Sampling with Gibbs Sampler   9%|██▏                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  10%|██▎                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  10%|██▎                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  10%|██▍                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  11%|██▌                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  11%|██▋                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  12%|██▊                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  12%|██▉                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  13%|██▉                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  13%|███                    |  ETA: 0:00:00
Sampling with Gibbs Sampler  14%|███▏                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  14%|███▎                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  15%|███▍                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  15%|███▌                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  15%|███▌                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  16%|███▋                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  16%|███▊                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  17%|███▉                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  17%|████                   |  ETA: 0:00:00
Sampling with Gibbs Sampler  18%|████▏                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  18%|████▏                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  19%|████▎                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  19%|████▍                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  20%|████▌                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  20%|████▋                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  20%|████▊                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  21%|████▊                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  21%|████▉                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  22%|█████                  |  ETA: 0:00:00
Sampling with Gibbs Sampler  22%|█████▏                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  23%|█████▎                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  23%|█████▍                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  24%|█████▍                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  24%|█████▌                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  25%|█████▋                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  25%|█████▊                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  25%|█████▉                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  26%|██████                 |  ETA: 0:00:00
Sampling with Gibbs Sampler  26%|██████▏                |  ETA: 0:00:00
Sampling with Gibbs Sampler  27%|██████▏                |  ETA: 0:00:00
Sampling with Gibbs Sampler  27%|██████▎                |  ETA: 0:00:00
Sampling with Gibbs Sampler  28%|██████▍                |  ETA: 0:00:00
Sampling with Gibbs Sampler  28%|██████▌                |  ETA: 0:00:00
Sampling with Gibbs Sampler  29%|██████▋                |  ETA: 0:00:00
Sampling with Gibbs Sampler  29%|██████▊                |  ETA: 0:00:00
Sampling with Gibbs Sampler  30%|██████▊                |  ETA: 0:00:00
Sampling with Gibbs Sampler  30%|██████▉                |  ETA: 0:00:00
Sampling with Gibbs Sampler  30%|███████                |  ETA: 0:00:00
Sampling with Gibbs Sampler  31%|███████▏               |  ETA: 0:00:00
Sampling with Gibbs Sampler  31%|███████▎               |  ETA: 0:00:00
Sampling with Gibbs Sampler  32%|███████▍               |  ETA: 0:00:00
Sampling with Gibbs Sampler  32%|███████▍               |  ETA: 0:00:00
Sampling with Gibbs Sampler  33%|███████▌               |  ETA: 0:00:00
Sampling with Gibbs Sampler  33%|███████▋               |  ETA: 0:00:00
Sampling with Gibbs Sampler  34%|███████▊               |  ETA: 0:00:00
Sampling with Gibbs Sampler  34%|███████▉               |  ETA: 0:00:00
Sampling with Gibbs Sampler  35%|████████               |  ETA: 0:00:00
Sampling with Gibbs Sampler  35%|████████               |  ETA: 0:00:00
Sampling with Gibbs Sampler  35%|████████▏              |  ETA: 0:00:00
Sampling with Gibbs Sampler  36%|████████▎              |  ETA: 0:00:00
Sampling with Gibbs Sampler  36%|████████▍              |  ETA: 0:00:00
Sampling with Gibbs Sampler  37%|████████▌              |  ETA: 0:00:00
Sampling with Gibbs Sampler  37%|████████▋              |  ETA: 0:00:00
Sampling with Gibbs Sampler  38%|████████▋              |  ETA: 0:00:00
Sampling with Gibbs Sampler  38%|████████▊              |  ETA: 0:00:00
Sampling with Gibbs Sampler  39%|████████▉              |  ETA: 0:00:00
Sampling with Gibbs Sampler  39%|█████████              |  ETA: 0:00:00
Sampling with Gibbs Sampler  40%|█████████▏             |  ETA: 0:00:00
Sampling with Gibbs Sampler  40%|█████████▎             |  ETA: 0:00:00
Sampling with Gibbs Sampler  40%|█████████▎             |  ETA: 0:00:00
Sampling with Gibbs Sampler  41%|█████████▍             |  ETA: 0:00:00
Sampling with Gibbs Sampler  41%|█████████▌             |  ETA: 0:00:00
Sampling with Gibbs Sampler  42%|█████████▋             |  ETA: 0:00:00
Sampling with Gibbs Sampler  42%|█████████▊             |  ETA: 0:00:00
Sampling with Gibbs Sampler  43%|█████████▉             |  ETA: 0:00:00
Sampling with Gibbs Sampler  43%|█████████▉             |  ETA: 0:00:00
Sampling with Gibbs Sampler  44%|██████████             |  ETA: 0:00:00
Sampling with Gibbs Sampler  44%|██████████▏            |  ETA: 0:00:00
Sampling with Gibbs Sampler  45%|██████████▎            |  ETA: 0:00:00
Sampling with Gibbs Sampler  45%|██████████▍            |  ETA: 0:00:00
Sampling with Gibbs Sampler  45%|██████████▌            |  ETA: 0:00:00
Sampling with Gibbs Sampler  46%|██████████▌            |  ETA: 0:00:00
Sampling with Gibbs Sampler  46%|██████████▋            |  ETA: 0:00:00
Sampling with Gibbs Sampler  47%|██████████▊            |  ETA: 0:00:00
Sampling with Gibbs Sampler  47%|██████████▉            |  ETA: 0:00:00
Sampling with Gibbs Sampler  48%|███████████            |  ETA: 0:00:00
Sampling with Gibbs Sampler  48%|███████████▏           |  ETA: 0:00:00
Sampling with Gibbs Sampler  49%|███████████▏           |  ETA: 0:00:00
Sampling with Gibbs Sampler  49%|███████████▎           |  ETA: 0:00:00
Sampling with Gibbs Sampler  50%|███████████▍           |  ETA: 0:00:00
Sampling with Gibbs Sampler  50%|███████████▌           |  ETA: 0:00:00
Sampling with Gibbs Sampler  50%|███████████▋           |  ETA: 0:00:00
Sampling with Gibbs Sampler  51%|███████████▊           |  ETA: 0:00:00
Sampling with Gibbs Sampler  51%|███████████▉           |  ETA: 0:00:00
Sampling with Gibbs Sampler  52%|███████████▉           |  ETA: 0:00:00
Sampling with Gibbs Sampler  52%|████████████           |  ETA: 0:00:00
Sampling with Gibbs Sampler  53%|████████████▏          |  ETA: 0:00:00
Sampling with Gibbs Sampler  53%|████████████▎          |  ETA: 0:00:00
Sampling with Gibbs Sampler  54%|████████████▍          |  ETA: 0:00:00
Sampling with Gibbs Sampler  54%|████████████▌          |  ETA: 0:00:00
Sampling with Gibbs Sampler  55%|████████████▌          |  ETA: 0:00:00
Sampling with Gibbs Sampler  55%|████████████▋          |  ETA: 0:00:00
Sampling with Gibbs Sampler  55%|████████████▊          |  ETA: 0:00:00
Sampling with Gibbs Sampler  56%|████████████▉          |  ETA: 0:00:00
Sampling with Gibbs Sampler  56%|█████████████          |  ETA: 0:00:00
Sampling with Gibbs Sampler  57%|█████████████▏         |  ETA: 0:00:00
Sampling with Gibbs Sampler  57%|█████████████▏         |  ETA: 0:00:00
Sampling with Gibbs Sampler  58%|█████████████▎         |  ETA: 0:00:00
Sampling with Gibbs Sampler  58%|█████████████▍         |  ETA: 0:00:00
Sampling with Gibbs Sampler  59%|█████████████▌         |  ETA: 0:00:00
Sampling with Gibbs Sampler  59%|█████████████▋         |  ETA: 0:00:00
Sampling with Gibbs Sampler  60%|█████████████▊         |  ETA: 0:00:00
Sampling with Gibbs Sampler  60%|█████████████▊         |  ETA: 0:00:00
Sampling with Gibbs Sampler  60%|█████████████▉         |  ETA: 0:00:00
Sampling with Gibbs Sampler  61%|██████████████         |  ETA: 0:00:00
Sampling with Gibbs Sampler  61%|██████████████▏        |  ETA: 0:00:00
Sampling with Gibbs Sampler  62%|██████████████▎        |  ETA: 0:00:00
Sampling with Gibbs Sampler  62%|██████████████▍        |  ETA: 0:00:00
Sampling with Gibbs Sampler  63%|██████████████▍        |  ETA: 0:00:00
Sampling with Gibbs Sampler  63%|██████████████▌        |  ETA: 0:00:00
Sampling with Gibbs Sampler  64%|██████████████▋        |  ETA: 0:00:00
Sampling with Gibbs Sampler  64%|██████████████▊        |  ETA: 0:00:00
Sampling with Gibbs Sampler  65%|██████████████▉        |  ETA: 0:00:00
Sampling with Gibbs Sampler  65%|███████████████        |  ETA: 0:00:00
Sampling with Gibbs Sampler  65%|███████████████        |  ETA: 0:00:00
Sampling with Gibbs Sampler  66%|███████████████▏       |  ETA: 0:00:00
Sampling with Gibbs Sampler  66%|███████████████▎       |  ETA: 0:00:00
Sampling with Gibbs Sampler  67%|███████████████▍       |  ETA: 0:00:00
Sampling with Gibbs Sampler  67%|███████████████▌       |  ETA: 0:00:00
Sampling with Gibbs Sampler  68%|███████████████▋       |  ETA: 0:00:00
Sampling with Gibbs Sampler  68%|███████████████▋       |  ETA: 0:00:00
Sampling with Gibbs Sampler  69%|███████████████▊       |  ETA: 0:00:00
Sampling with Gibbs Sampler  69%|███████████████▉       |  ETA: 0:00:00
Sampling with Gibbs Sampler  70%|████████████████       |  ETA: 0:00:00
Sampling with Gibbs Sampler  70%|████████████████▏      |  ETA: 0:00:00
Sampling with Gibbs Sampler  70%|████████████████▎      |  ETA: 0:00:00
Sampling with Gibbs Sampler  71%|████████████████▎      |  ETA: 0:00:00
Sampling with Gibbs Sampler  71%|████████████████▍      |  ETA: 0:00:00
Sampling with Gibbs Sampler  72%|████████████████▌      |  ETA: 0:00:00
Sampling with Gibbs Sampler  72%|████████████████▋      |  ETA: 0:00:00
Sampling with Gibbs Sampler  73%|████████████████▊      |  ETA: 0:00:00
Sampling with Gibbs Sampler  73%|████████████████▉      |  ETA: 0:00:00
Sampling with Gibbs Sampler  74%|████████████████▉      |  ETA: 0:00:00
Sampling with Gibbs Sampler  74%|█████████████████      |  ETA: 0:00:00
Sampling with Gibbs Sampler  75%|█████████████████▏     |  ETA: 0:00:00
Sampling with Gibbs Sampler  75%|█████████████████▎     |  ETA: 0:00:00
Sampling with Gibbs Sampler  75%|█████████████████▍     |  ETA: 0:00:00
Sampling with Gibbs Sampler  76%|█████████████████▌     |  ETA: 0:00:00
Sampling with Gibbs Sampler  76%|█████████████████▋     |  ETA: 0:00:00
Sampling with Gibbs Sampler  77%|█████████████████▋     |  ETA: 0:00:00
Sampling with Gibbs Sampler  77%|█████████████████▊     |  ETA: 0:00:00
Sampling with Gibbs Sampler  78%|█████████████████▉     |  ETA: 0:00:00
Sampling with Gibbs Sampler  78%|██████████████████     |  ETA: 0:00:00
Sampling with Gibbs Sampler  79%|██████████████████▏    |  ETA: 0:00:00
Sampling with Gibbs Sampler  79%|██████████████████▎    |  ETA: 0:00:00
Sampling with Gibbs Sampler  80%|██████████████████▎    |  ETA: 0:00:00
Sampling with Gibbs Sampler  80%|██████████████████▍    |  ETA: 0:00:00
Sampling with Gibbs Sampler  80%|██████████████████▌    |  ETA: 0:00:00
Sampling with Gibbs Sampler  81%|██████████████████▋    |  ETA: 0:00:00
Sampling with Gibbs Sampler  81%|██████████████████▊    |  ETA: 0:00:00
Sampling with Gibbs Sampler  82%|██████████████████▉    |  ETA: 0:00:00
Sampling with Gibbs Sampler  82%|██████████████████▉    |  ETA: 0:00:00
Sampling with Gibbs Sampler  83%|███████████████████    |  ETA: 0:00:00
Sampling with Gibbs Sampler  83%|███████████████████▏   |  ETA: 0:00:00
Sampling with Gibbs Sampler  84%|███████████████████▎   |  ETA: 0:00:00
Sampling with Gibbs Sampler  84%|███████████████████▍   |  ETA: 0:00:00
Sampling with Gibbs Sampler  85%|███████████████████▌   |  ETA: 0:00:00
Sampling with Gibbs Sampler  85%|███████████████████▌   |  ETA: 0:00:00
Sampling with Gibbs Sampler  85%|███████████████████▋   |  ETA: 0:00:00
Sampling with Gibbs Sampler  86%|███████████████████▊   |  ETA: 0:00:00
Sampling with Gibbs Sampler  86%|███████████████████▉   |  ETA: 0:00:00
Sampling with Gibbs Sampler  87%|████████████████████   |  ETA: 0:00:00
Sampling with Gibbs Sampler  87%|████████████████████▏  |  ETA: 0:00:00
Sampling with Gibbs Sampler  88%|████████████████████▏  |  ETA: 0:00:00
Sampling with Gibbs Sampler  88%|████████████████████▎  |  ETA: 0:00:00
Sampling with Gibbs Sampler  89%|████████████████████▍  |  ETA: 0:00:00
Sampling with Gibbs Sampler  89%|████████████████████▌  |  ETA: 0:00:00
Sampling with Gibbs Sampler  90%|████████████████████▋  |  ETA: 0:00:00
Sampling with Gibbs Sampler  90%|████████████████████▊  |  ETA: 0:00:00
Sampling with Gibbs Sampler  90%|████████████████████▊  |  ETA: 0:00:00
Sampling with Gibbs Sampler  91%|████████████████████▉  |  ETA: 0:00:00
Sampling with Gibbs Sampler  91%|█████████████████████  |  ETA: 0:00:00
Sampling with Gibbs Sampler  92%|█████████████████████▏ |  ETA: 0:00:00
Sampling with Gibbs Sampler  92%|█████████████████████▎ |  ETA: 0:00:00
Sampling with Gibbs Sampler  93%|█████████████████████▍ |  ETA: 0:00:00
Sampling with Gibbs Sampler  93%|█████████████████████▍ |  ETA: 0:00:00
Sampling with Gibbs Sampler  94%|█████████████████████▌ |  ETA: 0:00:00
Sampling with Gibbs Sampler  94%|█████████████████████▋ |  ETA: 0:00:00
Sampling with Gibbs Sampler  95%|█████████████████████▊ |  ETA: 0:00:00
Sampling with Gibbs Sampler  95%|█████████████████████▉ |  ETA: 0:00:00
Sampling with Gibbs Sampler  95%|██████████████████████ |  ETA: 0:00:00
Sampling with Gibbs Sampler  96%|██████████████████████ |  ETA: 0:00:00
Sampling with Gibbs Sampler  96%|██████████████████████▏|  ETA: 0:00:00
Sampling with Gibbs Sampler  97%|██████████████████████▎|  ETA: 0:00:00
Sampling with Gibbs Sampler  97%|██████████████████████▍|  ETA: 0:00:00
Sampling with Gibbs Sampler  98%|██████████████████████▌|  ETA: 0:00:00
Sampling with Gibbs Sampler  98%|██████████████████████▋|  ETA: 0:00:00
Sampling with Gibbs Sampler  99%|██████████████████████▋|  ETA: 0:00:00
Sampling with Gibbs Sampler  99%|██████████████████████▊|  ETA: 0:00:00
Sampling with Gibbs Sampler 100%|██████████████████████▉|  ETA: 0:00:00
Sampling with Gibbs Sampler 100%|███████████████████████| Time: 0:00:00
Sampling with Gibbs Sampler 100%|███████████████████████| Time: 0:00:00
  0.207177 seconds (196.13 k allocations: 146.852 MiB)

We can now visualize the results of both models

We first plot the latent function f (truth, the VI estimate, the samples)

p1 = plot(x, samples; label="", color=:black, alpha=0.01, lab="")
plot!(x, mean(mfull[1]); color=:blue, ribbon=sqrt.(var(mfull[1])), label="VI")
plot!(x_test, f_all[N_train+1:end]; color=:red, label="true f")

And we can also plot the predictions vs the data

p2 = plot_data(x, y; size=(600, 400))
μ_vi, σ_vi = proba_y(mfull, x_test)
plot!(x_test, μ_vi; ribbon=σ_vi, label="VI")
μ_mcmc, σ_mcmc = proba_y(mmcmc, x_test)
plot!(x_test, μ_mcmc; ribbon=σ_mcmc, label="MCMC")

This page was generated using Literate.jl.