Simulations
In addition to inference, the Gaussian process modelling framework allows simulations and predictions of the underlying process. This is done by conditioning the Gaussian process on some observations and then sampling from the conditioned distribution.
Sampling from the Gaussian process
Assuming a SingleBendingPowerLaw
we can draw realisations from the Gaussian process with a mean. First, we define the power spectral density function and plot it.
using Random
using Plots
using Pioran
rng = MersenneTwister(1234)
# Define power spectral density function
𝓟 = SingleBendingPowerLaw(0.3,1e-2,2.9)
f_min, f_max = 1e-3, 1e3
f0,fM = f_min/20.,f_max*20.
f = 10 .^ range(log10(f0), log10(fM), length=1000)
plot(f, 𝓟(f), xlabel="Frequency (day^-1)",ylabel="Power Spectral Density",legend=false,framestyle = :box,xscale=:log10,yscale=:log10,lw=2)
Then we approximate it to form a covariance function. We then build the Gaussian process and draw realisations from it using rand
.
variance = 15.2
# Approximation of the PSD to form a covariance function
𝓡 = approx(𝓟, f_min, f_max, 20, variance, basis_function="SHO")
μ = 1.3
# Build the GP
f = ScalableGP(μ, 𝓡) # Define the GP
T = 450 # days
t = range(0, stop=T, length=1000)
σ² = ones(length(t))*0.25
fx = f(t, σ²) # Gaussian process
realisations = [rand(rng,fx) for _ in 1:5]
plot(t, realisations, label="Realisations", xlabel="Time [days]", framestyle = :box,ylabel="Value")
Sampling from a Gaussian process built with semi-separable covariance functions is very efficient. The time complexity is O(N) where N is the number of data points.
Conditioning the Gaussian process
We can compute the conditioned or posterior distribution of the Gaussian process given some observations. Let's use a subset of the realisations to condition the Gaussian process and then sample from the conditioned distribution.
using StatsBase
idx = sort(sample(1:length(t), 50, replace = false));
t_obs = t[idx]
y_obs = realisations[1][idx]
yerr = 0.25*ones(length(t_obs))
fx = f(t_obs, yerr.^2) # Gaussian process
AbstractGPs.FiniteGP{ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Vector{Float64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}(
f: ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.ConstMean{Float64}(1.3), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.21179581122061306, 0.21179581122061306, 0.0002221441469079183, 0.0002221441469079183), Celerite(0.38542244831627853, 0.38542244831627853, 0.0006300497947253256, 0.0006300497947253256), Celerite(0.9038882037020084, 0.9038882037020084, 0.0017869601758986303, 0.0017869601758986303), Celerite(1.6067166095810141, 1.6067166095810141, 0.0050682131745472155, 0.0050682131745472155), Celerite(4.964666291074675, 4.964666291074675, 0.014374570362059993, 0.014374570362059993), Celerite(22.755812872815056, 22.755812872815056, 0.04076945187142277, 0.04076945187142277), Celerite(5.378025805869529, 5.378025805869529, 0.11563115724719714, 0.11563115724719714), Celerite(0.5417886639897328, 0.5417886639897328, 0.32795546451037977, 0.32795546451037977), Celerite(0.08070461417658467, 0.08070461417658467, 0.9301540282286335, 0.9301540282286335), Celerite(0.010848075938726597, 0.010848075938726597, 2.638121970376775, 2.638121970376775), Celerite(0.0015090806302273192, 0.0015090806302273192, 7.482295748198312, 7.482295748198312), Celerite(0.0002076692891940062, 0.0002076692891940062, 21.221440969050757, 21.221440969050757), Celerite(2.8676529647079397e-5, 2.8676529647079397e-5, 60.18868699641365, 60.18868699641365), Celerite(3.955543296634202e-6, 3.955543296634202e-6, 170.70839099171252, 170.70839099171252), Celerite(5.458037196298439e-7, 5.458037196298439e-7, 484.16664674402534, 484.16664674402534), Celerite(7.53042020169437e-8, 7.53042020169437e-8, 1373.203393562152, 1373.203393562152), Celerite(1.0389899331678754e-8, 1.0389899331678754e-8, 3894.70768540476, 3894.70768540476), Celerite(1.4338973791019396e-9, 1.4338973791019396e-9, 11046.249977144691, 11046.249977144691), Celerite(1.9677308815711412e-10, 1.9677308815711412e-10, 31329.601195702533, 31329.601195702533), Celerite(3.012647957419647e-11, 3.012647957419647e-11, 88857.65876316732, 88857.65876316732)], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732])), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.21179581122061306, 0.21179581122061306, 0.0002221441469079183, 0.0002221441469079183), Celerite(0.38542244831627853, 0.38542244831627853, 0.0006300497947253256, 0.0006300497947253256), Celerite(0.9038882037020084, 0.9038882037020084, 0.0017869601758986303, 0.0017869601758986303), Celerite(1.6067166095810141, 1.6067166095810141, 0.0050682131745472155, 0.0050682131745472155), Celerite(4.964666291074675, 4.964666291074675, 0.014374570362059993, 0.014374570362059993), Celerite(22.755812872815056, 22.755812872815056, 0.04076945187142277, 0.04076945187142277), Celerite(5.378025805869529, 5.378025805869529, 0.11563115724719714, 0.11563115724719714), Celerite(0.5417886639897328, 0.5417886639897328, 0.32795546451037977, 0.32795546451037977), Celerite(0.08070461417658467, 0.08070461417658467, 0.9301540282286335, 0.9301540282286335), Celerite(0.010848075938726597, 0.010848075938726597, 2.638121970376775, 2.638121970376775), Celerite(0.0015090806302273192, 0.0015090806302273192, 7.482295748198312, 7.482295748198312), Celerite(0.0002076692891940062, 0.0002076692891940062, 21.221440969050757, 21.221440969050757), Celerite(2.8676529647079397e-5, 2.8676529647079397e-5, 60.18868699641365, 60.18868699641365), Celerite(3.955543296634202e-6, 3.955543296634202e-6, 170.70839099171252, 170.70839099171252), Celerite(5.458037196298439e-7, 5.458037196298439e-7, 484.16664674402534, 484.16664674402534), Celerite(7.53042020169437e-8, 7.53042020169437e-8, 1373.203393562152, 1373.203393562152), Celerite(1.0389899331678754e-8, 1.0389899331678754e-8, 3894.70768540476, 3894.70768540476), Celerite(1.4338973791019396e-9, 1.4338973791019396e-9, 11046.249977144691, 11046.249977144691), Celerite(1.9677308815711412e-10, 1.9677308815711412e-10, 31329.601195702533, 31329.601195702533), Celerite(3.012647957419647e-11, 3.012647957419647e-11, 88857.65876316732, 88857.65876316732)], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732]))
x: [1.3513513513513513, 1.8018018018018018, 15.765765765765765, 44.14414414414414, 54.95495495495496, 58.55855855855856, 62.16216216216216, 85.58558558558559, 98.64864864864865, 100.45045045045045 … 396.39639639639637, 400.9009009009009, 401.8018018018018, 403.15315315315314, 415.31531531531533, 418.4684684684685, 422.52252252252254, 422.97297297297297, 425.22522522522524, 444.14414414414415]
Σy: [0.0625 0.0 … 0.0 0.0; 0.0 0.0625 … 0.0 0.0; … ; 0.0 0.0 … 0.0625 0.0; 0.0 0.0 … 0.0 0.0625]
)
We can compute the posterior distribution of the Gaussian process given the observations.
fp = posterior(fx, y_obs) # Posterior distribution
Pioran.PosteriorGP{AbstractGPs.FiniteGP{ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Vector{Float64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, Vector{Float64}}(AbstractGPs.FiniteGP{ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Vector{Float64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}(
f: ScalableGP{AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.GP{AbstractGPs.ConstMean{Float64}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.ConstMean{Float64}(1.3), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.21179581122061306, 0.21179581122061306, 0.0002221441469079183, 0.0002221441469079183), Celerite(0.38542244831627853, 0.38542244831627853, 0.0006300497947253256, 0.0006300497947253256), Celerite(0.9038882037020084, 0.9038882037020084, 0.0017869601758986303, 0.0017869601758986303), Celerite(1.6067166095810141, 1.6067166095810141, 0.0050682131745472155, 0.0050682131745472155), Celerite(4.964666291074675, 4.964666291074675, 0.014374570362059993, 0.014374570362059993), Celerite(22.755812872815056, 22.755812872815056, 0.04076945187142277, 0.04076945187142277), Celerite(5.378025805869529, 5.378025805869529, 0.11563115724719714, 0.11563115724719714), Celerite(0.5417886639897328, 0.5417886639897328, 0.32795546451037977, 0.32795546451037977), Celerite(0.08070461417658467, 0.08070461417658467, 0.9301540282286335, 0.9301540282286335), Celerite(0.010848075938726597, 0.010848075938726597, 2.638121970376775, 2.638121970376775), Celerite(0.0015090806302273192, 0.0015090806302273192, 7.482295748198312, 7.482295748198312), Celerite(0.0002076692891940062, 0.0002076692891940062, 21.221440969050757, 21.221440969050757), Celerite(2.8676529647079397e-5, 2.8676529647079397e-5, 60.18868699641365, 60.18868699641365), Celerite(3.955543296634202e-6, 3.955543296634202e-6, 170.70839099171252, 170.70839099171252), Celerite(5.458037196298439e-7, 5.458037196298439e-7, 484.16664674402534, 484.16664674402534), Celerite(7.53042020169437e-8, 7.53042020169437e-8, 1373.203393562152, 1373.203393562152), Celerite(1.0389899331678754e-8, 1.0389899331678754e-8, 3894.70768540476, 3894.70768540476), Celerite(1.4338973791019396e-9, 1.4338973791019396e-9, 11046.249977144691, 11046.249977144691), Celerite(1.9677308815711412e-10, 1.9677308815711412e-10, 31329.601195702533, 31329.601195702533), Celerite(3.012647957419647e-11, 3.012647957419647e-11, 88857.65876316732, 88857.65876316732)], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732])), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.21179581122061306, 0.21179581122061306, 0.0002221441469079183, 0.0002221441469079183), Celerite(0.38542244831627853, 0.38542244831627853, 0.0006300497947253256, 0.0006300497947253256), Celerite(0.9038882037020084, 0.9038882037020084, 0.0017869601758986303, 0.0017869601758986303), Celerite(1.6067166095810141, 1.6067166095810141, 0.0050682131745472155, 0.0050682131745472155), Celerite(4.964666291074675, 4.964666291074675, 0.014374570362059993, 0.014374570362059993), Celerite(22.755812872815056, 22.755812872815056, 0.04076945187142277, 0.04076945187142277), Celerite(5.378025805869529, 5.378025805869529, 0.11563115724719714, 0.11563115724719714), Celerite(0.5417886639897328, 0.5417886639897328, 0.32795546451037977, 0.32795546451037977), Celerite(0.08070461417658467, 0.08070461417658467, 0.9301540282286335, 0.9301540282286335), Celerite(0.010848075938726597, 0.010848075938726597, 2.638121970376775, 2.638121970376775), Celerite(0.0015090806302273192, 0.0015090806302273192, 7.482295748198312, 7.482295748198312), Celerite(0.0002076692891940062, 0.0002076692891940062, 21.221440969050757, 21.221440969050757), Celerite(2.8676529647079397e-5, 2.8676529647079397e-5, 60.18868699641365, 60.18868699641365), Celerite(3.955543296634202e-6, 3.955543296634202e-6, 170.70839099171252, 170.70839099171252), Celerite(5.458037196298439e-7, 5.458037196298439e-7, 484.16664674402534, 484.16664674402534), Celerite(7.53042020169437e-8, 7.53042020169437e-8, 1373.203393562152, 1373.203393562152), Celerite(1.0389899331678754e-8, 1.0389899331678754e-8, 3894.70768540476, 3894.70768540476), Celerite(1.4338973791019396e-9, 1.4338973791019396e-9, 11046.249977144691, 11046.249977144691), Celerite(1.9677308815711412e-10, 1.9677308815711412e-10, 31329.601195702533, 31329.601195702533), Celerite(3.012647957419647e-11, 3.012647957419647e-11, 88857.65876316732, 88857.65876316732)], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.21179581122061306, 0.38542244831627853, 0.9038882037020084, 1.6067166095810141, 4.964666291074675, 22.755812872815056, 5.378025805869529, 0.5417886639897328, 0.08070461417658467, 0.010848075938726597, 0.0015090806302273192, 0.0002076692891940062, 2.8676529647079397e-5, 3.955543296634202e-6, 5.458037196298439e-7, 7.53042020169437e-8, 1.0389899331678754e-8, 1.4338973791019396e-9, 1.9677308815711412e-10, 3.012647957419647e-11], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732], [0.0002221441469079183, 0.0006300497947253256, 0.0017869601758986303, 0.0050682131745472155, 0.014374570362059993, 0.04076945187142277, 0.11563115724719714, 0.32795546451037977, 0.9301540282286335, 2.638121970376775, 7.482295748198312, 21.221440969050757, 60.18868699641365, 170.70839099171252, 484.16664674402534, 1373.203393562152, 3894.70768540476, 11046.249977144691, 31329.601195702533, 88857.65876316732]))
x: [1.3513513513513513, 1.8018018018018018, 15.765765765765765, 44.14414414414414, 54.95495495495496, 58.55855855855856, 62.16216216216216, 85.58558558558559, 98.64864864864865, 100.45045045045045 … 396.39639639639637, 400.9009009009009, 401.8018018018018, 403.15315315315314, 415.31531531531533, 418.4684684684685, 422.52252252252254, 422.97297297297297, 425.22522522522524, 444.14414414414415]
Σy: [0.0625 0.0 … 0.0 0.0; 0.0 0.0625 … 0.0 0.0; … ; 0.0 0.0 … 0.0625 0.0; 0.0 0.0 … 0.0 0.0625]
)
, [5.049625421463632, 5.788046102707981, 3.253813167202233, -6.434451728276888, -2.153775909938405, 1.0264990265324576, 4.751459874105995, 4.31344119799946, 9.309942440653156, 9.295111340218552 … 2.5120861273426445, 0.07640695112513551, 0.7411102972888532, -0.19778064914068172, -8.3981535447054, -8.348219512676872, -8.300576243904722, -8.754491553645984, -8.714588428544848, 0.7835131871531228])
The mean and standard deviation of this distribution, can be computed using the mean
and std
functions. The posterior covariance matrix can be computed using the cov
function.
m = mean(fp,t);
s = std(fp,t);
1000-element Vector{Float64}:
0.8128885705176537
0.6108530914946403
0.4006365343352816
0.2120747150299907
0.21080303025566322
0.3912678115025849
0.5858127957039118
0.7648927567581641
0.9277798746859259
1.0747445973130927
⋮
1.2244987023111926
1.405812711700324
1.5784983878280654
1.7434754428873536
1.9013910043560451
2.052690740509442
2.1977093081230383
2.336733400913537
2.4700354318778466
We can plot the realisations, the observations and the posterior distribution.
plot(t, realisations[1], label="Realisation", xlabel="Time [days]", framestyle = :box,ylabel="Value")
plot!(t_obs, y_obs,yerr=yerr, label="Observations", seriestype=:scatter)
plot!(t,m, ribbon=s,label="Posterior distribution", lw=2)
plot!(t,m,ribbon=2*s, label=nothing)
The computation of the mean of the distribution is very efficient. The time complexity is O(N) where N is the number of data points. However, the computation of the covariance is very inefficient as the posterior covariance is not semi-separable.
Sampling from the conditioned distribution
We can draw realisations from the conditioned distribution using rand
.
samples_cond = rand(rng,fp,t,5);
1000×5 Matrix{Float64}:
5.53396 5.10929 6.40639 6.23271 4.27956
5.71575 4.95406 6.0142 6.00398 4.52487
5.63491 4.84978 5.81612 5.73105 4.71497
5.63823 4.9363 5.81306 5.62295 5.20891
5.39262 5.44513 5.93852 5.44746 5.3652
5.05934 5.54439 5.62429 5.56458 4.9885
4.78084 5.68507 5.5569 5.20994 4.46708
4.51815 5.60471 5.6951 4.94701 4.20013
4.01389 5.99095 6.0906 4.56868 3.97646
3.40165 6.32814 5.91974 4.38783 3.47988
⋮
0.189211 1.31306 -0.315965 2.20256 2.43656
0.498496 1.01516 -0.329997 2.28257 2.53727
0.82524 0.245869 -0.430343 2.72871 2.77738
0.994561 0.0734135 -0.297912 3.10073 2.47235
1.38801 -0.199538 -0.313239 2.93133 1.75012
1.75062 -0.584455 -0.355962 3.10093 1.47234
2.12952 -0.374646 -0.299121 3.11443 1.27404
2.65848 -0.501399 -0.178592 3.26339 1.48661
3.0586 -0.472663 -0.842167 3.62231 1.34132
We can plot the realisations, the observations and the posterior distribution.
plot(t_obs, y_obs,yerr=yerr, label="Observations", seriestype=:scatter, xlabel="Time [days]", framestyle = :box,ylabel="Value",color=:black,lw=2)
plot!(t, samples_cond, label=nothing)