gp.jl

t=0:0.25:10
n=length(t)
alpha=1.0
mu=zeros(n)
Sigma=zeros(n,n)

for i=1:n
    for j=1:n 
        Sigma[i,j]=exp(-abs2(t[i]-t[j]) / 2 * 0.5);
    end
end

using LinearAlgebra

function symfact(M, tol=1.0e-14)
    Asvd = svd(Hermitian(M))
    return Asvd.U * Diagonal(sqrt.(clamp.(Asvd.S,tol,Inf)))
end

A=symfact(Sigma)

using PyPlot
pygui(true)

for i in 1:10
    X=mu+A*randn(n)
    plot(t,X)
end

#simple inference
using ExtremalOptimization
initpoint(i) = randn(n)
function cost(x)
    f=A*x
    (f[1+end÷2]-0.0)^2+(f[1+end÷2+end÷4]-1.5)^2+(f[1+end÷4]-1.0)^2+0.5*(1-norm(x)/sqrt(length(x)))^2
end
using Statistics
SS=[begin
    local sol=optimize(cost, initpoint, 2*length(t))
    A*sol.x
end for i in 1:50]

μ=mean(SS)

subplot(1,2,1)
for i in 1:10
    X=mu+A*randn(n)
    plot(t,X)
end
subplot(1,2,2)
for so in SS
    plot(t,so)
end
plot(t,μ,lw=2,color="black")