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.