API Reference

approx(psd_model, f_min, f_max, n_components=20, norm=1.0,S_low::Real=20., S_high::Real=20. ; is_integrated_power::Bool = true, basis_function="SHO")

Approximate the PSD with a sum of basis functions to form a covariance function. The PSD model is approximated between f0=f_min/S_low and fM=f_min*S_high. By default it is normalised by its integral from f_min to f_max but it can also be normalised by its integral from 0 to infinity using the variance argument.

Arguments

• psd_model::PowerSpectralDensity: model of the PSD
• f_min::Real: the minimum frequency in the time series
• f_max::Real: the maximum frequency in the time series
• n_components::Integer=20: the number of basis functions to use
• norm::Real=1.0: normalisation of the PSD.
• S_low::Real=20.0: scaling factor for the lowest frequency in the approximation.
• S_high::Real=20.0: scaling factor for the highest frequency in the approximation.
• is_integrated_power::Bool = true: if the norm corresponds to integral of the PSD between f_min and f_max, if not it is the variance of the process, integral of the PSD from 0 to +inf.
• basis_function::String="SHO": the basis function to use, either "SHO" or "DRWCelerite"

Return

• covariance::SumOfSemiSeparable: the covariance function

Example

using Pioran
𝓟 = SingleBendingPowerLaw(1.0, 1.0, 2.0)
𝓡 = approx(𝓟, 1e-4, 1e-1, 30, 2.31,basis_function="SHO")

 get_covariance_from_psd(psd_features)

Get the covariance function from the PSD features

Abstract type for continuum power spectra, i.e. broadband noise power spectra with no features

 DoubleBendingPowerLaw(A, α₁, f₁, α₂, f₂, α₃)

Double bending power law model for the power spectral density

• A : the amplitude
• α₁: the first power law index
• f₁: the first bend frequency
• α₂: the second power law index
• f₂: the second bend frequency
• α₃: the third power law index

\[\mathcal{P}(f) = A\frac{(f/f₁)^{-α₁}}{1 + (f / f₁)^{α₂ - α₁}}\frac{1}{1 + (f / f₂)^{α₃ - α₂}}\]

Lorentzian(A, γ, f₀)

Lorentzian model for the power spectral density

• A: the amplitude
• γ: the width of the peak
• f₀: the central frequency

\[\mathcal{P}(f) = \frac{A}{4\pi^2 (f - f₀)^2 + γ^2}\]

 PowerLaw(α)

Power law model for the power spectral density

• A: the power law amplitude
• α: the power law index

\[\mathcal{P}(f) = A f^{-α}\]

Abstract type for power spectral densities

QPO(S₀, f₀, Q)

QPO model

• S₀: the amplitude at the peak
• f₀: the central frequency
• Q: quality factor

\[\mathcal{P}(f) = \frac{S_0 {f_0}^4 } {\left(f^ 2 -{f_0}^2\right)^ 2 + f^2 {f_0}^2 / Q^2 }\]

 SingleBendingPowerLaw(A, α₁, f₁, α₂)

Single bending power law model for the power spectral density

• A: the amplitude
• α₁: the first power law index
• f₁: the first bend frequency
• α₂: the second power law index

\[\mathcal{P}(f) = A \frac{(f/f₁)^{-α₁}}{1 + (f / f₁)^{α₂ - α₁}}\]

 approximated_psd(f, psd_model, f0, fM; n_components=20, norm=1.0, basis_function="SHO")

Return the approximated PSD. This is essentially to check that the model and the approximation are consistent.

Arguments

• f::AbstractVector{<:Real}: the frequencies at which to evaluate the PSD
• psd_model::PowerSpectralDensity: model of the PSD
• f0::Real: the lowest frequency
• fM::Real: the highest frequency
• n_components::Integer=20: the number of basis functions to use
• norm::Real=1.0: normalisation of the PSD
• basis_function::String="SHO": the basis function to use, either "SHO" or "DRWCelerite"
• individual::Bool=false: return the individual components

 build_approx(J, f0, fM; basis_function="SHO")

Prepare the approximation of a PSD

Arguments

• J::Integer: the number of basis functions
• f0::Real: the lowest frequency
• fM::Real: the highest frequency
• basis_function::String="SHO": the basis function to use, either "SHO" or "DRWCelerite"

Return

• spectral_points::Vector{Real}: the spectral points
• spectral_matrix::Matrix{Real}: the spectral matrix

 convert_feature(psd_feature)

Convert a PSD feature to a Celerite covariance function Only QPO is implemented

Arguments

• psd_feature::PowerSpectralDensity: the PSD feature

Return

• covariance::SemiSeparable: the covariance function

 get_approx_coefficients(psd_model, f0, fM; n_components=20, basis_function="SHO")

Get the coefficients of the approximated PSD

Arguments

• psd_model::PowerSpectralDensity: model of the PSD
• f0::Real: the lowest frequency
• fM::Real: the highest frequency
• n_components::Integer=20: the number of basis functions to use
• basis_function::String="SHO": the basis function to use, either "SHO" or "DRWCelerite"

Return

• amplitudes::Vector{Real}: the amplitudes of the basis functions

get_norm_psd(amplitudes,spectral_points,f_min,f_max,basis_function,is_integrated_power, cov_features=nothing)

Get the normalisation of the sum of basis functions.

Arguments

• amplitudes: amplitude of the basis function
• spectral_points: spectral points of the basis function
• f_min::Real: the minimum frequency in the time series
• f_max::Real: the maximum frequency in the time series
• basis_function::String="SHO": the basis function to use, either "SHO" or "DRWCelerite"
• is_integrated_power::Bool=true: if the norm corresponds to integral of the PSD between f_min and f_max or if it is the integral from 0 to infinity.
• cov_features: PSD features of the PSD, is nothing if there are no features.

 get_normalised_psd(psd_model::PowerSpectralDensity, spectral_points::AbstractVector{<:Real})

Get the PSD normalised at the lowest frequency

integral_celerite(a, b, c, d, x)

Computes the integral of the Celerite power spectrum:

integral_drwcelerite(a, c, x)

Computes the integral of the DRWCelerite basis function of amplitude a and width c for a given x.

The DRWCelerite basis function is defined as:

\[ \int \dfrac{a\, {d}x}{(x/c)^6+1} =\dfrac{ac}{3} \left[ \arctan{(x/c)} +\dfrac{\sqrt3}{4}\ln{\left(\dfrac{x^2+xc\sqrt3+c^2}{x^2-xc\sqrt3+c^2}\right)}+\dfrac{1}{2}\arctan{\left(\dfrac{x^2-c^2}{xc}\right)}\right]\]

integral_sho(a, c, x)

Computes the integral of the SHO basis function of amplitude a and width c for a given x.

This integral is obtained using Equation: 4.2.7.1.3 from the "Handbook of Mathematical Formulas and Integrals" 2009

\[ \int \dfrac{a\, {d}x}{(x/c)^4+1} =\dfrac{ac}{4\sqrt2} \left[\ln{\left(\dfrac{x^2+cx\sqrt2+c^2}{x^2-cx\sqrt2+c^2}\right)}+2\arctan{\left(\dfrac{cx\sqrt2}{c^2-x^2}\right)}\right]\]

integrate_basis_function(a,c,x₁,x₂,basis_function)

Computes the integral of the basis function between x₁ and x₂ for a given amplitude a and width c.

integrate_psd_feature(a, b, c, d, x₁, x₂)

Computes the integral of a celerite power spectral density with coefficients (a,b,c,d) between x₁ and x₂.

 psd_decomp(psd_normalised, spectral_matrix)

Get amplitudes of the basis functions by solving the linear system of the approximation

 separate_psd(psd::PowerSpectralDensity)

Separate the PSD into its ContinuumPowerSpectrum components and other components if it is a sum of PSDs

Arguments

• psd::PowerSpectralDensity: power spectral density or sum of PowerSpectralDensity objects

Return

• psd_continuum::Union{SumOfPowerSpectralDensity,PowerSpectralDensity,nothing}: continuum part of the psd
• psd_line::Union{PowerSpectralDensity,nothing,Vector{PowerSpectralDensity}}: non-continuum part of the psd

Exp(A, α)

Exponential covariance Function

• A: the amplitude of the covariance function
• α: the decay rate of the covariance function

\[k(τ) = A \exp(-α τ)\]

Example

Exp(1.0, 0.25)

evaluate(R::Exp, f)

This is the right formula but it disagrees with the Celerite implementation...

Evaluate the power spectral density at frequency f

 SHO(A, ω₀, Q)

Simple Harmonic Oscillator covariance Function

• A: the amplitude of the covariance function
• ω₀: the angular frequency of the simple harmonic oscillator
• Q: the quality factor of the simple harmonic oscillator

\[k(τ) = A \exp(-ω₀ τ / Q / 2) \left\{\begin{matrix} 2(1 + ω₀ τ) & Q = 1/2 \\ \cos(η ω₀ τ) + \frac{\sin(η ω₀ τ)}{2η Q} & Q < 1/2 \\ \cosh(η ω₀ τ) + \frac{\sinh(η ω₀ τ)}{2η Q} & Q \geq 1/2 \end{matrix}\right.\\ η = \sqrt{\left|1 - \frac{1}{4 Q^2}\right|}\]

See Foreman-Mackey et al. (2017) for more details.

 Celerite(a, b, c, d)

Celerite covariance Function

• a: the amplitude of the first term
• b: the amplitude of the second term
• c: the decay rate of the covariance function
• d: the period of the covariance function

\[k(τ) = \exp(-c τ) (a \cos(d τ) + b \sin(d τ))\]

See Foreman-Mackey et al. (2017) for more details.

evaluate(f, C::Celerite)

evaluate the power spectral density at frequency f

 +(::SemiSeparable, ::SemiSeparable)

Sum of semi-separable covariance functions

 celerite_coefs(covariance)

Get the celerite coefficients

ScalableGP(μ, 𝓡; solver=:celerite)

Scalable Gaussian Process for semi-separable covariance functions.

• μ : mean or mean function of the Gaussian Process
• 𝓡 : covariance function
• solver : indicates which solver to use for the likelihood computation.

A scalable Gaussian process has a covariance function formed of semi-separable kernels

Example

using Pioran
𝓟 = SingleBendingPowerLaw(1.0, 1.0, 2.0)
𝓡 = approx(𝓟, 1e-4, 1e-1, 30, 2.31,basis_function="SHO")
μ = 1.2

f = ScalableGP(𝓡) # zero-mean GP
f = ScalableGP(μ, 𝓡) # with mean μ

See Foreman-Mackey et al. (2017) for more details.

rand(rng::AbstractRNG, fp::PosteriorGP, N::Int=1)
rand(rng::AbstractRNG, fp::PosteriorGP, τ::AbstractVecOrMat{<:Real}, N::Int=1)
rand(fp::PosteriorGP, N::Int=1)

Sample N realisations from the posterior GP fp at the points τ.

rand(rng::AbstractRNG, f::ScalableGP)
rand(rng::AbstractRNG, f::ScalableGP, t::AbstractVecOrMat{<:Real})
rand(f::ScalableGP)
rand(f::ScalableGP, t::AbstractVecOrMat{<:Real})

Draw a realisation from the GP f at the points t.

logpdf(f::ScalableGP, Y::AbstractVecOrMat{<:Real})

Compute the log-likelihood of the data Y given the GP f.

_predict_cov(fp::PosteriorGP, τ::AbstractVecOrMat{<:Real})
Compute the posterior covariance of the GP at the points τ.

_predict_mean(fp::PosteriorGP, τ::AbstractVecOrMat{<:Real})

Compute the Posterior mean of the GP at the points τ.

posterior(f::ScalableGP, y::AbstractVecOrMat{<:Real})

Compute the posterior Gaussian process fp given the GP f and the data y.

cov(fp::PosteriorGP, τ::AbstractVecOrMat{<:Real})
cov(fp::PosteriorGP)

Compute the covariance of the posterior GP at the points τ.

mean(fp::PosteriorGP, τ::AbstractVecOrMat{<:Real})
mean(fp::PosteriorGP)

Compute the mean of the posterior GP at the points τ.

std(fp::PosteriorGP, τ::AbstractVecOrMat{<:Real})
std(fp::PosteriorGP)

Compute the standard deviation of the posterior GP at the points τ.

compute_nll(t, y, σ², a, b, c, d)

Compute the likelihood using the vectorised implementation in `get_values!`.

Still experimental.

get_values!(a, b, c, d, zp, U, V, P, D, t)

Compute the values of the matrices and vectors needed for the celerite algorithm. This is a vectorised version of the init_semi_separable! and solve_prec! functions. This function appears to be faster than the two previous functions when J > 16 but it also uses more memory.

More study of this implementation is needed.

log_likelihood(cov, τ, y, σ2; solver = :celerite)

Compute the log-likelihood of a semi-separable covariance function using the celerite algorithm.

Arguments

• cov::SumOfSemiSeparable or cov::CARMA or cov::SemiSeparable: the covariance function
• τ::Vector: the time points
• y::Vector: the data
• σ2::Vector: the measurement variances
• solver::Symbol: solver to use for the likelihood computation either :celerite for the classic celerite solver or :celerite_matrix for the solver using matrices algebra.

logl(a, b, c, d, τ, y, σ2)

Compute the log-likelihood of a GP with a semi-separable covariance function using the celerite algorithm.

Arguments

• a::Vector
• b::Vector
• c::Vector
• d::Vector
• τ::Vector: the time points
• y::Vector: the data
• σ2::Vector: the measurement variances

See Foreman-Mackey et al. (2017) for more details.

predict(cov, τ, t, y, σ2)

Compute the posterior mean of the GP at the points τ given the data y and time t.

Arguments

• cov::SumOfSemiSeparable or cov::CARMA or cov::SemiSeparable: the covariance function
• τ::Vector: the time points
• t::Vector: the data time points
• y::Vector: the data
• σ2::Vector: the measurement variances

simulate(rng, cov, τ, σ2)
simulate(cov, τ, σ2)

Draw a realisation from the GP with the covariance function cov at the points τ with the variances σ2.

Arguments

• rng::AbstractRNG: the random number generator
• cov::SumOfSemiSeparable or cov::CARMA or cov::SemiSeparable: the covariance function
• τ::Vector: the time points
• σ2::Vector: the measurement variances

solve_prec!(z, y, U, W, D, ϕ)

Forward and backward substitution of the celerite algorithm.

See Foreman-Mackey et al. (2017) for more details.

log_likelihood_direct(cov::KernelFunctions.SimpleKernel, t::Vector, y::Vector, σ²::Vector)

Compute the log-likelihood of the data Y given the GP f with the direct solver.

predict_direct(cov::KernelFunctions.SimpleKernel, τ::AbstractVector, t::AbstractVector, σ²::AbstractVector)

Compute the posterior covariance of the GP at the points τ given the times t and the noise variance σ².

predict_direct(cov::KernelFunctions.SimpleKernel, τ::AbstractVector, t::AbstractVector, y::AbstractVector, σ²::AbstractVector, with_covariance::Bool=false)

Compute the posterior mean of the GP at the points τ given the data (t, y) and the noise variance σ².

run_diagnostics(prior_samples, norm_samples, f_min, f_max, model, S_low=20.0, S_high=20.0; path="", basis_function="SHO", n_components=20)

Run the prior predictive checks for the model and the approximation of the PSD

Arguments

• prior_samples::Array{Float64, 2} : Model samples from the prior distribution
• norm_samples::Array{Float64, 1} : The samples of the normalisation of the PSD
• f_min::Float64 : The minimum frequency of the observed data
• f_max::Float64 : The maximum frequency of the observed data
• model::Function : The model
• S_low::Float64=20.0 : the scaling factor for the appoximation at low frequencies
• S_high::Float64=20.0 : the scaling factor for the appoximation at high frequencies
• path::String="" : The path to save the plots
• basis_function::String="SHO" : The basis function for the approximation of the PSD
• n_components::Int=20 : The number of components to use for the approximation of the PSD

run_posterior_predict_checks(samples, paramnames, paramnames_split, t, y, yerr, model, with_log_transform; S_low = 20, S_high = 20, is_integrated_power = true, plots = "all", n_samples = 100, path = "", basis_function = "SHO", n_frequencies = 1000, plot_f_P = false, n_components = 20)

Run the posterior predictive checks for the model and the approximation of the PSD

Arguments

• samples::Array{Float64, 2} : The samples from the posterior distribution
• paramnames::Array{String, 1} : The names of the parameters
• paramnames_split::Dict : Dictionary to map the name of the parameters to the components for instance: Dict( "psd" => ["α₁", "f₁", "α₂"], "norm" => "variance", "scaleerr" => "ν", "logtransform" => "c", "mean" => "μ" )
• t::Array{Float64, 1} : The time series
• y::Array{Float64, 1} : The values of the time series
• yerr::Array{Float64, 1} : The errors of the time series
• model::Function : The model or a string representing the model
• with_log_transform::Bool : If true, the flux is log-transformed
• S_low::Float64=20.0 : the scaling factor for the appoximation at low frequencies
• S_high::Float64=20.0 : the scaling factor for the appoximation at high frequencies
• plots::String or Array{String, 1} : The type of plots to make. It can be "all", "psd", "lsp", "timeseries" or a combination of them in an array
• n_samples::Int=200 : The number of samples to draw from the posterior predictive distribution
• path::String="all" : The path to save the plots
• basis_function::String="SHO" : The basis function for the approximation of the PSD
• is_integrated_power = true : if the norm corresponds to integral of the PSD between f_min and f_max or if it is the integral from 0 to infinity.
• n_frequencies::Int=1000 : The number of frequencies to use for the approximation of the PSD
• plot_f_P::Bool=false : If true, the plots are made in terms of f * PSD
• n_components::Int=20 : The number of components to use for the approximation of the PSD
• save_samples::Bool=false : Whether to save samples of the PSD, this can create a very large file

sample_approx_model(samples, norm_samples, f0, fM, model; n_frequencies=1_000, basis_function="SHO", n_components=20)

Check the approximation of the PSD by computing the residuals and the ratios of the PSD and the approximated PSD

Arguments

• samples::Array{Float64, n} : The model samples
• norm_samples::Array{Float64, 1} : The normalisation samples
• f0::Float64 : The minimum frequency for the approximation of the PSD
• fM::Float64 : The maximum frequency for the approximation of the PSD
• model::Function : The model
• n_frequencies::Int=1_000 : The number of frequencies to use for the approximation of the PSD
• basis_function::String="SHO" : The basis function for the approximation of the PSD
• n_components::Int=20 : The number of components to use for the approximation of the PSD
• is_integrated_power=true: Does the normalisation correspond to the integral of the PSD between two frequencies? If false, it corresponds to the true variance of the process.

Return

• psd::Array{Float64, 2} : The PSD
• psd_approx::Array{Float64, 2} : The approximated PSD
• residuals::Array{Float64, 2} : The residuals (psd-approx_psd)
• ratios::Array{Float64, 2} : The ratios (approx_psd/psd)

get_ppc_timeseries(samples, t, y, yerr, GP_model, with_log_transform; t_pred = nothing, n_samples = 1000)

Get the random posterior predictive time series from the model and samples

plot_boxplot_psd_approx(residuals, ratios; path="")

Plot the boxplot of the residuals and ratios for the PSD approximation

plot_lsp_ppc(samples, t, y, yerr, GP_model; plot_f_P=false, n_frequencies=1000, n_samples=1000, is_integrated_power = true, n_components=20, bin_fact=10, path="", basis_function="SHO",with_log_transform = true)

Plot the posterior predictive Lomb-Scargle periodogram of random time series from the model to compare with the one of the data.

Arguments

• samples::Array{Float64, n} : Posterior samples
• t::Array{Float64, 1} : Time series
• y::Array{Float64, 1} : The values of the time series
• yerr::Array{Float64, 1} : The errors of the time series
• GP_model::Function : The GP model is a function that returns a ScalableGP conditionned on t and yerr, its arguments are t,y,yerr,params
• S_low = 20.0: scaling factor for the lowest frequency in the approximation.
• S_high = 20.0: scaling factor for the highest frequency in the approximation.
• plot_f_P::Bool=false : If true, the plot is made in terms of f * PSD
• n_frequencies::Int=1000 : The number of frequencies to use for the approximation of the PSD
• n_samples::Int=1000 : The number of samples to draw from the posterior predictive distribution
• is_integrated_power::Bool = true : if the norm corresponds to integral of the PSD between f_min and f_max or if it is the integral from 0 to infinity.
• n_components::Int=20 : The number of components to use for the approximation of the PSD
• bin_fact::Int=10 : The binning factor for the LSP
• path::String="" : The path to save the plots
• basis_function::String="SHO" : The basis function for the approximation of the PSD
• with_log_transform = true : Apply a log transform to the data.

plot_mean_approx(f, residuals, ratios; path="")

Plot the frequency-averaged residuals and ratios

plot_ppc_timeseries(samples, t, y, yerr, GP_model, with_log_transform; t_pred = nothing, n_samples = 100, path = "")

Plot the posterior predictive time series and the residuals

Arguments

• samples::Matrix{Float64} : The samples of the model parameters
• t::Vector{Float64} : The time series
• y::Vector{Float64} : The values of the time series
• yerr::Vector{Float64} : The errors of the time series
• GP_model::PowerSpectralDensity : The model
• with_log_transform::Bool : If true, the flux was log-transformed for the inference
• t_pred::Vector{Float64}=nothing : The prediction times
• n_samples::Int=100 : The number of samples to draw from the posterior predictive distribution

plot_psd_ppc(samples_𝓟, samples_norm, samples_ν, t, y, yerr, model; S_low = 20.0, S_high = 20.0, plot_f_P = false, n_frequencies = 1000, path = "", n_components = 20, basis_function = "SHO", is_integrated_power = true, with_log_transform = false, save_samples = false)

Plot the posterior predictive power spectral density and the noise levels

Arguments

• samples_𝓟::Array{Float64, n} : The samples of the model parameters
• samples_norm::Array{Float64, 1} : The normalisation samples, either variance or integrated power between f_min or f_max.
• samples_ν::Array{Float64, 1} : The ν samples
• t::Array{Float64, 1} : The time series
• y::Array{Float64, 1} : The values of the time series
• yerr::Array{Float64, 1} : The errors of the time series
• model::Function : The model
• S_low = 20.0: scaling factor for the lowest frequency in the approximation.
• S_high = 20.0: scaling factor for the highest frequency in the approximation.
• plot_f_P::Bool = false : If true, the plot is made in terms of f * PSD
• n_frequencies::Int = 1000 : The number of frequencies to use for the approximation of the PSD
• path::String = "" : The path to save the plots
• n_components::Int = 20 : The number of components to use for the approximation of the PSD
• basis_function::String = "SHO" : The basis function for the approximation of the PSD
• is_integrated_power::Bool = true : if the norm corresponds to integral of the PSD between f_min and f_max or if it is the integral from 0 to infinity.
• with_log_transform::Bool = false : If true, the flux is log-transformed.
• save_samples::Bool = false : Save samples of the psd and approx

plot_quantiles_approx(f, f_min, f_max, residuals, ratios; path="")

Plot the quantiles of the residuals and ratios (with respect to the approximated PSD) of the PSD

plot_residuals_diagnostics(t, mean_res, res_quantiles; confidence_intervals=[95, 99], path="")

Plot the residuals and the autocorrelation function of the residuals

Arguments

• t::Vector{Float64} : The time series
• mean_res::Vector{Float64} : The mean of the residuals
• res_quantiles::Array{Float64, 2} : The quantiles of the residuals
• confidence_intervals::Array{Int, 1} : The confidence intervals
• path::String="" : The path to save the plot

plot_simu_ppc_timeseries(t_pred, ts_quantiles, t, y, yerr; path="")

Plot the posterior predictive simulated time series

Arguments

• t_pred::Vector{Float64} : The prediction times
• ts_quantiles::Array{Float64, 2} : The quantiles of the posterior predictive time series
• t::Vector{Float64} : The time series
• y::Vector{Float64} : The values of the time series
• yerr::Vector{Float64} : The errors of the time series
• path::String="" : The path to save the plot

check_conjugate_pair(r::Vector{Complex})

Check if the roots are complex conjugate pairs and negative real parts Returns true if the roots are complex conjugate pairs and false otherwise

check_order_imag_roots(r::Vector{Complex})

Check if the imaginary parts of the roots are in ascending order

check_roots_bounds(r::Vector{Complex},f_min::Float64,f_max::Float64)

Check if the roots are within the bounds of the frequency range

extract_subset(rng, prefix, t, y, yerr; n_perc=0.03, take_log=true, suffix="")
extract_subset(seed, prefix, t, y, yerr; n_perc=0.03, take_log=true)

Extract a subset of the data for the analysis and return initial guesses for the mean and variance. Either a random number generator or a seed can be provided.

Arguments

• seed::Int64 : Seed for the random number generator.
• rng::AbstractRNG : Random number generator.
• prefix::String : Prefix for the output files.
• t::Array{Float64,1} : Time array.
• y::Array{Float64,1} : Time series array.
• yerr::Array{Float64,1} : Time series error array.
• n_perc::Float64 : Percentage of the time series to extract.
• take_log::Bool : If true, log transform the time series for the estimation of the mean and variance.
• suffix::String : Suffix for the output files.

Return

• t_subset::Array{Float64,1} : Time array of the subset.
• y_subset::Array{Float64,1} : Time series array of the subset.
• yerr_subset::Array{Float64,1} : Time series error array of the subset.
• x̄::Float64 : Mean of the normal distribution for μ.
• va::Float64 : Variance of the normal distribution for μ.

simply put entries of paramnames_split into lists even when it's a string

separate_samples(samples, paramnames, with_log_transform::Bool)

Separate the samples into the parameters of the model and the parameters of the power spectral density.

ThreeUniformDependent(a, b, c, ϵ)
ThreeUniformDependent(a, b, c) (constructor with default ϵ = 1e-10)

Multivariate distribution to model three random variables where the first one x1 is given by U[a,b] and the second one x2 is given by U[x1,c] and the third one x3 is given by U[x2,c]. where a<b<c.

• a: lower bound of the first distribution
• b: upper bound of the first distribution
• c: upper bound of the second and third distribution
• ϵ: small value to make sure that the lower and upper bounds of each distribution are different

This means that the lower bound of the second distribution is dependent on the value of the first distribution and so on... This is implemented to overcome the limitations of the current Turing's implementation for dependent priors with dynamic support. See the following issues for more details: [1],[2],[3].

TwoLogUniformDependent(a, b, ϵ)
TwoLogUniformDependent(a, b) (constructor with default ϵ = 1e-10

Multivariate distribution to model three random variables where the first one x1 is given by log-U[a,b] and the second one x2 is given by log-U[x1,b].

• a: lower bound of the first distribution
• b: upper bound of the first distribution
• ϵ: small value to make sure that the lower and upper bounds of each distribution are different

This means that the lower bound of the second distribution is dependent on the value of the first distribution. This is implemented to overcome the limitations of the current Turing's implementation for dependent priors with dynamic support. See the following issues for more details: [1],[2],[3].

TwoUniformDependent(a, b, c, ϵ)
TwoUniformDependent(a, b, c) (constructor with default ϵ = 1e-10)

Multivariate distribution to model two random variables where the first one is given by U[a,b] and the second one is given by U[x,c], where x is a random variable sampled from the first distribution.

• a: lower bound of the first distribution
• b: upper bound of the first distribution
• c: upper bound of the second distribution
• ϵ: small value to make sure that the lower and upper bounds of each distribution are different

This means that the lower bound of the second distribution is dependent on the value of the first distribution.This is implemented to overcome the limitations of the current Turing's implementation for dependent priors with dynamic support. See the following issues for more details: [1],[2],[3].

Example

#```jldoctest julia> using Pioran, Distributions julia> d = TwoUniformDependent(0, 1, 2) TwoUniformDependent(0.0, 1.0, 2.0)

julia> rand(d) 2-element Array{Float64,1}: 0.123 1.234 #```

Bijectors.bijector(d::TwoUniformDependent)

Create a bijector for the TwoUniformDependent distribution. This is used to sample from the distribution using the Bijectors package. Adapted from the following issues [1],[2],[3].

CARMA(p, q, rα, β, norm, is_integrated_power=true)

Continuous-time AutoRegressive Moving Average (CARMA) model for the power spectral density

• p: the order of the autoregressive polynomial
• q: the order of the moving average polynomial
• rα: roots of the autoregressive polynomial length p+1
• β: the moving average coefficients length q+1
• norm: the normalisation of the process
• is_integrated_power: true if norm corresponds to the integral of the power spectrum between two frequencies (i.e. something like the sample variance). false if norm corresponds to the variance of the process, i.e. the integral of the PSD from 0 to infinity.

The power spectral density of the CARMA model is given by:

\[\mathcal{P}(f) = \sigma^2 \left|\dfrac{\sum\limits_{k=0}^q \beta_k \left(2\pi\mathrm{i}f\right)^k }{\sum\limits_{l=0}^p \alpha_l \left(2\pi\mathrm{i}f\right)^l}\right|^2\]

evaluate(model::CARMA, f)

evaluate the power spectral density of the CARMA model at frequency f.

quad2roots(quad)

Convert the coefficients of a quadratic polynomial to its roots.

Arguments

• quad::Vector{Real}: Coefficients of the quadratic polynomial.

Returns

• r::Vector{Complex}: Roots of the polynomial.

roots2coeffs(r)

Convert the roots of a polynomial to its coefficients.

Arguments

• r::Vector{Complex}: Roots of the polynomial.

Returns

• c::Vector{Complex}: Coefficients of the polynomial.