Modelling
Modelling the power spectral density
We use SingleBendingPowerLaw and DoubleBendingPowerLaw to model the power spectral density of the random process generating the time series data.
using Plots
using Pioran
𝓟1 = SingleBendingPowerLaw(1., .1, 3.4)
𝓟2 = SingleBendingPowerLaw(.4, 1e-2, 3.)
𝓟3 = SingleBendingPowerLaw(.1, 3., 2.4)
𝓟1d = DoubleBendingPowerLaw(1.2,1e-3,2.4,1.1,4.2)
𝓟2d = DoubleBendingPowerLaw(0.2,1e-2,1.3,24.3,3.1)
𝓟3d = DoubleBendingPowerLaw(0.5,1e-3,2.1,91.2,4.9)
f = 10 .^ range(-4, stop=3, length=1000)
l = @layout [a b]
p1 = plot(f,[f.*𝓟1(f),f.*𝓟2(f),f.*𝓟3(f)],xlabel="Frequency",ylabel="Frequency * Power",legend=false,framestyle = :box,xscale=:log10,yscale=:log10,ylims=(1e-15,1e1),lw=2)
p2 = plot(f,[f.*𝓟1d(f),f.*𝓟2d(f),f.*𝓟3d(f)],xlabel="Frequency",legend=false,framestyle = :box,xscale=:log10,yscale=:log10,ylims=(1e-15,1e1),lw=2)
plot(p1,p2,layout=l,size=(700,300),grid=false,left_margin=2Plots.mm,bottom_margin=20Plots.px,title=["Single bending power-law" "Double bending power-law"])
Approximating the power spectral density
To obtain the covariance function, we approximate the power spectral density by a sum of basis functions $\psi(f)$. At the moment, we use either of the following basis functions:
\[\begin{align}\begin{split} \psi_4(f) = \frac{1}{1+f^4}\\ \psi_6(f) = \frac{1}{1+f^6} \end{split} \end{align}\]
These basis functions have analytical Fourier transforms and are used to approximate the covariance function. The Fourier transform of the basis functions are given by
\[\begin{align}\begin{split} \phi_4(\tau) &= \dfrac{\pi}{\sqrt2} \exp\left(-\pi\sqrt{2}|\tau|\right) \left(\cos\left(\pi\sqrt{2}|\tau|\right)+\sin\left(\pi\sqrt{2}|\tau|\right)\right)\\ \phi_6(\tau) &=\dfrac{\pi}{3}\exp{\left(-2\pi |\tau|\right)}+\exp{\left(-\pi|\tau|\right)}\left(\pi/3\cos\left(\pi\sqrt{3}|\tau|\right)+\pi/\sqrt{3}\sin\left(\pi\sqrt{3}|\tau|\right)\right)\end{split} \end{align}\]
We need to specify the frequency range $f0=f_\mathrm{min}/S_\mathrm{low}$ and $fM=f_\mathrm{max}S_\mathrm{high}$ over which the approximation is performed. We also need to specify the number of basis functions $J$ to use. Once this is done the frequency grid is defined as:
\[f_j=f_0\left({f_M}/{f_0}\right)^{j/(J-1)}, ~~j=0,1,2,\dots,J-1\]
The approximation of the power spectral density is then given by:
\[\begin{align} \mathcal{P}(f)\simeq \tilde{\mathcal{P}}(f)&= \sum\limits_{j=0}^{J-1} a_j \psi(f/f_j)\\ \mathcal{R}(\tau)\simeq \tilde{\mathcal{R}}(\tau)&= \sum\limits_{j=0}^{J-1} a_j f_j \phi(\tau f_j) \end{align}\]
Adding the constraint that the approximation and the true power spectrum must be equal on the grid of frequencies gives a linear system of $J$ equations for the coefficients $a_j$. This system can be written with a Toeplitz matrix $B$ and a vector $\boldsymbol{a}$ as:
\[\begin{align} \boldsymbol{p} = \boldsymbol{a} B \quad \text{where } B_{ij}=\psi(f_i/f_j) \text{ and } p_j = \mathcal{P}(f_j)/\mathcal{P}(f_0) \end{align}\]
The values of $p_j$ are divided by $p_0$ so that the values of $a_j$ are not too high. Visually, the approximation can be seen as follows:
f0, fM = 1e-3, 1e3
𝓟 = SingleBendingPowerLaw(.4, 1e-1, 3.)
f = 10 .^ range(-3, stop=3, length=1000)
psd_approx = Pioran.approximated_psd(f, 𝓟, f0, fM, n_components=20,basis_function="SHO",individual=true)
plot(f,𝓟(f)/𝓟(f0),xscale=:log10,yscale=:log10,label="True PSD",xlabel="Frequency (day^-1)",ylabel="Power Spectral Density",lw=2,framestyle = :box,grid=false)
plot!(f,sum(psd_approx,dims=2),label="Approximated PSD",lw=2)
plot!(f,psd_approx,label=nothing,color=:black,alpha=.5,ylims=(1e-15,1e1),ls=:dot)
Limitations and diagnostics for the approximation
This approximation is limited by the steepness of the basis functions, this means that if the power spectrum you want to approximate is steeper than the basis functions, the approximation may fail. Equivalently, the basis functions are flat at low frequencies, modelling a rising power spectrum at low frequencies can also be difficult.
In order to check the quality of the approximation, we can compute the residuals and ratios between the true and approximated power spectrum. First, we need to define the range of allowed values for each parameter to check. As we adopt a Bayesian workflow, one can use the prior distribution to define the range of allowed values for each parameter. This can be done as follows:
using Distributions
using Random
rng = MersenneTwister(1234)
min_f_b, max_f_b = 1e-3, 1e3
function prior_transform(cube)
α₁ = quantile(Uniform(0.0, 1.25), cube[1])
f₁ = quantile(LogUniform(min_f_b, max_f_b), cube[2])
α₂ = quantile(Uniform(α₁, 4.0), cube[3])
norm = quantile(LogNormal(-1,2), cube[4])
return [α₁, f₁, α₂, norm]
end
P = 2000
unif = rand(rng, 4, P)
priors = mapreduce(permutedims, hcat, [prior_transform(unif[:, i]) for i in 1:P]')
l = @layout [a b ; c d]
p1 = histogram(priors[1,:],xlabel="α₁")
bins = 10.0 .^LinRange( log10(minimum(priors[2,:])),log10(quantile(priors[2,:],.99)),30)
p2 = histogram(priors[2,:],bins=bins,xaxis=(:log10,(bins[1],bins[end])),xlabel="f₁")
p3 = histogram(priors[3,:],xlabel="α₂")
bins = 10.0 .^LinRange( log10(minimum(priors[4,:])),log10(quantile(priors[4,:],1)),30)
p4 = histogram(priors[4,:],xlabel="norm",bins=bins,xaxis=(:log10,(bins[1],bins[end])))
plot(p1,p2,p3,p4,layout=l,size=(700,300),grid=false,left_margin=2Plots.mm,bottom_margin=20Plots.px,legend=false)
We can then use the function run_diagnostics to assess the quality of the approximation. The first argument is an array containing the parameters of the power spectral density, the second argument is the normalisation of the PSD of the process. f_min and f_max are the minimum and maximum frequencies of the time series, this is to show the window of observed frequencies in the plots.
using CairoMakie
CairoMakie.activate!(type = "png")
f_min, f_max = 1e-3 * 5, 1e3 / 5
S_low = S_high = 20
figs = run_diagnostics(priors[1:3, :], priors[4, :], f_min,f_max, SingleBendingPowerLaw, S_low,S_high, n_components=20, basis_function="SHO")(Makie.Figure[Scene (800px, 600px):
0 Plots
2 Child Scenes:
├ Scene (800px, 600px)
└ Scene (800px, 600px), Scene (800px, 600px):
0 Plots
3 Child Scenes:
├ Scene (800px, 600px)
├ Scene (800px, 600px)
└ Scene (800px, 600px), Scene (800px, 600px):
0 Plots
2 Child Scenes:
├ Scene (800px, 600px)
└ Scene (800px, 600px)], [0.30106248670934643 0.1278649888836782 … 3.5379278340998606 0.8659743643611062; 0.2973929534295234 0.1257745098299968 … 3.4792221135782238 0.8487424095569791; … ; 3.0248383150405235e-10 1.0567601383209985e-9 … 3.5413819671715603e-17 4.406571347697121e-11; 2.8977531055655343e-10 1.0303059450743374e-9 … 3.3842568858022605e-17 4.249996007083752e-11], [0.3010624867093463 0.12786498888367823 … 3.5379278340998606 0.8659743643611063; 0.29745503354466174 0.1258698016124998 … 3.482056090845377 0.850007323623472; … ; 3.0321348450309807e-10 1.0681387194608287e-9 … 3.5474804214016795e-17 4.42807740374178e-11; 2.8977531055655353e-10 1.0303059450743374e-9 … 3.384256885802262e-17 4.249996007083752e-11], [0.00024999999999999995, 0.0002541858320725946, 0.0002584417489057492, 0.00026276892395161777, 0.0002671685503098493, 0.0002716418410565517, 0.0002761900295787644, 0.0002808143699145311, 0.00028551613709866595, 0.0002902966275143097 … 3444.75238504349, 3502.429005105338, 3561.0713237515615, 3620.695509990588, 3681.3180035538608, 3742.955519428645, 3805.625052466724, 3869.3438820702645, 3934.1295769561357, 4000.000000000001])The following plots are produced: The mean of the residuals and ratios as a function of frequency.
3-element Vector{Makie.Figure}:
Figure()
Figure()
Figure()The quantiles of the residuals and ratios as a function of frequency.
1000×2000 adjoint(::Matrix{Float64}) with eltype Float64:
0.301062 0.127865 0.00374528 … 3.53793 0.865974
0.297393 0.125775 0.00372614 3.47922 0.848742
0.293768 0.123718 0.0037071 3.42148 0.831853
0.290188 0.121696 0.00368815 3.3647 0.8153
0.286651 0.119706 0.00366931 3.30885 0.799077
0.283157 0.117749 0.00365056 … 3.25392 0.783176
0.279705 0.115824 0.0036319 3.1999 0.767592
0.276296 0.11393 0.00361334 3.14676 0.752317
0.272928 0.112067 0.00359488 3.0945 0.737347
0.269602 0.110235 0.0035765 3.04311 0.722675
⋮ ⋱
4.08493e-10 1.2616e-9 1.05592e-12 4.86563e-17 5.67544e-11
3.91332e-10 1.23011e-9 1.00956e-12 4.64975e-17 5.47406e-11
3.74891e-10 1.19939e-9 9.65235e-13 4.44345e-17 5.27978e-11
3.59141e-10 1.16943e-9 9.22858e-13 4.2463e-17 5.09236e-11
3.44052e-10 1.1402e-9 8.82341e-13 … 4.0579e-17 4.91156e-11
3.29597e-10 1.11169e-9 8.43602e-13 3.87786e-17 4.73714e-11
3.1575e-10 1.08388e-9 8.06565e-13 3.7058e-17 4.56888e-11
3.02484e-10 1.05676e-9 7.71154e-13 3.54138e-17 4.40657e-11
2.89775e-10 1.03031e-9 7.37297e-13 3.38426e-17 4.25e-11The distribution of the mean, median and maximum values of the frequency-averaged residuals and ratios.
1000×2000 adjoint(::Matrix{Float64}) with eltype Float64:
0.301062 0.127865 0.00374528 … 3.53793 0.865974
0.297455 0.12587 0.00372505 3.48206 0.850007
0.293843 0.123874 0.00370476 3.42615 0.83404
0.290233 0.12188 0.00368446 3.37033 0.818105
0.286631 0.119893 0.00366418 3.31469 0.802233
0.283045 0.117917 0.00364396 … 3.25934 0.786454
0.279479 0.115954 0.00362382 3.20437 0.770796
0.27594 0.114008 0.00360379 3.14987 0.755288
0.272433 0.112082 0.0035839 3.09594 0.739954
0.268963 0.110179 0.00356418 3.04265 0.724818
⋮ ⋱
4.11932e-10 1.34224e-9 1.0616e-12 4.88857e-17 5.80715e-11
3.94668e-10 1.3023e-9 1.0153e-12 4.67328e-17 5.5951e-11
3.78026e-10 1.26258e-9 9.70808e-13 4.46655e-17 5.38832e-11
3.61977e-10 1.2231e-9 9.28034e-13 4.26795e-17 5.18664e-11
3.46495e-10 1.18389e-9 8.86896e-13 … 4.07709e-17 4.98989e-11
3.31554e-10 1.14497e-9 8.47319e-13 3.89361e-17 4.79795e-11
3.17134e-10 1.10637e-9 8.09235e-13 3.71717e-17 4.61071e-11
3.03213e-10 1.06814e-9 7.7258e-13 3.54748e-17 4.42808e-11
2.89775e-10 1.03031e-9 7.37297e-13 3.38426e-17 4.25e-11Using these three diagnostics we can see that a SingleBendingPowerLaw with the chosen prior distributions can be well approximated with 20 SHO basis functions.
Building the Gaussian process
Now that we have checked that approximation hold for our choice of priors we can build the Gaussian process.
Building the covariance function
The covariance function using the approximation of the power spectral density is obtained using the function approx. We need to specify the frequencies of the observations $f_\mathrm{min}$ and $f_\mathrm{max}$ and the scaling factors $S_\mathrm{low}$ and $S_\mathrm{high}$ extending the grid of frequencies over which the approximation is performed; Following the notations presented above: $f0=f_\mathrm{min}/S_\mathrm{low}$ and $fM = f_\mathrm{max}S_\mathrm{high}$, by default $S_\mathrm{low}=S_\mathrm{high}=20$. We also need to set the number of basis functions to use, and the normalisation of the power spectrum. One can also give the type of basis function to use, the default is SHO which corresponds to the basis function $\psi_4$, DRWCelerite corresponds to $\psi_6$.
𝓟 = SingleBendingPowerLaw(.4, 1e-1, 3.)
normalisation = 2.2
𝓡 = approx(𝓟, f_min, f_max, 20, normalisation, S_low, S_high, basis_function="SHO",is_integrated_power=true)Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.035006437502329235, 0.035006437502329235, 0.0011107207345395916, 0.0011107207345395916), Celerite(0.04492967946989494, 0.04492967946989494, 0.0026593107478368213, 0.0026593107478368213), Celerite(0.09895239548841843, 0.09895239548841843, 0.006366977255080993, 0.006366977255080993), Celerite(0.11268990937898554, 0.11268990937898554, 0.015243949733852687, 0.015243949733852687), Celerite(0.3216349779788861, 0.3216349779788861, 0.03649738238074974, 0.03649738238074974), Celerite(0.24807123977224216, 0.24807123977224216, 0.08738279408574272, 0.08738279408574272), Celerite(1.4789472424502454, 1.4789472424502454, 0.20921370805646405, 0.20921370805646405), Celerite(2.3219807349474735, 2.3219807349474735, 0.5009038231918576, 0.5009038231918576), Celerite(0.4186781892965387, 0.4186781892965387, 1.1992743803408137, 1.1992743803408137), Celerite(0.05713094740707749, 0.05713094740707749, 2.87132773348978, 2.87132773348978), Celerite(0.010622782994912401, 0.010622782994912401, 6.874592743959563, 6.874592743959563), Celerite(0.001796488667512047, 0.001796488667512047, 16.4592933241592, 16.4592933241592), Celerite(0.00031736662641762477, 0.00031736662641762477, 39.40718335187898, 39.40718335187898), Celerite(5.507263749964364e-5, 5.507263749964364e-5, 94.34950025765693, 94.34950025765693), Celerite(9.62865926712638e-6, 9.62865926712638e-6, 225.89354127095058, 225.89354127095058), Celerite(1.6781828357034303e-6, 1.6781828357034303e-6, 540.8390277487408, 540.8390277487408), Celerite(2.9285172131350254e-7, 2.9285172131350254e-7, 1294.8880799799056, 1294.8880799799056), Celerite(5.112875978017468e-8, 5.112875978017468e-8, 3100.2480472859197, 3100.2480472859197), Celerite(8.803006707828108e-9, 8.803006707828108e-9, 7422.678533614501, 7422.678533614501), Celerite(1.8072802445520616e-9, 1.8072802445520616e-9, 17771.531752633466, 17771.531752633466)], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466])info "Normalisation of the power spectrum"
As can be observed in the code, we also specify the optical argument is_integrated_power which tells if the normalisation corresponds to the integrated power between $f_\mathrm{min}$ and $f_\mathrm{max}$. If not, the normalisation corresponds to $\mathcal{R}(0)$ the variance of the process which is the integral of the power spectrum between $0$ and $+\infty$. Both integrals can be computed analytically as described in: Integral of the basis functions.
Building the Gaussian process
The Gaussian process is built using the type ScalableGP. If the mean of the process $\mu$ is known, it can be given as a first argument. Otherwise, the mean is assumed to be zero.
μ = 1.3
f = ScalableGP(μ, 𝓡)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.035006437502329235, 0.035006437502329235, 0.0011107207345395916, 0.0011107207345395916), Celerite(0.04492967946989494, 0.04492967946989494, 0.0026593107478368213, 0.0026593107478368213), Celerite(0.09895239548841843, 0.09895239548841843, 0.006366977255080993, 0.006366977255080993), Celerite(0.11268990937898554, 0.11268990937898554, 0.015243949733852687, 0.015243949733852687), Celerite(0.3216349779788861, 0.3216349779788861, 0.03649738238074974, 0.03649738238074974), Celerite(0.24807123977224216, 0.24807123977224216, 0.08738279408574272, 0.08738279408574272), Celerite(1.4789472424502454, 1.4789472424502454, 0.20921370805646405, 0.20921370805646405), Celerite(2.3219807349474735, 2.3219807349474735, 0.5009038231918576, 0.5009038231918576), Celerite(0.4186781892965387, 0.4186781892965387, 1.1992743803408137, 1.1992743803408137), Celerite(0.05713094740707749, 0.05713094740707749, 2.87132773348978, 2.87132773348978), Celerite(0.010622782994912401, 0.010622782994912401, 6.874592743959563, 6.874592743959563), Celerite(0.001796488667512047, 0.001796488667512047, 16.4592933241592, 16.4592933241592), Celerite(0.00031736662641762477, 0.00031736662641762477, 39.40718335187898, 39.40718335187898), Celerite(5.507263749964364e-5, 5.507263749964364e-5, 94.34950025765693, 94.34950025765693), Celerite(9.62865926712638e-6, 9.62865926712638e-6, 225.89354127095058, 225.89354127095058), Celerite(1.6781828357034303e-6, 1.6781828357034303e-6, 540.8390277487408, 540.8390277487408), Celerite(2.9285172131350254e-7, 2.9285172131350254e-7, 1294.8880799799056, 1294.8880799799056), Celerite(5.112875978017468e-8, 5.112875978017468e-8, 3100.2480472859197, 3100.2480472859197), Celerite(8.803006707828108e-9, 8.803006707828108e-9, 7422.678533614501, 7422.678533614501), Celerite(1.8072802445520616e-9, 1.8072802445520616e-9, 17771.531752633466, 17771.531752633466)], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466])), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.035006437502329235, 0.035006437502329235, 0.0011107207345395916, 0.0011107207345395916), Celerite(0.04492967946989494, 0.04492967946989494, 0.0026593107478368213, 0.0026593107478368213), Celerite(0.09895239548841843, 0.09895239548841843, 0.006366977255080993, 0.006366977255080993), Celerite(0.11268990937898554, 0.11268990937898554, 0.015243949733852687, 0.015243949733852687), Celerite(0.3216349779788861, 0.3216349779788861, 0.03649738238074974, 0.03649738238074974), Celerite(0.24807123977224216, 0.24807123977224216, 0.08738279408574272, 0.08738279408574272), Celerite(1.4789472424502454, 1.4789472424502454, 0.20921370805646405, 0.20921370805646405), Celerite(2.3219807349474735, 2.3219807349474735, 0.5009038231918576, 0.5009038231918576), Celerite(0.4186781892965387, 0.4186781892965387, 1.1992743803408137, 1.1992743803408137), Celerite(0.05713094740707749, 0.05713094740707749, 2.87132773348978, 2.87132773348978), Celerite(0.010622782994912401, 0.010622782994912401, 6.874592743959563, 6.874592743959563), Celerite(0.001796488667512047, 0.001796488667512047, 16.4592933241592, 16.4592933241592), Celerite(0.00031736662641762477, 0.00031736662641762477, 39.40718335187898, 39.40718335187898), Celerite(5.507263749964364e-5, 5.507263749964364e-5, 94.34950025765693, 94.34950025765693), Celerite(9.62865926712638e-6, 9.62865926712638e-6, 225.89354127095058, 225.89354127095058), Celerite(1.6781828357034303e-6, 1.6781828357034303e-6, 540.8390277487408, 540.8390277487408), Celerite(2.9285172131350254e-7, 2.9285172131350254e-7, 1294.8880799799056, 1294.8880799799056), Celerite(5.112875978017468e-8, 5.112875978017468e-8, 3100.2480472859197, 3100.2480472859197), Celerite(8.803006707828108e-9, 8.803006707828108e-9, 7422.678533614501, 7422.678533614501), Celerite(1.8072802445520616e-9, 1.8072802445520616e-9, 17771.531752633466, 17771.531752633466)], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466]))At the moment, the GP does not include the measurement variance σ² and the time values t. This is done in the next step.
σ² = yerr .^ 2
fx = f(t, σ²)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.035006437502329235, 0.035006437502329235, 0.0011107207345395916, 0.0011107207345395916), Celerite(0.04492967946989494, 0.04492967946989494, 0.0026593107478368213, 0.0026593107478368213), Celerite(0.09895239548841843, 0.09895239548841843, 0.006366977255080993, 0.006366977255080993), Celerite(0.11268990937898554, 0.11268990937898554, 0.015243949733852687, 0.015243949733852687), Celerite(0.3216349779788861, 0.3216349779788861, 0.03649738238074974, 0.03649738238074974), Celerite(0.24807123977224216, 0.24807123977224216, 0.08738279408574272, 0.08738279408574272), Celerite(1.4789472424502454, 1.4789472424502454, 0.20921370805646405, 0.20921370805646405), Celerite(2.3219807349474735, 2.3219807349474735, 0.5009038231918576, 0.5009038231918576), Celerite(0.4186781892965387, 0.4186781892965387, 1.1992743803408137, 1.1992743803408137), Celerite(0.05713094740707749, 0.05713094740707749, 2.87132773348978, 2.87132773348978), Celerite(0.010622782994912401, 0.010622782994912401, 6.874592743959563, 6.874592743959563), Celerite(0.001796488667512047, 0.001796488667512047, 16.4592933241592, 16.4592933241592), Celerite(0.00031736662641762477, 0.00031736662641762477, 39.40718335187898, 39.40718335187898), Celerite(5.507263749964364e-5, 5.507263749964364e-5, 94.34950025765693, 94.34950025765693), Celerite(9.62865926712638e-6, 9.62865926712638e-6, 225.89354127095058, 225.89354127095058), Celerite(1.6781828357034303e-6, 1.6781828357034303e-6, 540.8390277487408, 540.8390277487408), Celerite(2.9285172131350254e-7, 2.9285172131350254e-7, 1294.8880799799056, 1294.8880799799056), Celerite(5.112875978017468e-8, 5.112875978017468e-8, 3100.2480472859197, 3100.2480472859197), Celerite(8.803006707828108e-9, 8.803006707828108e-9, 7422.678533614501, 7422.678533614501), Celerite(1.8072802445520616e-9, 1.8072802445520616e-9, 17771.531752633466, 17771.531752633466)], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466])), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.035006437502329235, 0.035006437502329235, 0.0011107207345395916, 0.0011107207345395916), Celerite(0.04492967946989494, 0.04492967946989494, 0.0026593107478368213, 0.0026593107478368213), Celerite(0.09895239548841843, 0.09895239548841843, 0.006366977255080993, 0.006366977255080993), Celerite(0.11268990937898554, 0.11268990937898554, 0.015243949733852687, 0.015243949733852687), Celerite(0.3216349779788861, 0.3216349779788861, 0.03649738238074974, 0.03649738238074974), Celerite(0.24807123977224216, 0.24807123977224216, 0.08738279408574272, 0.08738279408574272), Celerite(1.4789472424502454, 1.4789472424502454, 0.20921370805646405, 0.20921370805646405), Celerite(2.3219807349474735, 2.3219807349474735, 0.5009038231918576, 0.5009038231918576), Celerite(0.4186781892965387, 0.4186781892965387, 1.1992743803408137, 1.1992743803408137), Celerite(0.05713094740707749, 0.05713094740707749, 2.87132773348978, 2.87132773348978), Celerite(0.010622782994912401, 0.010622782994912401, 6.874592743959563, 6.874592743959563), Celerite(0.001796488667512047, 0.001796488667512047, 16.4592933241592, 16.4592933241592), Celerite(0.00031736662641762477, 0.00031736662641762477, 39.40718335187898, 39.40718335187898), Celerite(5.507263749964364e-5, 5.507263749964364e-5, 94.34950025765693, 94.34950025765693), Celerite(9.62865926712638e-6, 9.62865926712638e-6, 225.89354127095058, 225.89354127095058), Celerite(1.6781828357034303e-6, 1.6781828357034303e-6, 540.8390277487408, 540.8390277487408), Celerite(2.9285172131350254e-7, 2.9285172131350254e-7, 1294.8880799799056, 1294.8880799799056), Celerite(5.112875978017468e-8, 5.112875978017468e-8, 3100.2480472859197, 3100.2480472859197), Celerite(8.803006707828108e-9, 8.803006707828108e-9, 7422.678533614501, 7422.678533614501), Celerite(1.8072802445520616e-9, 1.8072802445520616e-9, 17771.531752633466, 17771.531752633466)], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.035006437502329235, 0.04492967946989494, 0.09895239548841843, 0.11268990937898554, 0.3216349779788861, 0.24807123977224216, 1.4789472424502454, 2.3219807349474735, 0.4186781892965387, 0.05713094740707749, 0.010622782994912401, 0.001796488667512047, 0.00031736662641762477, 5.507263749964364e-5, 9.62865926712638e-6, 1.6781828357034303e-6, 2.9285172131350254e-7, 5.112875978017468e-8, 8.803006707828108e-9, 1.8072802445520616e-9], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466], [0.0011107207345395916, 0.0026593107478368213, 0.006366977255080993, 0.015243949733852687, 0.03649738238074974, 0.08738279408574272, 0.20921370805646405, 0.5009038231918576, 1.1992743803408137, 2.87132773348978, 6.874592743959563, 16.4592933241592, 39.40718335187898, 94.34950025765693, 225.89354127095058, 540.8390277487408, 1294.8880799799056, 3100.2480472859197, 7422.678533614501, 17771.531752633466]))
x: [0.0, 5.5, 7.25, 7.5, 7.75, 10.25, 10.5, 13.5, 16.5, 20.0 … 721.0, 721.25, 724.0, 725.0, 725.25, 725.75, 726.75, 727.75, 728.5, 729.5]
Σy: [0.0007769427202587575 0.0 … 0.0 0.0; 0.0 0.0008096922961612274 … 0.0 0.0; … ; 0.0 0.0 … 0.0007682230802508954 0.0; 0.0 0.0 … 0.0 0.0006084949733623623]
)
The log-likelihood of the Gaussian process given the data y can be computed using the function logpdf from the Distributions package.
logpdf(fx, y)-450.2383788400836