API

BayesianDiscovery

Here is my package.

Posterior Struct

Initiates a structure of class Posterior to hold posterior samples

Arguments

DEtection(Y, SpaceStep, TimeStep, νS, νT, bufferSpace, bufferTime, batchSpace, batchTime, learning_rate, v0, v1, Λ, Λnames)

DEtection sampler function. Can accept missing values in the input data argument. Returns the model, parameters, and posterior values. The model are the model settings. The parameters are the final value of the parameters in from the sampler. The posterior are the saved posterior values.

Consider the function $U_t = M(U, U_x, U_y, U_{xy}, ...)$. DEtection() is used to determine $M$ given a library of potential values, Λ(U, ∂U), to search over. Within the function, partial derivatives are denoted as $ΔU_t, ΔU_x, ΔU_y, ΔU_xx, , ΔU_xy, ...$, so a potential function could be

function Λ(U, ∂U){
  u = U[1]
  u_x = ∂U[1]
  u_y = ∂U[2]

  return [u, u_x, u_y, u*u_x, u*u_y]
}

To make the function identify the correct partial derivatives, the argument that are passed into ∂U, Λnames, are required. For the example above, Λnames = [ΔU_x, ΔU_y] because the function Λ uses the partial derivatives $U_x$ and $U_t$.

If you want to add covariates to the function, a possible Λ(U, ∂U, X) is

function Λ(U, ∂U, X){
  u = U[1]
  u_x = ∂U[1]
  u_y = ∂U[2]
  x1 = X[1]
  x2 = X[2]

  return [u, u_x, u_y, u*u_x, u*u_y, x1, x2, u*x1, u_x*x2]
}

where Λnames = [ΔU_x, ΔU_y].

Required Arguments (in order)

• Y: Input data. Needs to be Array{Float64, 3} or Array{Union{Missing, Float64}, 3} where the dimensions are Space, Time, Components
• SpaceStep: of type StepRangeLen (range function). For example, range(-1, 1, step = 0.1) for 1 dimension and [range(-1, 1, step = 0.1), range(-1, 1, step = 0.1)] for 2.
• TimeStep: of type StepRangeLen (range function). For example, range(-1, 1, step = 0.1).
• νS::Int or νS::Vector{Int}: Number of spatial basis functions.
• νT::Int: Number of temporal basis functions
• batchSpace::Int:
• batchTime::Int:
• learning_rate::Float64:
• v0::Float64:
• v1::Float64:
• Λ::Function:
• Λnames::Vector{String}:

Optional Arguments

• response = "ΔUt": Order of the temporal derivative (default first order). Use "ΔUtt" for second and so on.
• degree = 4: Degree of the B-spline. Must be at least one order higher than the highest order partial derivative.
• orderTime = 1: Order of the highest order temporal derivative (default ∂U_t)
• orderSpace = 3: Order of the highest order spatial derivative (default ∂Uxxx and ∂Uyyy)
• latent_dim = size(Y, 3): Dimension of the latent space. Default is same as data dimension.
• covariates = nothing: Additional covariates.
• nits = 2000: Number of samples for the Gibbs sampler.
• burnin = nits / 2: Number of samples to discard as burnin (default is half of nits).
• learning_rate_end = learning_rate: End learning rate (default is same as initial learning rate).

bspline(x, knot_locs, degree, derivative)

Creates a matrix of B-Splines evaluated at each x with specified knot locations and degree. Dimension is (length(x), length(xknotlocs)). Returns the basis matrix and the derivative of the basis matrix.

Arguments

• x: data
• knot_locs: knot locations
• degree: degree of the B-Spline
• derivative: order of the derivative for Psi_x

Examples

x = 1:1:30
knot_locs = 5:5:25
degree = 3
derivative = 1

Psi, Psi_x = bspline(x, knot_locs, degree, derivative)

create_pars()

Constructs the parameter and model classes.

Arguments

• Y::Array{Float64, 3}: Space x Time x Component data array
• SpaceStep::Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}: Vector of spatial locations [x, y]
• TimeStep::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}: Time samples
• νS::Vector{Int32}: number of spatial basis functions [x, y]
• νT::Int: number of temporal basis functions
• bufferSpace::Vector{Int}: spatial buffer [x, y]
• bufferTime::Int: temporal buffer
• batchSpace::Int: batch size for space [x, y]
• batchTime::Int: batch size for time
• learning_rate::Float64: learning rate
• v0::Float64: ssvs not included variance
• v1::Float64: ssvs included variance
• Λ::Function: function library
• Λnames::Vector{String}: names of arguments in function library

fold3(Y)

Used to construct a tensor from the mode-3 matrix Y.

Examples

Y3 = unfold3(Y) I = size(Y, 1) J = size(Y, 2) K = size(Y, 3)

Y = fold3(Y3, I, J, K)

spatial_bspline(x, y, x_knot_locs, y_knot_locs, degree)

Creates a matrix of B-Splines evaluated at each x and y with specified knot locations and degree. Dimension is (length(x)length(y), length(xknotlocs)length(yknotlocs)). Returns the basis matrix and the derivative of the basis matrix.

Arguments

• x: data in the x direction
• y: data in the x direction
• xknotlocs: knot locations for x
• yknotlocs: knot locations for y
• degree: degree of the B-Spline

Examples

x = 1:1:30
y = 1:1:30
x_knot_locs = 5:5:25
y_knot_locs = 5:5:25
degree = 3

Psi, Psi_x, Psi_y, Psi_xy = spatial_bspline(x, y, x_knot_locs, y_knot_locs, degree)

# plot the functions
Plots.pyplot()
contour(x, y, reshape(Psi, length(x), length(y)), fill = true)
contour(x, y, reshape(Psi_x, length(x), length(y)), fill = true)
contour(x, y, reshape(Psi_y, length(x), length(y)), fill = true)
contour(x, y, reshape(Psi_xy, length(x), length(y)), fill = true)

tensor_mult(G, A, B, C)

Used to multiply a tensor G ∈ R(I1, I2, I3) with matrices A ∈ R(P, I1), B ∈ R(Q, I2), C ∈ R(R, I3)

unfold3(Y)

Used to get the mode-3 matrix of the tensor Y.

update_A!(pars)

Used within DEtection() function with a spike-and-slab prior. Updates 𝐀 with the elastic net prior (Li 2010) from

$𝐳(𝐬, t) = ℋ(𝐬,t) 𝚯 𝐀 (𝛗phi_{t^{(0)}}(t) ⊗ 𝛙(𝐬))' + 𝛜(𝐬, t)$

$𝚯 𝐀 (𝛗_{t^{(i)}}(t) ⊗ g(𝛙(𝐬)))' = 𝐌 𝐟(⋅) + 𝛈(𝐬, t)$

$p(𝐀) ∝ exp{-λ₁‖𝐀‖₁ - λ₂‖𝐀‖₂²}$

using Stochastic Gradient Desent with a Constant Learning Rate (Mandt 2016). λ₁ and λ₂ are set to 1/100.

update_M!(pars)

Used within DEtection() function with a spike-and-slab prior. Updates 𝐌 and 𝚺ᵤ from

$𝚯 𝐀 (𝛗_{t^{(i)}}(t) ⊗ g(𝛙(𝐬)))' = 𝐌 𝐟(⋅) + 𝛈(𝐬, t)$

where $𝛈(𝐬, t) ∼ N_N(0, 𝚺ᵤ)$and $𝚺ᵤ = diag(σ²ᵤ1, ..., σ²ᵤN)$

update_gamma!(pars)

Used within DEtection() function with a spike-and-slab prior. Updates gamma and pi from the spike-and-slab prior from

$p(𝐌⎹ 𝛄) = p-slab(𝐌 ᵧ)∏ⱼ p-spike(Mⱼ)$

where $p(γⱼ = 1⎹ π) = π$, $π ∼ ℬ(a, b)$, $p-slab(𝐌 ᵧ) = N(𝐌 ᵧ⎹ 𝐚ᵧ, 𝚺ᵤ𝐀ᵧ)$, and $𝐀ᵧ = c𝐈$.

update_ΣZ!(pars)

Used within DEtection() function with a spike-and-slab prior. Updates 𝚺z from

$𝐳(𝐬, t) = ℋ(𝐬,t) 𝚯 𝐀 (𝛗phi_{t^{(0)}}(t) ⊗ 𝛙(𝐬))' + 𝛜(𝐬, t)$

where $𝛜(𝐬, t) ∼ N_N(0, 𝚺z)$, $𝚺z = diag(σ²z1, ..., σ²zm)$, and $σ²z ∼ Half-t()$.

ΔL(z, H, ψ, gψ, ϕ, ϕ_t, Θ, ΣZinv, ΣUinv, A, M, fcurr, fprime)

Used within DEtection() function. Calculates the gradient of the log likelihood.