Forecasting sales volume and sales revenue is an essential part of the operations of a business. In Part 1, we analyzed a paper to understand some principles supporting the end goal:
We want to make use of a bayesian reactive message passing approach, and have an eventual implementation in the RxInfer framework. In Part 2 we implemented a rough foundation by simulating 1,000 sales data points as a standardized normal distribution, making use of a 5th order autoregression process. We then used RxInfer to learn the order as well as the parameters \(\boldsymbol\theta\) of the process. The inference process also learned the value of the transition precision parameter \(\lambda_x\) as well as the \(x_t\) state component.
versioninfo() ## Julia versionJulia Version 1.12.7
Commit 6d172b025e4 (2026-08-15 08:05 UTC)
Build Info:
Official https://julialang.org release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 12 × Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
WORD_SIZE: 64
LLVM: libLLVM-18.1.7 (ORCJIT, skylake)
GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 12 virtual cores)import Pkg
Pkg.add(Pkg.PackageSpec(;name="RxInfer", version="3.6.0"))
Pkg.add(Pkg.PackageSpec(;name="Distributions", version="0.25.108"))
Pkg.add(Pkg.PackageSpec(;name="Random"))
Pkg.add(Pkg.PackageSpec(;name="Plots", version="1.41.6"))
Pkg.add(Pkg.PackageSpec(;name="LaTeXStrings", version="1.3.1"))
Pkg.add(Pkg.PackageSpec(;name="BenchmarkTools", version="1.5.0"))
Pkg.add(Pkg.PackageSpec(;name="Parameters", version="0.12.3"))
using RxInfer, Distributions, LinearAlgebra, Random, Plots, LaTeXStrings, BenchmarkTools, Parameters Updating registry at `~/.julia/registries/General.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`Pkg.status()Status `~/.julia/environments/v1.12/Project.toml`
⌃ [6e4b80f9] BenchmarkTools v1.5.0
⌃ [31c24e10] Distributions v0.25.108
⌃ [b964fa9f] LaTeXStrings v1.3.1
⌃ [d96e819e] Parameters v0.12.3
⌃ [91a5bcdd] Plots v1.41.6
⌃ [86711068] RxInfer v3.6.0
[9a3f8284] Random v1.11.0
Info Packages marked with ⌃ have new versions available and may be upgradable.## macro h(x) ## for help
# quote
# display("text/markdown", @doc $x)
# end
# end
# @h Normal() ## help on Normal()Our model for the Variational Bayesian Inference of a latent autoregressive model can be represented mathematically as follows:
\[\begin{aligned} p(\lambda_x) &= \mathcal{Gam}(\lambda_x;\ a, b)\\ p(\boldsymbol{\theta}) &= \mathcal{N}(\boldsymbol{\theta};\ \mathbf{m}_{\theta}, \mathbf{V}_{\theta})\\ p(\mathbf{x}_t \mid \mathbf{x}_{t-1}, \boldsymbol{\theta}) &= \mathcal{N}\left(\mathbf{x}_t;\ \mathbf{A}\left(\boldsymbol{\theta}\right) \mathbf{x}_{t-1},\, \mathbf{V}(\lambda_x)\right)\\ p(y_t \mid \mathbf{x}_t) &= \mathcal{N}(y_t;\ \mathbf{x}_t, \lambda_y^{-1}) \end{aligned}\]
where
\(M^{AR}\) is the order of the autoregression process
\(\boldsymbol{\theta}=(\theta_1, \theta_2, ..., \theta_{M^{AR}})^{\top}\) is a vector of transition coefficients
\(\mathbf{x}_{t-1} = (x_{t-1}, x_{t-2}, ..., x_{t-M^{AR}})^\top\) is the previous state
\(\mathbf{x}_t = (x_{t}, x_{t-1}, ..., x_{t-M^{AR}+1})^\top\) is the current state
\(y_t\) is the current observation
\(\lambda_x\) is the transition precision
\(\lambda_y\) is the observation precision
\(\mathbf{v}_u = (1,0,...,0)^\top\) is an \(M^{AR}\)-dimensional unit vector selecting the first component of the vector that follows
\(\mathbf{V}(\lambda_x) = (1/ \lambda_x)\mathbf{v}_u\mathbf{v}_u^\top\) is the transition covariance
\(\mathbf{A}(\boldsymbol{\theta}) = \begin{bmatrix} \boldsymbol{\qquad\theta}^\top \\ \mathbf{I}_{M^{AR}-1} & \mathbf{0} \end{bmatrix}\) is the transition matrix
## Unicode character for ꞉ to indicate time sequences
## Copy and paste this printed form to the rest of the code identifier
## The previously used \lmrk (triangular) colon silently failed
## in GraphPPL.
print(Char(0xA789)) ## prints ꞉
repr(UInt32(last("_xˣ꞉"))) ## verify "0x0000a789"꞉"0x0000a789"To generate some synthetic sales data, we use a \(M^{AR}=5\) order autoregressive (AR) process with the following coeffcients (corresponding to stable poles):
_θ⁵ = [0.10699399235785655, -0.5237303489793305, 0.3068897071844715, -0.17232255282458891, 0.13323964347539288];function generate_sales_data(
rng, ##. random number generator
N, ##. number of data points
_θˣ, ##. AR parameters
λˣₓ, ##. transition precision
λˣᵧ) ##. observation precision
Mᴬᴿ = length(_θˣ) ##. order of process (use Mᴹᴬ for MA process)
_xˣ꞉ = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
y꞉ = Vector{Float64}(undef, N + 3Mᴬᴿ)
λˣₓ_std = sqrt(inv(λˣₓ)) ##. transition standard dev
λˣᵧ_std = sqrt(inv(λˣᵧ)) ##. observation standard dev
_xˣ꞉[1] = randn(rng, Mᴬᴿ) ##. Mᴬᴿ random numbers
## +++ Matrix implementation
## Build _AIθₜI, the transition matrix
_Iₘ₋₁ = Diagonal(ones(Mᴬᴿ - 1))
_V0 = Vector(zeros(Mᴬᴿ - 1)) ## vector of zeros
_Aᵇᵒᵗᵗᵒᵐ = hcat(_Iₘ₋₁, _V0)
_θₜ = _θˣ
_AIθₜI = vcat(transpose(_θₜ), _Aᵇᵒᵗᵗᵒᵐ)
## Build _VIλₓI, the transition covariance matrix
_vᵤ = Vector(vcat(1.0, zeros(Mᴬᴿ - 1)))
_VIλₓI = (1/λˣₓ)*_vᵤ*transpose(_vᵤ)
for i in 2:Mᴬᴿ
_VIλₓI[i, i] = 1e-12*tiny ##tiny values on remaining off-diagonals to avoid PosDefException
end
## ---
for t in 2:(N + 3Mᴬᴿ) #.
## _xˣ꞉[t] = vcat(
## rand( ##. noised current value is pushed in at the top of state vector
## rng,
## Normal(dot(θˣ, _xˣ꞉[t-1]), λˣₓ_std)
## ),
## _xˣ꞉[t-1][1:end-1] ##. rest shifted down, last element pushed out & lost
## ) #.
_xˣ꞉[t] = rand(rng, MvNormal(_AIθₜI*_xˣ꞉[t-1], _VIλₓI)) ##. Eq (2b)
## y꞉[t] = rand(rng, Normal(xˣː[t][1], λˣᵧ_std)) #.
y꞉[t] = rand(rng, Normal(transpose(_vᵤ)*_xˣ꞉[t], λˣᵧ_std)) ##. Eq (2c)
end
return _xˣ꞉[1+3Mᴬᴿ:end], y꞉[1+3Mᴬᴿ:end] #.
endgenerate_sales_data (generic function with 1 method)## Seed for reproducibility
seed = 123
rng = MersenneTwister(seed)
## Number of observations in synthetic dataset
N = 1000
## AR process parameters
λˣₓ = 1.0 ##. real λₓ (transition precision)
λˣᵧ = 0.5 ##. real λᵧ (observation precision)
_θˣ = _θ⁵ ##. real θ (AR parameters)
_xˣ꞉, y꞉ = generate_sales_data( #.
rng, N, _θˣ, λˣₓ, λˣᵧ); #._xˣ꞉[1:7]7-element Vector{Vector{Float64}}:
[-1.5244304669899376, 0.7117529427153132, -0.6818253883618433, -1.347929820049962, 2.334795350741872]
[-0.4065895890454423, -1.524430466990528, 0.7117529427157929, -0.6818253883624347, -1.3479298200510053]
[0.941894974914772, -0.4065895890438795, -1.5244304669902413, 0.7117529427167336, -0.6818253883624831]
[-0.94636366590129, 0.9418949749154402, -0.4065895890428543, -1.5244304669908302, 0.711752942718009]
[0.24712361685519496, -0.9463636659021836, 0.9418949749150592, -0.4065895890408091, -1.5244304669903614]
[2.564623886118386, 0.24712361685545106, -0.9463636659016156, 0.9418949749164531, -0.4065895890410718]
[0.8139262873419921, 2.564623886119036, 0.24712361685480724, -0.9463636659019005, 0.9418949749154859]y꞉[1:7]7-element Vector{Float64}:
0.07169900096517967
-0.6385825951402568
0.5606040423489563
1.6638830079733458
-0.30057340406379157
2.1719951691142283
-1.101725570531378Let’s plot our synthetic dataset:
plot(
first.(_xˣ꞉),
title="Simulated standardized sales data",
xlabel=L"$\mathrm{time}, t$",
color="red",
label=L"$\mathrm{Hidden\ state\ component}, \mathbf{x}^*_t$",
) #.
scatter!(y꞉, label=L"$\mathrm{Observation}, y_t$", color="orange")
Now we specify a probabilistic model, the inference constraints and then run an inference procedure with RxInfer. We need two different models for:
Multivariate AR, with order \(M^{AR} > 1\)
Univariate AR (reduces to a simple State-Space-Model), with order \(M^{AR}\) = 1
@model function lar_model(
y꞉, ##. data/observations
𝚃ᴬᴿ, ##. Uni/Multi variate
Mᴬᴿ, ##. AR order
_vᵤ, ##. unit vector
λᵧ) ##. observation precision
## Priors
λₓ ~ Gamma(α = 1.0, β = 1.0) ##. for transition precision
if 𝚃ᴬᴿ === Multivariate
_θ ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
_x₀ ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
else ## Univariate
_θ ~ Normal(μ = 0.0, γ = 1.0)
_x₀ ~ Normal(μ = 0.0, γ = 1.0)
end
_xₜ₋₁ = _x₀
for t in eachindex(y꞉)
_x꞉[t] ~ AR(_xₜ₋₁, _θ, λₓ) ##.Eq (2b)
if 𝚃ᴬᴿ === Multivariate
y꞉[t] ~ Normal(μ = dot(_vᵤ, _x꞉[t]), γ = λᵧ) ##.Eq (2c)
else
y꞉[t] ~ Normal(μ = _vᵤ*_x꞉[t], γ = λᵧ) ##.Eq (2c)
end
_xₜ₋₁ = _x꞉[t]
end
end@constraints function ar_constraints()
q(_x₀, _x꞉, _θ, λₓ) = q(_x₀, _x꞉)q(_θ)q(λₓ)
endar_constraints (generic function with 1 method)@meta function ar_meta(𝚃ᴬᴿ, Mᴬᴿ, stype)
AR() -> ARMeta(𝚃ᴬᴿ, Mᴬᴿ, stype) #.
endar_meta (generic function with 1 method)@initialization function ar_init(𝚃ᴬᴿ, Mᴬᴿ)
q(λₓ) = GammaShapeRate(1.0, 1.0)
if 𝚃ᴬᴿ === Multivariate
q(_θ) = MvNormalMeanPrecision(zeros(Mᴬᴿ), diageye(Mᴬᴿ)) #.
else ## Univariate
q(_θ) = NormalMeanPrecision(0.0, 1.0) #.
end
endar_init (generic function with 1 method)function run_inference(Nmodels)
results = Vector{Any}(undef, Nmodels)
for m = 1:Nmodels
println("m = $m")
if m > 1 𝚃ᴬᴿ = Multivariate elseif m === 1 𝚃ᴬᴿ = Univariate
else error("m=$m which is not an allowed value.") end
results[m] = infer(
model = lar_model(
𝚃ᴬᴿ=𝚃ᴬᴿ, ##. Uni/Multi variate
Mᴬᴿ=m, ##. AR order
_vᵤ=ReactiveMP.ar_unit(𝚃ᴬᴿ, m), ##. unit vector
λᵧ=λˣᵧ), ##. observation precision
data = (y꞉=y꞉,), #.
constraints = ar_constraints(),
meta = ar_meta(𝚃ᴬᴿ, m, ARsafe()), #.
options = (limit_stack_depth = 100,),
initialization = ar_init(𝚃ᴬᴿ, m),
returnvars = (_x꞉=KeepLast(), λₓ=KeepEach(), _θ=KeepEach()),
free_energy = true,
iterations = 100,
showprogress = false
)
end
return(results)
end run_inference (generic function with 1 method)Now we run the inference for \(M^{AR}=1,2, ..., 5\):
results = run_inference(5)m = 1
m = 2
m = 3
m = 4
m = 55-element Vector{Any}:
Inference results:
Posteriors | available for (λₓ, _x꞉, _θ)
Free Energy: | Real[2000.24, 1998.79, 1997.69, 1996.65, 1995.67, 1994.76, 1993.92, 1993.15, 1992.46, 1991.84 … 1988.34, 1988.34, 1988.34, 1988.34, 1988.34, 1988.34, 1988.34, 1988.34, 1988.34, 1988.34]
Inference results:
Posteriors | available for (λₓ, _x꞉, _θ)
Free Energy: | Real[2014.05, 2009.48, 2005.92, 2001.9, 1997.7, 1993.1, 1988.4, 1983.81, 1979.68, 1976.42 … 1967.94, 1967.95, 1967.96, 1967.98, 1967.95, 1967.96, 1967.97, 1967.97, 1967.98, 1967.99]
Inference results:
Posteriors | available for (λₓ, _x꞉, _θ)
Free Energy: | Real[2026.8, 2020.77, 2017.93, 2014.62, 2010.94, 2006.75, 2002.33, 1997.5, 1992.68, 1988.14 … 1970.96, 1970.91, 1970.96, 1970.93, 1971.01, 1971.0, 1970.95, 1971.02, 1970.96, 1970.98]
Inference results:
Posteriors | available for (λₓ, _x꞉, _θ)
Free Energy: | Real[2038.5, 2030.06, 2027.7, 2025.15, 2022.27, 2018.96, 2015.2, 2011.2, 2006.9, 2002.27 … 1968.65, 1968.66, 1968.53, 1968.64, 1968.44, 1968.69, 1968.79, 1968.67, 1968.65, 1968.71]
Inference results:
Posteriors | available for (λₓ, _x꞉, _θ)
Free Energy: | Real[2049.95, 2037.96, 2036.03, 2033.91, 2031.5, 2028.88, 2025.89, 2022.59, 2018.9, 2014.87 … 1969.73, 1969.86, 1969.77, 1969.91, 1969.61, 1969.79, 1969.77, 1969.69, 1969.56, 1969.66]function plot_free_energy(results)
p = plot(title="Bethe Free Energy", legend=:topright, xlabel="iterations")
for m = 1:length(results)
p = plot!(p, results[m].free_energy, label="AR$m")
end
return p
endplot_free_energy (generic function with 1 method)plot_free_energy(results)
function plot_converged_free_energy(results)
values = [results[m].free_energy[end] for m in 1:length(results)]
bar(
1:length(results),
values,
legend=:none,
ylim=(1960, 2000),
title="Converged Free Energy vs AR Order"
)
xlabel!("AR order")
ylabel!("Minimum Free Energy")
endplot_converged_free_energy (generic function with 1 method)plot_converged_free_energy(results)
values = [results[m].free_energy[end] for m in 1:length(results)]5-element Vector{Float64}:
1988.3356732027921
1967.988103133217
1970.9795806202876
1968.7069285195203
1969.6558516038376function plot_transition_precision(results)
p = plot(title="Inferred transition precision", legend=:topright, xlabel="iterations")
for m = 1:length(results)
p = plot!(p, mean.(results[m].posteriors[:λₓ]), label="AR$m")
end
plot!(
p, [λˣₓ], seriestype=:hline, color="red", linestyle=:dash,
label=L"Real transition precision $\lambda_x^*$"
)
return p
endplot_transition_precision (generic function with 1 method)plot_transition_precision(results)
best_order = 2
@unpack _x꞉, _θ, λₓ = results[best_order].posteriorsDict{Symbol, Vector} with 3 entries:
:λₓ => GammaShapeRate{Float64}[GammaShapeRate{Float64}(a=501.0, b=176.311), …
:_x꞉ => MvNormalWeightedMeanPrecision{Float64, Vector{Float64}, Matrix{Float6…
:_θ => MvNormalWeightedMeanPrecision{Float64, Vector{Float64}, Matrix{Float6…p1 = plot(
first.(_xˣ꞉), ## first.() gives xˣₜ, not the delayed xˣₜ₋₁, xˣₜ₋₂, ..., xˣₜ₋ₘ₊₁
title="AR, " * L"$M^{AR}=$" * "$best_order",
xlabel=L"time $t$",
color="red",
label=L"$\mathrm{Hidden\ state\ component}, \mathbf{x}^*_t$")
p1 = scatter!(p1, y꞉, label=L"$\mathrm{Observation}, y_t$", color="orange")
p1 = plot!(
p1, first.(mean.(_x꞉)),
ribbon=first.(std.(_x꞉)),
fillalpha=0.2,
label=L"$\mathrm{Inferred\ state\ component}, \mathbf{x}_t$",
color="green",
legend=:bottomright)
p2 = plot(
mean.(λₓ), ribbon=std.(λₓ),
fillalpha=0.2,
label=L"Inferred transition precision $\lambda_x$",
color="green",
legend=:topright,
xlabel="iterations",
)
p2 = plot!(
[λˣₓ],
seriestype=:hline,
color="red", #.
label=L"Real transition precision $\lambda_x^*$"
)
p3 = plot(results[best_order].free_energy, label="Bethe Free Energy", xlabel="iterations")
plot(p1, p2, p3, layout=@layout([ a; b c ]))
Let’s look at a smaller range of data for the sake of interest:
subrange = 700 : 750
plot(
subrange,
first.(_xˣ꞉)[subrange],
title=L"Subrange of: AR, $M^{AR}=$"*"$best_order",
xlabel=L"time $t$",
color="red",
label=L"$\mathrm{Hidden\ state\ component}, \mathbf{x}^*_t$") #.
scatter!(subrange, y꞉[subrange], label=L"$\mathrm{Observation}, y_t$", color="orange")
plot!(
subrange,
first.(mean.(_x꞉))[subrange],
ribbon=sqrt.(first.(var.(_x꞉)))[subrange],
fillalpha=0.2,
label=L"$\mathrm{Inferred\ state\ component}, \mathbf{x}_t$",
color="green",
legend=:bottomright)
Finally, we look at how the learned AR coefficients converge:
let
pθ = plot()
θms = mean.(_θ) ##. _θ means
θvs = var.(_θ) ##. _θ variances
l = length(θms)
edim(e) = (a) -> map(r -> r[e], a)
for t in 1:length(first(θms)) #.
tmp = latexstring("\$ \\theta_{$(t),t} \$")
pθ = plot!(
pθ,
θms |> edim(t),
ribbon=θvs |> edim(t) .|> sqrt,
fillalpha=0.2,
label="Inferred $tmp")
end
for t in 1:length(_θˣ)
tmp = latexstring("\$ \\theta_{$(t),t} \$")
pθ = plot!(
pθ,
[ _θˣ[t] ],
seriestype=:hline,
label="Real $tmp")
end
plot(
pθ,
title=L"AR coefficients, $M^{AR}=$"*"$best_order",
xlabel=L"time $t$",
legend=:outertopright,
size=(800, 450)) #.
end