How to compute simulations efficiently?

Modeling
1

Hi all,

I got a (maybe stupid) question about how to run simulate efficiently…

I have written a model with the structure

functions {
real[,] run_simulation(argument1; argument2) {
 [....]
return(results);
}

data {
int N1;
int N2;
real experimental_results[N1,N2];
}

parameters{
real<lower=0> parameter_1;
real<lower=0> parameter_2;
real<lower=0> sigma;
}

transformed parameters {
real simulations[N1, N2];
simulation = run_simulation(parameter_1, parameter_2);
}

model {
//priors
parameter_1 ~ lognormal (1,2);
parameter_2 ~ lognormal (1,2);

//likelihood
for (i in 1:N1) {
    for (j in 1:N2) {
        experimental_results[i,j] ~ normal (simulation[i,j] + sigma )
   }
}

I’m wondering if the way of writing the likelyhood in for loops like this

//likelihood
for (i in 1:N1) {
    for (j in 1:N2) {
        experimental_results[i,j] ~ normal (simulation[i,j] + sigma )
   }
}

Does this mean that my run_simulation is run N1*N2 times in each iteration, or is it run once? Are there better ways of writing this?
Thanks!

2

Because you have run_simulation evaluated in the transformed parameters block it will be run every step of the leapfrog integrator (many times per iteration of the Markov chain), but it is not run additional times because of your loop in the model block (the loop is just making use of the result computed in transformed parameters).

To speed it up you can use vectorization:

to_vector(experimental_results) ~ normal (to_vector(simulation), sigma);

I assume this is a typo and you meant ~ normal (simulation[i,j], sigma ), i.e comma not plus, right?

3

Actually since you’re using 2-D real arrays instead of matrices you’d need to_array_1d instead of to_vector. If you change them to matrices you could use to_vector.