ReactiveBayes
-
|
I would like to use RxInfer for parameter estimation from stochastic differential equations. Here is my attempt to do that using a simple Ornstein-Uhlenbeck process. I am sharing my code, but it fails. I am wondering how one would formulate the problem in RxInfer. here is the stacktrace: |
Beta Was this translation helpful? Give feedback.
Replies: 2 comments 7 replies
-
|
Hey, thanks for trying out The original problem in your stacktrace is due to the specification of the Even with the correct using ExponentialFamilyProjection
using Plots
using Statistics
using StatsPlots
using RxInfer
using Distributions
using Random
using BayesBase
function ou(A,τ,N,Δt)
B = exp(-Δt/τ)
x = [rand(Normal(0,sqrt(A)))]
for i in 2:N
push!(x,rand(Normal(x[end]*B,sqrt(A*(1-B^2)))))
end
return x
end
# Fuse both deterministic functions into a brand-new distribution
struct TransformedGaussianDistribution{A, B, C, D}
x_prev::A
delta_t::B
τ::C
ampl::D
end
function mean_transform(x_prev, delta_t, t)
return x_prev * exp(-delta_t/t)
end
function std_transform(ampl, delta_t, τ)
return ampl*sqrt(1-exp(-delta_t/τ)^2)
end
function BayesBase.logpdf(d::TransformedGaussianDistribution, out)
transformed_dist = Normal(
mean_transform(d.x_prev, d.delta_t, d.τ),
std_transform(d.ampl, d.delta_t, d.τ)
)
return logpdf(transformed_dist, out)
end
BayesBase.insupport(d::TransformedGaussianDistribution, out) = true
@node TransformedGaussianDistribution Stochastic [ out, xprev, deltat, τ, ampl ]
@model function ou_noise(x, delta_t)
ampl ~ Gamma(1.0, 1.0)
τ ~ Gamma(1.0, 1.0)
x[1] ~ Normal(mean = 0.0, precision = 1.0)
for i in 2:length(x)
x[i] ~ TransformedGaussianDistribution(x[i - 1], delta_t, τ, ampl)
end
end
xx = ou(1.0,1.0,1000,0.5);
constraints = @constraints begin
q(ampl) :: ProjectedTo(Gamma)
q(τ) :: ProjectedTo(Gamma)
q(ampl, τ) = MeanField()
end
initialization = @initialization begin
q(ampl) = Gamma(1, 1)
q(τ) = Gamma(1, 1)
end
result = infer(
model = ou_noise(delta_t=0.5),
data = (x = xx,),
constraints = constraints,
initialization = initialization,
iterations = 20,
showprogress = true,
options = (
limit_stack_depth = 500,
rulefallback = NodeFunctionRuleFallback(),
)
)
using StatsPlots
@gif for (i, q) in enumerate(zip(result.posteriors[:ampl], result.posteriors[:τ]))
q_ampl = q[1]
q_τ = q[2]
p1 = plot(q_ampl, label = "q(ampl)", fill = 0, fillalpha = 0.2)
p1 = vline!([ mean(q_ampl) ], label = "mean = $(round(mean(q_ampl); digits = 2))", fill = 0, fillalpha = 0.2)
p2 = plot(q_τ, label = "q(τ)", fill = 0, fillalpha = 0.2)
p2 = vline!([ mean(q_τ) ], label = "mean = $(round(mean(q_τ); digits = 2))", fill = 0, fillalpha = 0.2)
plot(p1, p2; title = "Iteration $i")
endI made some changes to the model structure and replaced the This example uses Ideally, you will be able to try this yourself next week. I will ping you when we release this new functionality.
|
Beta Was this translation helpful? Give feedback.
-
|
With reduced number of samples and number of VMP iterations I managed to reduce the computation time from 1minute to about 1 second constraints = @constraints begin
parameters = ProjectionParameters(
strategy = ExponentialFamilyProjection.ControlVariateStrategy(
nsamples = 500
)
)
q(ampl) :: ProjectedTo(Gamma; parameters = parameters)
q(τ) :: ProjectedTo(Gamma; parameters = parameters)
q(ampl, τ) = MeanField()
end
|
Beta Was this translation helpful? Give feedback.
-
|
Thank you. That was very helpful. I could not find an example like this in your tutorials. |
Beta Was this translation helpful? Give feedback.
-
|
When I run your code, I get: I could not find this function in any of your packages. |
Beta Was this translation helpful? Give feedback.
-
|
Thank you for providing us with a use-case! It's always fascinating to see how people are using As a side note, I experimented with different input parameters for the xx = ou(0.5, 5.0, 1000, 0.5);
xx = ou(2.5, 0.5, 1000, 0.5);
In the first example, |
Beta Was this translation helpful? Give feedback.
-
|
@hstrey as I mentioned in my original message:
and
This functionality has not been released yet, but it will be very soon. We plan to release it early next week. |
Beta Was this translation helpful? Give feedback.
-
|
Thanks. We are looking into RxInfer as a potential tool to perform parameter inference on Stochastic Differential Equations. |
Beta Was this translation helpful? Give feedback.
-
|
@hstrey The new version 3.5.0 is now publicly available, please let us know if you have some issues with the new release. |
Beta Was this translation helpful? Give feedback.
-
|
Thanks. The code is working for me. I will explore further. I consider my question answered. |
Beta Was this translation helpful? Give feedback.
Hey, thanks for trying out
RxInfer!The original problem in your stacktrace is due to the specification of the
Normaldistribution. You should use eitherNormal(mean = ..., precision = ...)orNormal(mean = ..., variance = ...).Even with the correct
Normalspecification, this model is challenging for the typical use of RxInfer.I tried running it on an unreleased version of the package (which will be available soon), and here is my attempt: