continious Active Inference on the RxInfer platform

[Charelvanhoof]

Hi Team,

Thank you for our excellent meeting/discussion yesterday. Looking forward to see how we can make (continuous) active inference work on the RxInfer platform. It looks promising. I guess we might need to add some features to the platform to capture the specifics of continuous active inference.

As requested I opened a topic over here to get things going.

As discussed, as a first step, lets recreate the IWAI 2024 tutorial (Basic active inference capsule) on RxInfer. Here is the link: https://www.kaggle.com/code/charel/active-inference-example-iwai2024 (including working python code, so we can double check we get the same results).

The example is set in a simulation environment where an early evolutionary imaginary aquatic ancestor (Hydra vulgaris) must preserve its physical integrity to survive. For example, it needs to keep in a certain temperature range. It is a basic case of a single active inference capsule with one sensor, one hidden state and one causal state. And Hydar can take action to move up or down.

I propose to build the example step by step:

Step 1.1 Inference without generalized coordinates
Inference of the hidden state, causal state (without generalized coordinates). In the example the function of sensory mapping g is non-linear. If this is an extra challenge we can take a first step with a linear g function
Step 1.2 Add the action
Realtime/online inference including the inference of the reflex action u, which is input for the next timestep (it is not discrete/EFE active inference)
Step 1.3 Add the generalized coordinates
Need to express the states in generalized coordinates e.g [position, speed, acceleration,..] and inference all generalized coordinates. (see the generalized coordinates as a kind of short-term-prediction). Need to be able to generate colored noise and the precision matrix of the colored noise.

Thank for the help because I am new to RxInfer/Julia.

[bvdmitri]

@ThijsvdLaar @ofSingularMind do you want to collaborate on this?

[ThijsvdLaar]

@Charelvanhoof Thanks for opening this discussion and for sharing your notebooks, they are a very nice resource. I like the idea of an education-oriented narrative for explaining continuous-time active inference (CT AIF), with open-source notebooks that exemplify crucial design decisions in the framework.

@Charelvanhoof @ofSingularMind If you're interested, we can start to take CT AIF apart so to speak, and then explain the ingredients that make it tick. I think that alternative design choices for CT AIF can be made (along the way), and RxInfer can play a role here (even without much modification).

For example, I've been trying to implement Friston and Ao's 2012 paper on "Free Energy, Value and Attractors" with RxInfer, and see the following connections:

1. Desired dynamics are described by the generative model, which can be represented by a factor graph;
2. Veridical dynamics (depending on action) are not explicitly represented by a model, but by an equation (Friston & Ao, Eq. H.3) that directly computes the action gradient of the generative process. I see this equation as a result of inference on a (deterministic) "veridical" model that represents the function y(a). Then, this model can also be represented by a factor graph;
3. Generalized coordinates of motion can be included in these models. However, I think this is a design choice (namely, the one that leads to a biologically plausible inference scheme). In contrast, inference from an engineering viewpoint may lead to an (extended) Kalman filter, which (with or without generalized coordinates) can be implemented with RxInfer.

I've been meaning to collect examples of reactive AIF implementations into a public repository, https://github.com/biaslab/RxAIF, so others may learn and take inspiration (but I haven't gotten very far yet). I can already share some materials on CT AIF in this repo (throughout the coming weeks, because I'm currently on leave), and we can collaborate if you like.

[bertdv]

Hi all-- Thanks @Charelvanhoof for opening this discussion. Perhaps I can add my view on this. I think there are two different threads in this issue. One is the simulation of continuous-time dynamical systems, and the other is the simulation of active inference (AIF) processes. I recommend that we consider these issues separately.

In principle, all dynamical systems are simulated using discrete time steps, regardless of whether they are described on paper by differential equations (i.e., continuous-time systems) or by difference equations (i.e., discrete-time systems). The default in RxInfer is to define and simulate systems by difference equations. @Ismail Senoz has worked on the simulation of CT dynamical systems in RxInfer, which essentially comes down to simulating by adaptive time step durations. We can ask @ismailsenoz for more input here.

The simulation of an AIF system (agent + environment) is nothing but the simulation of a dynamical system. By itself, transitioning from discrete-time to continuous-time simulation (or vice versa) is conceptually unrelated to AIF—it's merely a change in the simulation paradigm. So, going forward, I recommend using difference equations for both model specification and for simulation of the inference process, as RxInfer far better supports this. Once that works, we can consider whether the dynamical AIF system can also be specified by differential equations and if the continuous-time inference process can be simulated in RxInfer using the techniques that Ismail has worked on.

As for the second thread, how to simulate (discrete-time) AIF processes in RxInfer, I would consider the following. The expected free energy (EFE) is not a variational free energy (VFE), but rather a cost function for policies. In almost all of the AIF literature, it is suggested to compute EFE for a "smartly selected" set of policy candidates (represented by: u) and then assign a prior "softmax(-G(u))" to those candidates. This is a method that (1) is not consistent with the FEP idea that FE minimization is the only driving process, and (2) does not scale well (since the "smartly selected"-issue leads to tree-search algorithms that explode for large-dimensional state spaces and/or deep horizons). I recommend following Theorem 1 in https://arxiv.org/abs/2504.14898, which reframes EFE-based planning as an emergent property of regular variational free energy (VFE) minimization in a generative model-rollout into the future, terminated by both target priors and epistemic priors. This approach is, in principle, directly realizable in RxInfer.

As an example, consider the slide below, which shows some key code fragments of multiple agents that need to move to known end points without collisions.

The code already includes a generative model and target priors. Inference in this model leads to the agents being moved by a "KL control" driving force. In order to turn this agent into an AIF agent, the code only needs to be extended by adding epistemic priors (eqs 9 in https://arxiv.org/abs/2504.14898). This is not entirely straightforward, but we are working on a nice user interface for RxInfer to support this. @Wouter Nuijten is the most informed person about this work. Going forward, I recommend syncing with @wouterwln on when and how epistemic priors are supported in RxInfer.

In summary, I recommend proceeding as follows:

1. Transfer the example of the Hydra vulgaris as an AIF agent to a discrete-time (DT) model specification (by difference equations). The outcome of this exercise should be a generative model (that specifies how the agent believes its observations will be generated as a function of its actions), and target priors (that specify the preferences over future states as a probability distribution).
2. Implement the DT version of the AIF agent in Python, and check that the behavior is similar to the CT version. (for this check, it is necessary to develop a few benchmark tests).
3. Code the generative model, including the epistemic and target priors, in RxInfer (following the example of the multi-agents above). This should be less than 1 page of code.
4. Run the inference process in RxInfer, and check that the Hydra vulgaris AIF agent in RxInfer behaves similarly (on the benchmark tests) as the AIF agent in Python.

Generally speaking, BIASlab and Lazy Dynamics members can help with tasks 3 and 4 if you encounter problems. In particular, @wouterwln is the expert in helping with task 3. @bvdmitri and other Lazy Dynamics members are experts if task 4 (the inference process) fails. Everybody else in BIASlab can also help, either by themselves or in finding appropriate help.

If the above works, then the same recipe can be used for all other agents. As discussed, the automatic simulation of continuous-time inference processes in RxInfer is a separate issue that needs the involvement of @ismailsenoz. I recommend not considering it until we have many AIF agents working smoothly by discrete-time simulations.

[Charelvanhoof]

First of all, my apologies, I missed these replies. All notifications went to my old Volvocars account, where I worked in the past. And I did not have in my routine to check the github discussion pages. I have repaired the settings, so I can react next time.
I am not sure if I understand the reply about simulation of continuous-time dynamical systems, vs simulation of active inference (AIF) processes. Maybe a reference because in the notebooks there is a generative model and a generative process?
In my mind RxInfer example notebook should have the focus on the active inference part with the simple generative model, left side of the slide. It has discrete time-steps, I selected dt = 0.005. The generative process/simulated world is also a simple model simulated using the same discrete time steps. I presume RxInfer can do something similar or if not, we might find a solution that the world simulation is a python program. Note that it is not multi-agent.
I do not want to change the Hydra example, because it is the canonical active inference basic example in the continuous domain. And I do suspect that it will simply fit in RxInfer (at least without generalized coordinates). However, i must first look at myself because I did not yet spend sufficient time on it. I will do so coming months. A sparring partner would still be helpful to get into RxInfer (and at the same time create a notebook so others can get into it easier)

[wmkouw]

Summary of our pair programming session:

• We debugged the current notebook that ran batch inference on a linearized state transition, nonlinear likelihood and 1-step ahead planning. This involved defining factorization constraints, initializing the variational distributions and limiting stack depth.
• We cast the model to a recursive form, where we can step over time in a discrete-time simulation, feed observations one-by-one, and call infer to update the state estimate.
• We added the predictvars keyword argument so RxInfer produces a predictive distribution for the observation y.

It should now be possible to calculate action with the prediction error formula. That action can then be fed to the environment at each timestep, which would generate the next observation (thus completing the loop).

Next steps are:

• Longer planning horizon. I have shared a notebook with RxInfer code for the Bayesian Thermostat in Thijs' paper. At BIASlab, we will update this code to the current version of RxInfer and send it to you.
• Setting up a reactive agent-environment loop, which replaces the for loop that steps over time with an event-based interaction.
• Increasing the depth of the generalized coordinate vector.

If there's anything I missed or needs to be corrected, let me know.

[Charelvanhoof]

My thanks and good meeting! I have been testing/reviewing/playing with the code and see if I can reproduce the typical active inference graphs. I have build the same typical example in RxInfer as I presented in IWAI

If I run the code I get this result:

while i was expecting this (type of) result:

The hidden state estimate will move between the prior and sensory observations, closer to the prior if the prior has a higher precision, and closer to the sensory observations if these have a higher precision. (and is this example both precisions given as input). With RxInfer code it looks like the prior is 'ignored' and the sensory observation drives the hidden state estimate. Note that in both examples there is no action (yet), it is pure hidden state inference.

I suspect the issue it is in the model definition, probably because x_init is a Gaussian and mu_v a real?
`
@model function Hydar_model(y, x_init, mu_v, dt, x_s2, y_s2)

x_prev ~ x_init
f = (mx::Real, mv::Real) -> -mx + mv
f_step = (mx::Real, mv::Real) -> mx + dt*f(mx, mv)
g = (x::Real) -> 25 - 16 / (1 + exp(5 - x / 5))


for i in 1:length(y)
    x[i] ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=x_s2) # x_t+1=x_t + dt*f(x_t, v_t) = xt+1 = xt -dt*x_t + dt*v_t
    #x[i] ~ Normal(mean= x_prev - mu_v , var=x_s2) 
    #xk[i] ~ f_step(x_prev,mu_v) where { meta = DeltaMeta(method = Linearization()) }
    #x[i]  ~ Normal(mean = xk[i], var=x_s2) 
    gk[i] ~ g(x[i]) where { meta = DeltaMeta(method = Linearization()) }
    y[i]  ~ Normal(mean = gk[i], var=y_s2)

    x_prev = x[i]
end

end
`
eg if I use the second line in the for loop (x[i] ~ Normal(mean= x_prev - mu_v , var=x_s2) ) I get an error but somehow the first line in the for loop doesn't give an error. the error is : Error: We encountered an error during inference, here are some helpful resources to get you back on track...

[Charelvanhoof]

Thank you for the discussion today, one more step forward. By skewing the precisions/variance of the prior the graph starts to be appearing. Increased the precision of the prior with a factor 100.

gp_x_init = 25.0   # actual depth of hydar
gp_x_p = 1.0       # actual precision of function of motion
gp_y_p = 1.0       # actual precision of function of sensory observations
gm_x_init = 15.0   # Hyadr's belief of depth
gm_x_target = 15.0 # Prior belief of Hydar
gm_x_p = 200.0       # Hyadrs'belief of precision of function of motion
gm_y_p = 1.0       # Hydars belief of precision of function of sensory observations

Resulting in:

However engineering the precisions with factor 100... And also at the end of the trail there is a "funny" decline.

After our discussion I had a brainwave that it was maybe in the function g because the input parameter was specified as a ::real iso of a distribution. But that didn't make the difference.

Model code now looks:

@model function Hydar_model(y, x_init, mu_v, dt, x_s2, y_s2)

    x_prev ~ x_init
    #f = (mx::Real, mv::Real) -> -mx + mv
    #f_step = (mx::Real, mv::Real) -> mx + dt*f(mx, mv)
    g = (x) -> 25 - 16 / (1 + exp(5 - x / 5))


    for i in 1:length(y)
        x[i] ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=x_s2) # x_t+1=x_t + dt*f(x_t, v_t) = xt+1 = xt -dt*x_t + dt*v_t
        # tmp ~ mu_v - x_prev
        # x[i] ~ Normal(mean= tmp , var=x_s2) 
        #xk[i] ~ f_step(x_prev,mu_v) where { meta = DeltaMeta(method = Linearization()) }
        #x[i]  ~ Normal(mean = xk[i], var=x_s2) 
        gk[i] ~ g(x[i]) where { meta = DeltaMeta(method = Linearization()) }
        y[i]  ~ Normal(mean = gk[i], var=y_s2)

        x_prev = x[i]
    end

end

and

constraints = @constraints begin
    q(x,gk,x_prev) = q(x,gk)q(x_prev)
    #q(mu_x,gk,xk,x_prev) = q(mu_x)q(x_prev)q(gk)q(xk)
end

inits = @initialization begin
    q(x) = vague(NormalMeanVariance)
    q(x_prev) = vague(NormalMeanVariance)
    q(gk) = vague(NormalMeanVariance)
    #q(xk) = vague(NormalMeanVariance)
    μ(x) = vague(NormalMeanVariance)
end

result = infer(
    model = Hydar_model(x_init=x_init, mu_v=mu_v, dt=dt, x_s2=1/gm_x_p, y_s2=1/gm_y_p),
    data = (y = y,),
    free_energy = false,
    constraints=constraints,
    initialization=inits,
    iterations=10,
    showprogress=true,
    returnvars=(x = KeepLast(),),
    options = (limit_stack_depth = 100,),
);

Next I simplify the issue and first try a non-linear model.

[wmkouw]

Hi Charel, I promised to post a bug report of the error you encountered. However, I can't seem to reproduce it on my end. The code below runs without errors.

using RxInfer
using Distributions

g(x) = 25 - 16 / (1 + exp(5 - x / 5))

@model function mwe(y,u,prior, s2_x, s2_y)

    x_prev ~ prior
    for i in 1:length(y)
        
        x[i]   ~ Normal(mean=x_prev-u, var=s2_x)
        gk[i]  ~ g(x[i]) where { meta = DeltaMeta(method = Linearization()) }
        y[i]   ~ Normal(mean=gk[i], var=s2_y)

        x_prev = x[i]
    end
end

x0 = Normal()
u  = 1.0
y  = randn(3)

cs = @constraints begin
    q(x_prev,x,gk) = q(x_prev)q(x,gk)
end

is = @initialization begin
    q(x_prev) = vague(NormalMeanVariance)
    q(x)      = vague(NormalMeanVariance)
    q(gk)     = vague(NormalMeanVariance)
    μ(x)      = vague(NormalMeanVariance)
end

res = infer(
    model          = mwe(prior=x0, u=u),
    data           = (y=y,),
    free_energy    = true,
    constraints    = cs,
    initialization = is,
    iterations     = 10,
    returnvars     = (x = KeepLast(),),
)

I tried a few variants (e.g., other constraints or initialization) but it always runs for free_energy=true. Let's have another look at this error next session.

In the meantime, I will compare Friston's solution with our message passing solution and see if I find any discrepancies that explain the difference in prior variances and estimated state means.

[wmkouw]

I have an hypothesis: the difference is due to filtering vs smoothing. Can you post the generative process where you simulate data? Then I can test this hypothesis.

[Charelvanhoof]

In this basic/simple Hydar is not (yet) moving, so it has the same temperature reading every cycle, probably not the rootcause.
in code:

execute_world, observe_world,f,g = initializeWorld(x_p=gp_x_p, y_p=gp_y_p, x_init=gp_x_init, dt=dt);
y = [g(gp_x_init, 0) for _ in 1:N]
# y = [g(rand(Normal(gp_x_init,1/gp_y_p)), 0) for _ in 1:N] variant with added sensory noise

code of the generative process::

# The generative process " the world"
function initializeWorld(; x_p, y_p, x_init, dt)
    # input parameters
    # x_p = precision of the generative process noise on the function of motion
    # y_p =  precision of the generative process noise on the function of sensory mapping
    #x_init
        
    # generative process
    f_motion(x,v,u) = u
    g_temp(x,v) = 25 - 16 / (1 + exp(5 - x / 5))

    x_t_last = x_init 
    y_t_temp = g_temp(x_init,0)
    
    function execute(a_t::Float64)
        # Original python code
        # x_dot= self.hydar.f(self.x[i],self.v, action) + np.random.randn(1) * np.sqrt(self.Sigma2_x) where f= 0*x + 1*u
        # self.x[i+1]= self.x[i] + self.dt*x_dot
        # self.y_temp[i+1] = self.hydar.g_temp(self.x[i+1],self.v) + np.random.randn(1) * np.sqrt(self.Sigma2_y_temp) where g=25 -16 / (1 + np.exp(5-x/5))
        x_t_w = rand(Normal(0, sqrt(1/x_p)))
        y_t_w = rand(Normal(0, sqrt(1/y_p)))
        x_dot = f_motion(0,0,a_t) + x_t_w
        x_t = x_t_last + dt*x_dot 
        x_t_last = x_t # Memory for next step
        y_t_temp = g_temp(x_t,0) + y_t_w
                
        return y_t_temp
    end

    observe() = y_t_temp # Temperature is observed

    return (execute, observe, f_motion, g_temp)
end;

[Charelvanhoof]

Tested it with a linear model, same results.

# Linear test

dt = 0.005; # integration time step, average neuron resets 200 times per second
T = 2+dt; # maximum time considered
N = Int(round(T/dt));
println("Amount of steps: ",N);

@model function Hydar_model(y, x_init, mu_v, dt, x_s2, y_s2)

    x_prev ~ x_init

    for i in 1:length(y)
        x[i] ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=x_s2)      # x_t+1=x_t + dt*f(x_t, v_t) = xt+1 = xt -dt*x_t + dt*v_t
        y[i]  ~ Normal(mean = x[i], var=y_s2)
        x_prev = x[i]
    end

end

# experiment input parameters
# generative process, actual world
gp_x_init = 30.0   # actual depth of hydar

# generative model, belief of Hydar

gm_x_prior = 25.0  # Prior belief of Hydar
# gm_x_init = gm_x_prior Hyadr initial belief of depth is the same as the prior
gm_x_p = 1.0       # Hyadrs' belief of precision of function of motion, trust in internal model
gm_y_p = 1.0       # Hydars' belief of precision of function of sensory observations, trust in sensory observations


x_init = Normal(gm_x_prior , 1/gm_x_p); # belief of Hydar of initial depth (mean, variance)
y = [gp_x_init for _ in 1:N] # sensory readings, all the same = actual depth/temperature of hydar

constraints = @constraints begin
    q(x,x_prev) = q(x)q(x_prev)
end

inits = @initialization begin
    q(x) = vague(NormalMeanVariance)
    q(x_prev) = vague(NormalMeanVariance)
    μ(x) = vague(NormalMeanVariance)
end

result = infer(
    model = Hydar_model(x_init=x_init, mu_v=gm_x_prior, dt=dt, x_s2=1/gm_x_p, y_s2=1/gm_y_p),
    data = (y = y,),
    free_energy = true,
    constraints=constraints,
    initialization=inits,
    iterations=10,
    showprogress=true,
    returnvars=(x = KeepLast(),),
    options = (limit_stack_depth = 100,),
);

xmarginals  = result.posteriors[:x]
convert.(NormalMeanVariance,xmarginals)

# print graph
mu_x= [mean(m) for m in xmarginals]
prior = fill(gm_x_prior, N)
TL = 1:N
plot(TL, y, lw = 2,xlabel = "Time", ylabel = "Depth/temp", label="generative process (hidden state)", color = :orange)
plot!(TL, prior, lw = 2, label="prior belief mu_v", color = :black)
plot!(TL, mu_x, lw = 2, label="belief mu_x", color = :blue)

Result:

Expected result with python active inference code (same prior, init, precisions, learning rates 1, embedding order: 1 derivative, etc):

Heart of the Python active inference code:

def ai_step (self, i, mu_v, y):
        """
        Perform active inference step   

        INPUTS:
            i       - tic, timestamp
            mu_v    - Hierarchical prior input signal (mean) at timestamp,  in generalized coordinates of motion
            y       - sensory input signal at timestamp, in generalized coordinates of motion

        INTERNAL:
            mu_x      - Belief or hidden state estimation, in generalized coordinates of motion

        """

        #-------------------------------#
        #                               #
        #           Prediction          #
        #                               #
        #-------------------------------#
        
        # Note that the predictions are in generalized coordinates of motion  
        mu_x_hat = self.std_vec*self.f(self.mu_x[0],mu_v[0]) + self.std_vec_inv * self.df(self.mu_x[0],mu_v[0]) * self.mu_x
        mu_y_hat = self.std_vec*self.g(self.mu_x[0],mu_v[0]) + self.std_vec_inv * self.dg(self.mu_x[0],mu_v[0]) * self.mu_x
        
        #-------------------------------#
        #                               #
        #        Prediction Error       #
        #                               #
        #-------------------------------#

        # Note that the predictions erros are in generalized coordinates of motion 
        
        eps_x = self.D.dot(self.mu_x) - mu_x_hat  # prediction error hidden state
        eps_y = y - mu_y_hat #prediction error sensory observation  
        
        # In case a linear state space model is used the below calculation is equivalent
        #eps_x = (self.D-self.Atilde).dot(self.mu_x) - mu_v  # prediction error hidden state
        #eps_y = y - self.Ctilde.dot(self.mu_x) #prediction error sensory observation

        #-------------------------------#
        #                               #
        # Prediction Error minimization #
        #                               #
        #-------------------------------#            
     
        # Calculate Free Energy to report out
        F= 0.5*( eps_x.T.dot(self.Pi_w).dot(eps_x) + eps_y.T.dot(self.Pi_z).dot(eps_y)).item(0) 
        # Note, the item(0) is needed to translate the python matrix result to a scaler, e.g. [[1]] to 1
        
        # Gradient descent inference/perception 
        Atilde=self.a * self.I
        dFdmu_x = (self.D-Atilde).T.dot(self.Pi_w).dot(eps_x) - (self.dg(self.mu_x[0],mu_v[0]) * self.I ).T.dot(self.Pi_z).dot(eps_y)
        dmu_x = np.dot(self.D,self.mu_x) - self.alpha_mu * dFdmu_x  
        self.mu_x = self.mu_x + self.dt * dmu_x
        
        # In case a linear state space model is used with linear equation of sensory mapping the below calculation is equivalent
        #Ctilde=self.c * self.I
        #dFdmu_x = (self.D-Atilde).T.dot(self.Pi_w).dot(self.eps_x) - Ctilde.T.dot(self.Pi_z).dot(self.eps_y)
        
        # Gradient descent action 
        if i>self.actiontime:
            fwd=fwd_model(self.p, self.mu_x[0], mu_v[0], self.u, self.fwd_method)          
            dFdu = fwd.T.dot(self.Pi_z).dot(eps_y).item(0)
            du = -self.alpha_u * dFdu  
            self.u = self.u + self.dt * du
        else:
            self.u=0
        

        return self.u, F, self.mu_x[0] , self.g(self.mu_x[0],mu_v[0])

[Charelvanhoof]

Basic active inference code in a diagram (in generalised coordinates of motion)

[wmkouw]

I've been testing filtering vs smoothing a bit, and I don't think that's the issue necessarily. However, my current implementation produces ...

Okay, so that didn't work. As discussed today, the difference between RxInfer and Friston's code may have to do with the time discretization of the generative model. For the state transition, we have

$$\dot{\mu}_x = f(\mu_x, \mu_v) + w$$

where, for now, $$f(\mu_x, \mu_v) = -\mu_x + \mu_v$$ and $$w$$ is a Wiener process with variance $$\sigma_x^2$$.

I double-checked the steps today and realized a mistake. For an approximate (first-order) discretization of the state SDE, the discretization should be:

$$\frac{\mu_{x,t+1} - \mu_{x,t}}{\Delta t} = -\mu_{x,t} + \mu_v + w_t$$ $$\mu_{x,t+1} = (1-\Delta t)\mu_{x,t} + \Delta t \mu_v + \Delta w_t$$

where $$\Delta w_t$$ is a Gaussian distribution with variance $$\sigma^2_x \Delta t$$. This is the Euler-Marayuma method and is accurate for small $$\Delta t$$. I actually think it's possible to perform an exact discretization here but I need some time for that.

Anyway, if I change the state transition in the batch inference (aka smoothing) model, I get

It's odd that it jumps up so quickly instead of the graceful exponential decay in Friston's code. But it does reach the $$\pm$$24 value that Friston's code levels off to (the tail at the end has to do with the batch inference process). There must be another discrepancy between the two codebases. I will investigate further.

Code:

"""Generative process"""

dt = 0.005
T = 10+dt
tsteps = range(start=dt,step=dt,stop=T)
len_time = length(tsteps)

# f(mx::Real, mv::Real) = -mx + mv
# f_step(mx::Real, mv::Real) = mx + dt*f(mx, mv)
# g(x) = 25 - 16 / (1 + exp(5 - x / 5))

gp_x_init = 25.0   # actual depth of hydar
gp_x_p = 1.0       # actual precision of function of motion
gp_y_p = 1.0       # actual precision of function of sensory observations
gm_x_init = 15.0   # Hyadr's belief of depth
gm_x_target = 15.0 # Prior belief of Hydar
gm_x_p = 1.0       # Hydars'belief of precision of function of motion
gm_y_p = 1.0       # Hydars belief of precision of function of sensory observations
v_t = 0.0

# The generative process " the world"
function initializeWorld(; x_p, y_p, x_init, dt)
    # input parameters
    # x_p = precision of the generative process noise on the function of motion
    # y_p =  precision of the generative process noise on the function of sensory mapping
    #x_init
        
    # generative process
    f_motion(x,v,u) = u
    g_temp(x,v) = 25 - 16 / (1 + exp(5 - x / 5))

    x_t_last = x_init 
    y_t_temp = g_temp(x_init,0)
    
    function execute(a_t::Float64)
        # Original python code
        # x_dot= self.hydar.f(self.x[i],self.v, action) + np.random.randn(1) * np.sqrt(self.Sigma2_x) where f= 0*x + 1*u
        # self.x[i+1]= self.x[i] + self.dt*x_dot
        # self.y_temp[i+1] = self.hydar.g_temp(self.x[i+1],self.v) + np.random.randn(1) * np.sqrt(self.Sigma2_y_temp) where g=25 -16 / (1 + np.exp(5-x/5))
        x_t_w = rand(Normal(0, sqrt(1/x_p)))
        y_t_w = rand(Normal(0, sqrt(1/y_p)))
        x_dot = f_motion(0,0,a_t) + x_t_w
        x_t = x_t_last + dt*x_dot 
        x_t_last = x_t # Memory for next step
        y_t_temp = g_temp(x_t,0) + y_t_w
                
        return y_t_temp
    end

    observe() = y_t_temp # Temperature is observed

    return (execute, observe, f_motion, g_temp)
end;

execute_world, observe_world,f,g = initializeWorld(x_p=gp_x_p, y_p=gp_y_p, x_init=gp_x_init, dt=dt);
y = [g(gp_x_init, 0) for _ in 1:len_time]
# y = [g(rand(Normal(gp_x_init,1/gp_y_p)), 0) for _ in 1:len_time] #variant with added sensory noise

"""Model specification"""


@model function Hydar_model(y, x_init, mu_v, dt, x_s2, y_s2)

    x_prev ~ x_init

    for i in 1:length(y)
        x[i] ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=dt*x_s2) # x_t+1=x_t + dt*f(x_t, v_t) = xt+1 = xt -dt*x_t + dt*v_t
        # tmp ~ mu_v - x_prev
        # x[i] ~ Normal(mean= tmp , var=x_s2) 
        #xk[i] ~ f_step(x_prev,mu_v) where { meta = DeltaMeta(method = Linearization()) }
        #x[i]  ~ Normal(mean = xk[i], var=x_s2) 
        gk[i] ~ g(x[i],0) where { meta = DeltaMeta(method = Linearization()) }
        y[i]  ~ Normal(mean = gk[i], var=y_s2)

        x_prev = x[i]
    end

end

"""Inference"""

x_init = NormalMeanVariance(gm_x_init, 1. / gm_x_p)
mu_v = 1.0

constraints = @constraints begin
    q(x,gk,x_prev) = q(x,gk,x_prev)
    #q(mu_x,gk,xk,x_prev) = q(mu_x)q(x_prev)q(gk)q(xk)
end

inits = @initialization begin
    q(x) = vague(NormalMeanVariance)
    q(x_prev) = vague(NormalMeanVariance)
    q(gk) = vague(NormalMeanVariance)
    #q(xk) = vague(NormalMeanVariance)
    μ(x) = vague(NormalMeanVariance)
end

result = infer(
    model = Hydar_model(x_init=x_init, mu_v=mu_v, dt=dt, x_s2=1/gm_x_p, y_s2=1/gm_y_p),
    data = (y = y,),
    free_energy = false,
    constraints=constraints,
    initialization=inits,
    iterations=10,
    showprogress=true,
    returnvars=(x = KeepLast(),),
    options = (limit_stack_depth = 100,),
)

estimates = mean.(result.posteriors[:x])
var_state = var.(result.posteriors[:x])

plot(xlabel="time", ylabel="depth")
hline!([gp_x_init], label="Generative process", color="orange")
# scatter!(tsteps, y, label="Observations")
plot!([0.; tsteps], [gm_x_init; estimates], label="State belief", color="blue")
hline!([gm_x_target], label="Prior belief", color="black")

[Charelvanhoof]

Many thanks, investigating the difference between RxInfer and Friston's code (time discretization) is definitively the right track. Your logic
$\dot \mu_x(t) = -\mu_x(t) +\mu_v(t) + w(t)$ with $\mu_x(t+1)= \mu_x(t)+dt*\dot\mu_x(t)$ results in $\mu_x(t+1)= (1-dt)\mu_x(t)+dt*\mu_v + dt*w(t)$ makes absolute sense. I think the solution direction goes back to the earlier post of Bert, let's discuss the afternoon.

The Code above contains an copy-paste error, the "mu_v = 1.0" appears and that sets the prior to 1.0. If you remove it and run the model with prior=gm_x_target with "model = Hydar_model(x_init=x_init, mu_v=gm_x_target, dt=dt, x_s2=1/gm_x_p, y_s2=1/gm_y_p)," you get rid of the first suspicious bump.

next i corrected the precision of the Hydar's initial belief of the depth as well x_init = NormalMeanVariance(gm_x_init, dt*1. / gm_x_p):

even a better curve. However still not yet getting the same results as Friston. Best illustrated by running the linear batch example with "x[i] ~ Normal(mean= (1-dt)x_prev + dtmu_v , var=dt*x_s2) " I get:

while the expected result is (exactly in the middle, both precisions 1):

[Charelvanhoof]

I can confirm "the tail at the end has to do with the batch inference process", below the linear recursive graph
linear example:

And some testing to find the delta to get close the the Friston graphs, around a factor 1/70 (maybe gives a clue where to find the delta):

factor=1/70

@model function Hydar_model_linear_recursive(y, x_prev_p, mu_v, dt, x_s2, y_s2)

    x_prev ~ x_prev_p

    x ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=dt*factor*x_s2)
    y ~ Normal(mean= x, var=y_s2) 

end

resulting in:

[wmkouw]

Ok, so I compared the two estimators (as currently specified) for the linear model of generalized coordinate depth 1 by hand.

Gradient descent on Free Energy (GDFE)

The generative model for a single generalized coordinate is specified to be:

$$\begin{align} \mu_x' &= f(\mu_x, \mu_v) + w \\ y &= g(\mu_x, \mu_v) + z \end{align}$$

where $\mu_x$ is the state, $\mu_x'$ is the time-derivative of the state, $\mu_x$, $\mu_v$ is the velocity input, $w \sim \mathcal{N}(0, \sigma_x^2)$ and $z \sim \mathcal{N}(0, \sigma_y^2)$. For the linear case and depth 1, $f(\mu_x, \mu_v) = -\mu_x + \mu_v$, and $g(\mu_x, \mu_v) = \mu_x$.

The free energy under this model, with $\Pi_x = \Pi_y = 1$, is:

$$\begin{align} \mathcal{F}[y,\mu_x, \mu_v] &= \frac{1}{2}\big( (\mu_x' - f(\mu_x, \mu_v))^{\top} \Pi_x (\mu_x' - f(\mu_x, \mu_v)) + (y - g(\mu_x, \mu_v)^{\top} \Pi_y (y - g(\mu_x, \mu_v))) \big) \\ &= \frac{1}{2}\big( (\mu_x' + \mu_x - \mu_v)^2 + (y - \mu_x)^2 \big) \end{align}$$

The gradient of the free energy with respect to $\mu_x$ is:

$$\begin{align} \frac{\partial}{\partial \mu_x} \mathcal{F}[y,\mu_x, \mu_v] = (\mu_x' + \mu_x - \mu_v) + (y - \mu_x) \end{align}$$

The change in the state estimate is given as:

$$\begin{align} \dot{\mu}_x &= \mu_x' - \frac{\partial}{\partial \mu_x} \mathcal{F}[y,\mu_x, \mu_v] \\ &= -\mu_v + y \end{align}$$

If we then discretize this continuous-time change $\dot{\mu}_x$ with a Forward Euler, we get:

$$\begin{align} \mu_x(t+\Delta t) = \mu_x(t) + \Delta t( -\mu_v + y) \end{align}$$

Note: I sense a mistake somewhere because this equation would diverge over time. And your simulation converges.

Message passing (MP)

In this case, the continuous-time state transition has to be discretized at the start. With Euler-Marayuma, this is:

$$\begin{align} \mu_x' &= -\mu_x + \mu_v + w \implies \mu_{x,t+1} = \mu_{x,t} + \Delta t (-\mu_{x,t} + \mu_v + w) \end{align}$$

which implies that the discrete-time stochastic state transition is $p(\mu_{x,t+1} \mid \mu_{x,t+1}) = \mathcal{N}( \mu_{x,t} \mid (1-\Delta t)\mu_{x,t} + \Delta t \mu_v, \Delta t^2 \sigma_x^2 )$.

Thus, we have the generative model:

$$\begin{align} p(\mu_{x,t-1}) &= \mathcal{N}(\mu_{x,t-1} \mid m_{x,t-1}, s^2_{x,t-1}) \\ p(\mu_{x,t} | \mu_{x,t-1}) &= \mathcal{N}(\mu_{x,t} \mid (1-\Delta t)\mu_{x,t-1} + \Delta t \mu_v, \Delta t^2 \sigma_x^2 ) \\ p(y_t | \mu_{x,t}) &= \mathcal{N}(y_t \mid \mu_{x,t}, \sigma_y^2) \end{align}$$

with the accompanying recursive factor graph:

              \mu_v
\mu_x,t-1       |     \mu_x,t
      --------[ N ]----[=]---- 
                        |
                      [ N ]
                        |
                       y_t

There are 4 messages:

0. The message from the state prior.
1. The prediction message from the state transition node.
2. The correction message from the likelihood node.
3. The state posterior node (which will serve as the prior for the next time-step).

Message 0 is the state prior itself:

$$\nu_0(\mu_{x,t-1}) = \mathcal{N}(\mu_{x,t-1} \mid m_{x,t-1}, s^2_{x,t-1}) $$

Message 1 is the integral of the state transition over the incoming message 0:

$$\begin{align} \nu_1(\mu_{x,t}) &= \int p(\mu_{x,t} | \mu_{x,t-1}) \nu_0(\mu_{x,t-1}) \mathrm{d}\mu_{x,t-1} \\ &= \int \mathcal{N}(\mu_{x,t} \mid (1-\Delta t)\mu_{x,t-1} + \Delta t \mu_v, \Delta t^2 \sigma_x^2 ) \mathcal{N}(\mu_{x,t-1} \mid m_{x,t-1}, s^2_{x,t-1}) \mathrm{d}\mu_{x,t-1} \\ &\propto \mathcal{N}(\mu_{x,t} \mid (1-\Delta t)m_{x,t-1} + \Delta t \mu_v, \Delta t^2 \sigma_x^2 + s^2_{x,t-1}) \end{align}$$

Message 2 is the integral of the likelihood node over the observed data point:

$$\begin{align} \nu_2(\mu_{x,t}) &= \int \mathcal{N}(y_t \mid \mu_{x,t}, \sigma_y^2) \delta(y_t - \hat{y}_t) dy \\ &= \mathcal{N}( \mu_{x,t} \mid \hat{y}_t, \sigma_y^2) \end{align}$$

Message 3 is the product of messages 1 and 2:

$$\begin{align} \nu_3(\mu_{x,t}) &= \nu_1(\mu_{x,t}) \cdot \nu_2(\mu_{x,t}) \\ &= \mathcal{N}(\mu_{x,t} \mid (1-\Delta t)m_{x,t-1} + \Delta t \mu_v, \ \Delta t^2 \sigma_x^2 + s^2_{x,t-1}) \cdot \mathcal{N}( \mu_{x,t} \mid \hat{y}_t, \sigma_y^2) \\ &= \mathcal{N}(\mu_{x,t} \mid m_{x,t}, s^2_{x,t}) \end{align}$$

where

$$\begin{align} s^2_{x,t} &= ( (\Delta t^2 \sigma_x^2 + s^2_{x,t-1})^{-1} + \sigma^{-2}_{y})^{-1} \\ m_{x,t} &= \frac{(\Delta t^2 \sigma_x^2 + s^2_{x,t-1})^{-1} ( (1-\Delta t)m_{x,t-1} + \Delta t \mu_v) + \sigma_y^{-2} \hat{y}_t }{(\Delta t^2 \sigma_x^2 + s^2_{x,t-1})^{-1} + \sigma^{-2}_{y}} \\ &= ( (1-\Delta t)m_{x,t-1} + \Delta t \mu_v) + \underbrace{\frac{\sigma_y^{-2}}{(\Delta t^2 \sigma_x^2 + s^2_{x,t-1})^{-1} + \sigma^{-2}_{y}}}_{\text{Kalman gain}} \underbrace{(\sigma_y^{-2} \hat{y}_t - ( (1-\Delta t)m_{x,t-1} + \Delta t \mu_v))}_{\text{innovation}} \end{align}$$

Comparison

The two approaches lead to different results. But we can perhaps find a corresponding case. If we assume $\sigma_x^2 = \sigma_y^2 = 1$, then the mean estimate is:

$$\begin{align} m_{x,t} &= ( (1-\Delta t)m_{x,t-1} + \Delta t \mu_v) + \frac{1}{(\Delta t^2 + s^2_{x,t-1})^{-1} + 1} (\hat{y}_t - ( (1-\Delta t)m_{x,t-1} + \Delta t \mu_v)) \end{align}$$

Let's use the same notation on the GPFE estimate:

$$\begin{align} \mu_{x, t} = \mu_{x,t-1} + \Delta t(\hat{y}_t -\mu_v ) \end{align}$$

It seems to me that the time step sizes are in entirely different places, and that the GDFE method is missing the previous state uncertainty $s^2_{x,t-1}$. Mmh, I'm not sure we can reconcile this.

[Charelvanhoof]

Interesting puzzel indeed. I cannot look "under the hood" of RxInfer so i can't judge, but maybe I can help to make the problem even simpler, to try to find a hint. By going back to an even more simple case. Also no motion in the water. Hydar imply believes its depth is the prior/desired depth, no generalized coordinates
$\mu_x(t+\Delta t) =\mu_x(t) + \Delta t* \dot\mu_x$
$\mu_x(t+\Delta t) =\mu_x(t) + \Delta t*- \frac{\delta F}{\delta \ mu_x}$
$\mu_x(t+\Delta t) =\mu_x(t) - \Delta t*( \frac{\epsilon_x}{ \delta mu_x}\frac{1}{ \sigma_x^2}\epsilon_x+ \frac{\epsilon_y}{ \delta mu_x} \frac{1}{ \sigma_y^2}\epsilon_y)$
$\epsilon_x=\mu_x-\mu_v$ and $\epsilon_y=y-\mu_x$
$\mu_x(t+\Delta t) =\mu_x(t) - \Delta t*( \frac{1}{ \sigma_x^2}(\mu_x-\mu_v)- \frac{1}{ \sigma_y^2}(y-\mu_x))$
in python code

        # Calculate prediction errors
        self.eps_x = self.mu_x - mu_v  # prediction error hidden state
        self.eps_y = y - self.mu_x #prediction error sensory observation
        # Free energy gradient
        dFdmu_x = self.eps_x/self.Sigma_w - self.eps_y/self.Sigma_z
        # Perception dynamics
        dmu_x   = 0 - dFdmu_x  # Note that this is an example without generalised coordinates of motion hence u'=0
        # motion of mu_x 
        self.mu_x = self.mu_x + self.dt * dmu_x

with:
actual_depth x_init=30 (sensory signal y = 30)
prior mu_v= 25
precision prior = 1
precision sensor =1
resulting in posterior=27.5:

Also typical behaviour is if the precisions are lowered the graph takes longer to convert (eg wih both variances set on 10 = presion 0.1) the graph looks like this below. With the current rXInfer code ( x[i] ~ Normal(mean= (1-dt)x_prev + dtmu_v , var=dt*x_s2)) I dont get that effect, nearly the same graph compared to both precisions set on 1. Maybe that is also a relevant detail.

[wmkouw]

I think I have it now. As discussed Friday, the issue is probably that Friston is implicitly utilizing a mean-field factorization over states, $q(x) = \prod_{k} q(_k)$. The model in RxInfer has been using a structured factorization, $q(x) = q(x_0) \prod_k q(x_k \mid x_{k-1})$. If I re-derive the message passing rules by hand under a mean-field factorization for RxInfer, I get the update rules:

$$\begin{align} m_k &= \frac{1}{\frac{1}{\Delta t} + 1} \Big(\frac{1}{\Delta t} \big( (1-\Delta t)m_{k-1} + \Delta t \mu_v \big) + y_k \Big)\\ v_k &= \frac{1}{\frac{1}{\Delta t} + 1} \end{align}$$

where the variance is constant over states due to the mean-field. Running it for y_k = 30 (for all k), gm_x_init = 25 and \mu_v = 25, gives an estimate of 27.5:

The code has become a whole lot simpler for this as well. For the recursive model, it's not actually even necessary to apply constraints and initialization;

"""Model specification"""


@model function recursive_Hydar_model(y, x_prev, mu_v, dt, x_s2, y_s2)

    x    ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=dt*x_s2)
    y    ~ Normal(mean = x, var=y_s2)

end

"""Inference"""

gm_x_init = 25.0
gm_x_p = 1.0     
gm_y_p = 1.0     
mu_v = 25

estimates = zeros(len_time)
var_state = zeros(len_time)

x_prev = gm_x_init
for t in 1:len_time

    result = infer(
        model = recursive_Hydar_model(x_prev=x_prev, mu_v=mu_v, dt=dt, x_s2=1/gm_x_p, y_s2=1/gm_y_p),
        data = (y = y[t],),
    )
    x_prev = mean(result.posteriors[:x])

    estimates[t] = mean.(result.posteriors[:x])
    var_state[t] = var.(result.posteriors[:x])
end

plot(xlabel="time", ylabel="depth", size=(800,300))
hline!([gp_x_init], label="Generative process", color="orange")
plot!([0.; tsteps], [gm_x_init; estimates], label="State belief", color="blue")
hline!([gm_x_init], label="Prior belief", color="black")

savefig("figures/hydar_linear_filtering.png")

[Charelvanhoof]

We were so close the first time, interesting detour but super nice we got matching graphs.
I have been testing with the precisions.
Experiment setup:
gp_x_init = 30.0 # actual depth of hydar
gm_x_prior = 25.0 # Prior belief of Hydar
gm_x_init = 25.0 # Hyadr's initial belief of depth

Exactly the same results as Friston (with /without generalized coordinates) when changing the precision of the sensory prediction error:

• with gm_x_p = 1.0 and gm_y_p = 10.0 both Friston and rxinfer result in: 29.545454
• with gm_x_p = 1.0 and gm_y_p = 0.1 both Friston and rxinfer result in: 25.454537

The same results as Friston without generalized coordinates when changing the precision of the motion prediction error :

• with gm_x_p = 10.0 and gm_y_p = 1.0 rxinfer result in: 25.454537, Friston without generalized coordinates: 25.454545
• with gm_x_p = 0.1 and gm_y_p = 1.0 rxinfer result in: 29.545454, Friston without generalized coordinates: 29.54538
  Note that it are symmetric results.

However different results compared to Friston with generalized coordinates:

• with gm_x_p = 10.0 and gm_y_p = 1.0 Friston with generalized coordinates: 27.17475 (1 embedding order, approx same with 3 embedding orders)
• with gm_x_p = 0.1 and gm_y_p = 1.0 Friston with generalized coordinates: 28.08948 ( 1 embedding order, 28.6015 with 3 embedding orders)
  Friston with generalized coordinates leads to different results behaviour, but let's compare once we have also RxInfer expressed in generalized coordinates.

[Charelvanhoof]

Today I have been testing a non-linear model, and the short conclusion is that it works.

I did the same experiment as experiment 1.01 in notebook: https://www.kaggle.com/code/charel/active-inference-code-by-example-1
With the same settings, eg precisions 1, prior 25, init 30, Learning rate 1 (a learning rate is introduced in the typical continious active inference set-up), etc.

The RxInfer generative model

"""Generative model specification"""

@model function recursive_Hydar_model(y, x_prev, mu_v, dt, x_s2, y_s2)
    
    g = (x::Real) -> 25 - 16 / (1 + exp(5 - x / 5))

    x    ~ Normal(mean= (1-dt)*x_prev + dt*mu_v , var=dt*x_s2)
    gk   ~ g(x) where { meta = DeltaMeta(method = Linearization()) }
    y    ~ Normal(mean = gk, var=y_s2)

end

The RxInfer results:

The Frison/Python results:

Highly simular graphs. the Python code converges to μₓ last 5 values: [26.776059 26.776059 26.776059 26.776059 26.776059] while the RxIbnfer code to [26.776391, 26.776391, 26.776391, 26.776391, 26.776391, 26.776391]. Which can be easily explained because in the Python code the exact derivative of the non-linear function g has to be given as input to the model (return -16/5* np.exp(5-x/5) / (np.exp(5-x/5)+1)**2) while in RxInfer is is approximated (gk ~ g(x) where { meta = DeltaMeta(method = Linearization()) }) and it does not have be calculated/given upfront. RxInfer has done a remarkable good job to approximate the non-linear function.

Next I did some experiments with different precisions. Also here the results are highly comparable

• With high preicion of the sensory observation (10) compared to the generative model (1) the python code converges to μₓ trust sensor last 5 values: [29.083478 29.083478 29.083478 29.083478 29.083478] and the RxInfer code to 29.084549, 29.084549, 29.084549, 29.084549, 29.084549, 29.084549]
• With high precision of the generative model (10) compared to the sensory observations (1) the Python code converges to μₓ trust model 5 values: [25.277768 25.277768 25.277768 25.277768 25.277768] and the RxInfer code to [25.277764, 25.277764, 25.277764, 25.277764, 25.277764, 25.277764]

Giving me the confidence that this continuous active setup on RxInfer is working correctly (non-generalized coordinates, excluding action, etc)

[wmkouw]

Nice! That's excellent! Great to see.

So, the major pain point still is the generalized coordinates setup? In hindsight, the major issue was that RxInfer defaults to a structured factorization of the variational model (for Gaussians) while Friston assumes mean-field over states. It could be that there is still a mismatch between the factorization of each coordinate in the generalized coordinates vector between RxInfer and Friston. Do you happen to know what Friston assumes there? Because if you construct the generalized coordinates as a state vector in RxInfer, then I think it defaults to a structured factorization as well (I have to double-check).

[bvdmitri]

Great progress @Charelvanhoof ! A small tip: it is possible to provide the inverse function to the DeltaMeta, which might improve the accuracy even further (in the situations where you have access to the inverse of course).

g = (x::Real) -> 25 - 16 / (1 + exp(5 - x / 5))

# has singularities at y=9 and at y=25
inverse_g = (y::Real) -> 25 - 5 * log((y - 9) / (25 - y))

DeltaMeta(method = Linearization(), inverse = inverse_g)

In this particular case the inverse function has singularities, but still would be interesting to test. It should also in theory improve the performance.

[Charelvanhoof]

Started to program the support functions for colored noise, eg the typical precision matrices in case of colored noise.
Also a first shot from a mathematical perspective:

Same formulation as without generalized coordinates, but now on each level of the generalized coordinates.

However I was expecting to see the different levels of the generalized coordinates influence each other, like for example this equation when friston claculates the motion. $\dot{\tilde{\mu}}_x = D \tilde\mu_x - \partial{\tilde{\mu}_x} F(\tilde{y}, \mu)$

Now it looks like that each level (position, speed, acceleration, etc) is estimated independently of each other (or is that done "under the hood" in RxInfer)

[bvdmitri]

Could you provide the model specification with the model macro? RxInfer does not assume independence "under the hood" unless its explicitly specified in the model macro