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.

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.

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).

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.

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

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.