Define distribution functional constraints

[mssaw]

My model describes ODE system (more than 30 eqs.) defined by 12 parameters.
I want to sample 8 parameters (say theta) where other 4 (say phi) are defined by functional constraints which depend on theta parameters. All parameters have to be positive as they represent some biological processes. Some of them have to by less then 1.

For simplicity say I have following model:

functions{
real fun(real a, real b){
 real c;
 ...some complicated algebraic expressions depending on a and b
 return c; 
}
}
data{
int NoObs;
real[NoObs] d;
}
parameters{
real a;
real b;
}
transformed parameters{
real y[NoObs];
real par[3];
real<lower = 0 > c;

c = fun(a,b);
par[1] = a;
par[2] = b;
par[3] = c;

y = integrate_ode_rk45(some_ode_model, ..., par )
}
model {
  a ~ normal(0,1) T[0, ];
  b ~ normal(0,1) T[0, ];

  d ~ normal(y, 1); 
} 

I would like to sample in model block from the space where my c=fun(a,b)>0.
Is it possible to do something like that

model {
  a ~ normal(0,1) T[0, ];
  b ~ normal(0,1) T[0, ];
  c ~ fun(a,b) T[0, ];

  d ~ normal(y, 1); 
} 

and treat fun(a,b) as fake distribution returning fun(a,b) value? I just want to operate into space where c>0. Otherwise sampler is taking ages (I mean 2 weeks for 2000 chain) to run even very few iterations (days) and have problems with convergence.

Thanks in advance for any suggestions.

[wds15]

First thing is to declare Parameter a and b with lower limits equal 0 instead of using truncation as you do.

Next you should probably do the declaration of c with the lower limit as you do in the transformed parameters block, but move the integrate call into the model block. I think that that the constraints get checked at the end of a block. Alternatively, you could include an if statement which rejects the draw whenever c is not greater than 0. You don’t want to integrate when c is off here.

[andre.pfeuffer]

All good @wds15, except for the if statement, I’d suggest to change fun by return exp(c)
instead of c.

real fun(real a, real b){
 real c;
 ...some complicated algebraic expressions depending on a and b
 return exp(c); 
}

or write:

log(c) ~ fun(a, b);
target += -log(c);

[wds15]

Solving odes is “evil” in stan…if statements can even be worse, but don’t need to be if the typical set is not running often against the constraint limit. In any case, you are right in that avoiding a hard if statement is rather important for good performance of nuts. So expressing the function for the log of c and then exponentiating this would be the best way (or any other way which keeps c positive by construction).

[nhuurre]

@mssaw

I want to sample 8 parameters (say theta) where other 4 (say phi) are defined by functional constraints which depend on theta parameters. All parameters have to be positive as they represent some biological processes. Some of them have to by less then 1.

So your parameter space is an 8-dimensional submanifold of 12-dimensional space, defined by 4 equations and more than 12 inequalities. That’s sounds quite difficult geometry. There’s no general solution, we’ll need to see the exact functional constraints.

But regardless of that, ODE system with more than 30 equations is going to perform super slow no matter how it’s parametrized. See this earlier thread.

Costs/benefits of solving large ODE systems as systems of difference equations/in discrete time

I have been working with large systems of ordinary differential equations in a systems modelling context, similar to the ‘system dynamics’ literature of Jay Forrester, John Sterman et al. System dynamics models often include a large number of ODEs (20+ is not uncommon) and parameters, and data for some of the state variables may never be available/have a high proportion of missingness. For example, see the C-ROADS model (https://www.climateinteractive.org/tools/c-roads/). Solving large systems …

[mssaw]

@andre.pfeuffer

I’d suggest to change fun by return exp(c)
instead of c .

What I would gain doing that?

[mssaw]

@andre.pfeuffer

log(c) ~ fun(a, b); target += -log(c);

Also this does not work as fun is not a distribution.

[andre.pfeuffer]

@mssaw

Also this does not work as fun is not a distribution.

@mssaw

I would like to sample in model block from the space where my c=fun(a,b)>0.

This confuses me. How to sample from a non-distribution? (No offense, I don’t understand.)

@mssaw

I’d suggest to change fun by return exp(c)
instead of c .

What I would gain doing that?

You asked for a positive c , applying exp make it positive.
You have to convert the constraint space to the unconstrained either explicit or implicit, because Stan works on the unconstrained space. There is an explanation in the Stan manual.