Using custom mean function

In some cases, we want to use a mean function that is not constant over time. For instance, this can be useful when one wants to model a periodic signal embedded in red noise.

Let's say we want to model a mean function of the form:

\[m(t ; A, T_0, \phi,\mu) = A \sin{(2\pi t / T_0 +\phi)} +\mu\]

and the broadband noise is modelled with a power-law power spectrum, typical for the variability of active galaxies.

Implementation

using Plots
using Pioran
using Random
using Distributions

To use a custom mean function, first define the function like:

mean_function(x, A=A, ϕ=ϕ, T₀=T₀, μ=μ) = @. A * sin(2π * x / T₀ + ϕ) + μ

mean_function (generic function with 5 methods)

Then, create a CustomMean struct with the function as an argument:

μ_fun = CustomMean(mean_function)

CustomMean{typeof(Main.mean_function)}(Main.mean_function)

And use it in the ScalableGP struct like any mean value and that's it! This whole process can be summarised in a function as follows:

function GP_model(t, y, σ, params)

    α₁, f₁, α₂, variance, ν, μ, A, ϕ, T₀ = params

    T = (t[end] - t[1]) # duration of the time series
    Δt = minimum(diff(t)) # min time separation
    f_min, f_max = 1 / T, 1 / Δt / 2

    # Rescale the measurement variance
    σ² = ν .* σ .^ 2

    # Define the mean
    mean_function(x, A=A, ϕ=ϕ, T₀=T₀, μ=μ) = @. A * sin(2π * x / T₀ + ϕ) + μ
    μ_fun = CustomMean(mean_function)

    # Define the power spectral density function
    𝓟 = SingleBendingPowerLaw(α₁, f₁, α₂)

    # Approximate  the PSD to form a covariance function
    𝓡 = approx(𝓟, f_min, f_max, 20, variance, basis_function="SHO")

    # Build the GP using the mean function and (auto)covariance function
    f = ScalableGP(μ_fun, 𝓡)

    # Condition on the times and errors of the observation
    fx = f(t, σ²)
    return fx
end

GP_model (generic function with 1 method)

Sampling from the GP

We first define the GP to sample from:

t = LinRange(0,2000,200)
σ = 0.5 * ones(length(t))

params = [0.3,1e-2,2.9,1.03,1.0,0.2,2.3,0.3,320]
GP = GP_model(t, nothing, σ, params)

AbstractGPs.FiniteGP{ScalableGP{AbstractGPs.GP{CustomMean{Main.var"#mean_function#1"{Float64, Float64, Float64, Float64}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, LinRange{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}(
f: ScalableGP{AbstractGPs.GP{CustomMean{Main.var"#mean_function#1"{Float64, Float64, Float64, Float64}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(AbstractGPs.GP{CustomMean{Main.var"#mean_function#1"{Float64, Float64, Float64, Float64}}, Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}}(CustomMean{Main.var"#mean_function#1"{Float64, Float64, Float64, Float64}}(Main.var"#mean_function#1"{Float64, Float64, Float64, Float64}(2.3, 0.3, 320.0, 0.2)), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.005679255783333235, 0.005679255783333235, 0.00011107207345395915, 0.00011107207345395915), Celerite(0.005056818936868288, 0.005056818936868288, 0.0001939563435589432, 0.0001939563435589432), Celerite(0.010016973412992233, 0.010016973412992233, 0.0003386905640358681, 0.0003386905640358681), Celerite(0.011051151646984954, 0.011051151646984954, 0.0005914284424117007, 0.0005914284424117007), Celerite(0.022740132105458825, 0.022740132105458825, 0.0010327645338725382, 0.0010327645338725382), Celerite(0.022446339643001457, 0.022446339643001457, 0.001803434711519142, 0.001803434711519142), Celerite(0.05279759973104132, 0.05279759973104132, 0.0031491948571440134, 0.0031491948571440134), Celerite(0.04481525644925057, 0.04481525644925057, 0.005499188955894198, 0.005499188955894198), Celerite(0.13122035327935017, 0.13122035327935017, 0.009602797078124968, 0.009602797078124968), Celerite(0.1260602178395168, 0.1260602178395168, 0.016768602145377818, 0.016768602145377818), Celerite(0.6217736149319335, 0.6217736149319335, 0.029281678621586957, 0.029281678621586957), Celerite(0.8770572162926112, 0.8770572162926112, 0.05113227062485024, 0.05113227062485024), Celerite(0.2882612435546154, 0.2882612435546154, 0.0892882246622794, 0.0892882246622794), Celerite(0.08219135777553364, 0.08219135777553364, 0.15591693789297742, 0.15591693789297742), Celerite(0.028538171673376533, 0.028538171673376533, 0.2722653699731647, 0.2722653699731647), Celerite(0.009632907700215713, 0.009632907700215713, 0.4754353996966427, 0.4754353996966427), Celerite(0.0033491594631901243, 0.0033491594631901243, 0.8302150923820582, 0.8302150923820582), Celerite(0.001190259180458787, 0.001190259180458787, 1.4497387027948254, 1.4497387027948254), Celerite(0.00035134464676607236, 0.00035134464676607236, 2.5315635979958984, 2.5315635979958984), Celerite(0.0002191600720433031, 0.0002191600720433031, 4.420668523467678, 4.420668523467678)], [0.005679255783333235, 0.005056818936868288, 0.010016973412992233, 0.011051151646984954, 0.022740132105458825, 0.022446339643001457, 0.05279759973104132, 0.04481525644925057, 0.13122035327935017, 0.1260602178395168, 0.6217736149319335, 0.8770572162926112, 0.2882612435546154, 0.08219135777553364, 0.028538171673376533, 0.009632907700215713, 0.0033491594631901243, 0.001190259180458787, 0.00035134464676607236, 0.0002191600720433031], [0.005679255783333235, 0.005056818936868288, 0.010016973412992233, 0.011051151646984954, 0.022740132105458825, 0.022446339643001457, 0.05279759973104132, 0.04481525644925057, 0.13122035327935017, 0.1260602178395168, 0.6217736149319335, 0.8770572162926112, 0.2882612435546154, 0.08219135777553364, 0.028538171673376533, 0.009632907700215713, 0.0033491594631901243, 0.001190259180458787, 0.00035134464676607236, 0.0002191600720433031], [0.00011107207345395915, 0.0001939563435589432, 0.0003386905640358681, 0.0005914284424117007, 0.0010327645338725382, 0.001803434711519142, 0.0031491948571440134, 0.005499188955894198, 0.009602797078124968, 0.016768602145377818, 0.029281678621586957, 0.05113227062485024, 0.0892882246622794, 0.15591693789297742, 0.2722653699731647, 0.4754353996966427, 0.8302150923820582, 1.4497387027948254, 2.5315635979958984, 4.420668523467678], [0.00011107207345395915, 0.0001939563435589432, 0.0003386905640358681, 0.0005914284424117007, 0.0010327645338725382, 0.001803434711519142, 0.0031491948571440134, 0.005499188955894198, 0.009602797078124968, 0.016768602145377818, 0.029281678621586957, 0.05113227062485024, 0.0892882246622794, 0.15591693789297742, 0.2722653699731647, 0.4754353996966427, 0.8302150923820582, 1.4497387027948254, 2.5315635979958984, 4.420668523467678])), Pioran.SumOfCelerite{StructArrays.StructArray{<:Pioran.SemiSeparable}}(Celerite[Celerite(0.005679255783333235, 0.005679255783333235, 0.00011107207345395915, 0.00011107207345395915), Celerite(0.005056818936868288, 0.005056818936868288, 0.0001939563435589432, 0.0001939563435589432), Celerite(0.010016973412992233, 0.010016973412992233, 0.0003386905640358681, 0.0003386905640358681), Celerite(0.011051151646984954, 0.011051151646984954, 0.0005914284424117007, 0.0005914284424117007), Celerite(0.022740132105458825, 0.022740132105458825, 0.0010327645338725382, 0.0010327645338725382), Celerite(0.022446339643001457, 0.022446339643001457, 0.001803434711519142, 0.001803434711519142), Celerite(0.05279759973104132, 0.05279759973104132, 0.0031491948571440134, 0.0031491948571440134), Celerite(0.04481525644925057, 0.04481525644925057, 0.005499188955894198, 0.005499188955894198), Celerite(0.13122035327935017, 0.13122035327935017, 0.009602797078124968, 0.009602797078124968), Celerite(0.1260602178395168, 0.1260602178395168, 0.016768602145377818, 0.016768602145377818), Celerite(0.6217736149319335, 0.6217736149319335, 0.029281678621586957, 0.029281678621586957), Celerite(0.8770572162926112, 0.8770572162926112, 0.05113227062485024, 0.05113227062485024), Celerite(0.2882612435546154, 0.2882612435546154, 0.0892882246622794, 0.0892882246622794), Celerite(0.08219135777553364, 0.08219135777553364, 0.15591693789297742, 0.15591693789297742), Celerite(0.028538171673376533, 0.028538171673376533, 0.2722653699731647, 0.2722653699731647), Celerite(0.009632907700215713, 0.009632907700215713, 0.4754353996966427, 0.4754353996966427), Celerite(0.0033491594631901243, 0.0033491594631901243, 0.8302150923820582, 0.8302150923820582), Celerite(0.001190259180458787, 0.001190259180458787, 1.4497387027948254, 1.4497387027948254), Celerite(0.00035134464676607236, 0.00035134464676607236, 2.5315635979958984, 2.5315635979958984), Celerite(0.0002191600720433031, 0.0002191600720433031, 4.420668523467678, 4.420668523467678)], [0.005679255783333235, 0.005056818936868288, 0.010016973412992233, 0.011051151646984954, 0.022740132105458825, 0.022446339643001457, 0.05279759973104132, 0.04481525644925057, 0.13122035327935017, 0.1260602178395168, 0.6217736149319335, 0.8770572162926112, 0.2882612435546154, 0.08219135777553364, 0.028538171673376533, 0.009632907700215713, 0.0033491594631901243, 0.001190259180458787, 0.00035134464676607236, 0.0002191600720433031], [0.005679255783333235, 0.005056818936868288, 0.010016973412992233, 0.011051151646984954, 0.022740132105458825, 0.022446339643001457, 0.05279759973104132, 0.04481525644925057, 0.13122035327935017, 0.1260602178395168, 0.6217736149319335, 0.8770572162926112, 0.2882612435546154, 0.08219135777553364, 0.028538171673376533, 0.009632907700215713, 0.0033491594631901243, 0.001190259180458787, 0.00035134464676607236, 0.0002191600720433031], [0.00011107207345395915, 0.0001939563435589432, 0.0003386905640358681, 0.0005914284424117007, 0.0010327645338725382, 0.001803434711519142, 0.0031491948571440134, 0.005499188955894198, 0.009602797078124968, 0.016768602145377818, 0.029281678621586957, 0.05113227062485024, 0.0892882246622794, 0.15591693789297742, 0.2722653699731647, 0.4754353996966427, 0.8302150923820582, 1.4497387027948254, 2.5315635979958984, 4.420668523467678], [0.00011107207345395915, 0.0001939563435589432, 0.0003386905640358681, 0.0005914284424117007, 0.0010327645338725382, 0.001803434711519142, 0.0031491948571440134, 0.005499188955894198, 0.009602797078124968, 0.016768602145377818, 0.029281678621586957, 0.05113227062485024, 0.0892882246622794, 0.15591693789297742, 0.2722653699731647, 0.4754353996966427, 0.8302150923820582, 1.4497387027948254, 2.5315635979958984, 4.420668523467678]))
x: LinRange{Float64}(0.0, 2000.0, 200)
Σy: [0.25 0.0 … 0.0 0.0; 0.0 0.25 … 0.0 0.0; … ; 0.0 0.0 … 0.25 0.0; 0.0 0.0 … 0.0 0.25]
)

We now sample a few realisations from the GP and see that the realisations are indeed periodic.

rng = MersenneTwister(12)
y = [rand(rng,GP) for i in 1:3]
Plots.scatter(t,y,yerr=σ,xlabel="Time",ylabel="Value",legend=false,framestyle = :box,ms=3)

Inference

To use such process for inference it is as easy as before. You need to get the loglikelihood using the logpdf function as follows:

GP = GP_model(t,y, σ, params)
logpdf(GP,y[1])

-310.86672106653333