Gaussian Process Sampling

gaussianprocess.jl

using LinearAlgebra
using Distributions
using Plots
plotly()

m(x) = 0.25*x^2
k(x,y) = exp(-0.5 * (x-y)^2)

xmin = -5
xmax = 5
ns = 51
keps = 1e-9
xs = range(xmin, xmax, length=ns)

mu = m.(xs)
sigma = k.(xs, xs') + keps*Matrix{Float64}(I, ns, ns)

d = MvNormal(mu, sigma)

display(scatter(xs, rand(d)))