Pystan with a user defined function (model)

Dear Stan community, hello!
I have a python agent based model named “flee” that I’m using to simulate the movement of agents from different nodes. I run the flee model for a specific scenario from the terminal via the following command:
fabsim localhost flee: “simulation name”,simulation_period=“number of time steps”

The agents have a specified “movespeed”, “awarelev” (awareness of the # of next nodes) and “movechance” (probability of moving from one node to another according to the node’s classification, there are 3 classifications) - all are specified in input files read by flee.

I would like to link up the flee model with Stan (pystan in python) to explore the parameter space of movespeed, awarelev and movechance. I thought of applying flee as a user function to Stan. I am VERY NEW to Stan so I’m not sure if the syntex below is correct, but most importantly, I am not sure how to perform the link between flee and Stan.

N is the number of end nodes
t is the number of time steps
x is the number of agents in the origin nodes
yhat is the number of agend in the end nodes

model = “”"
data {
int<lower=0> N;
int<lower=0> t;
matrix[N,t] x;
matrix[N,t] yhat;
parameters {
real movespeed;
real awarelev;
vector[3] movechance;
model {
yhat ~ flee(movespeed, awarelev, movechance, x)

I am stumbling around in the dark here, hope you guys can help.
Happy to provide more information!

unfortunately, Stan has quite strong restrictions on what functions can be connected to it. The more problematic is that you need to be able to evaluate not only the value of the function, but also the gradient of the function with respect to all of the parameters. To use it as a probability distribution, you need to return the density, not a stochastic outcome. The second is that you need to call the function from C++.

I fear that neither of those are satisfied in your use case and while having a C++ wrapper might end up feasible, the gradient part is usually not achievable for agent-based simulations.

Inference for your model might be possible with approximate Bayesian computation.

Unfortunately we can’t really provide support for ABC here as it is out of scope of this forum, but I’ll quote @avehtari who previously advised this:

If you think your case fits ABC framework, then I recommend using ELFI with BOLFI (and HMC/NUTS)
If need more information you can ask at

If the function is too complicated to program as a Stan function, you can plumb it in through C++. If it (a) is fully templated, and (b) uses Eigen and std::vector types, it’ll map directly into Stan types and be fairly straightforward. If not, you need gradients with respect to all the input parameters. Either way, the top-level overview is in:

and there are more details in

We’re still consolidating our wiki pages!