I’m not sure whether I should make a new thread, given that the initial problem was “solved” (in the sense that it is not expected to use time that depend on the parameters as I was trying to do).
Please let me know if you think that I should make the previous comment by yizhang as a solution as well.
Anyways, I’m trying the approach you mentioned earlier to turn this problem into a canonical ODE, but the code is giving me an unknown error (literally) and I have no idea on what could be wrong (I’ve double checked my math and my code).
The model
First of all, I’ve omitted minor things during the previous description of my model, and I shall make the model clear here.
In my use case, the parameters in my system are
\theta = (h, \Omega_b, \Omega_c, \Omega_r)
and a quantity \lambda which can be derived from the parameters which is given by
\lambda = \frac{1}{2} + W_0\left( -\frac{\Omega_b + \Omega_c + \Omega_r}{2e^{1/2}} \right)
where W_0 is the 0-th branch of the Lambert function.
As for the data, the following two quantities will be used as datasets:
R = \sqrt{\Omega_b + \Omega_c} \int_0^{z_*} u(z) dx
w_b = \Omega_b*h^2
where u(z) \equiv 1/E(z), as you defined previously and z_* depends only on the parameters (the expression to compute its value from the parameters is fairly complicated and does not bring any meaningful information).
The experiment gave us that
R = 1.710 \pm 0.019
w_b = 0.0224 \pm 0.0001
Using the equations you have so kindly provided I obtained:
u' = \frac{e^{-\lambda u^2}}{\lambda u (u^{-2} - 2 \lambda) - u^{-3}} \left[\frac{3}{2}(\Omega_b + \Omega_c)(1+z)^2 + 2 \Omega_r(1+z)^3) \right]
To compute the vale of R I solve the previous differential equation using one of Stan’s ODE solver, to obtain u(z), and integrate it using one of Stan’s integration routines, to finally obtain \int_0^{z_*} u(z) dx.
Stan model file
/* block for user defined functions */
functions {
// auxiliary function to solve for u(z) ≡ 1/E(z) using ODE solver
vector ode(real z, vector y, array[] real faketheta) {
real Omega_b = faketheta[2];
real Omega_c = faketheta[3];
real Omega_r = faketheta[4];
real lambda = faketheta[5];
real u = y[1];
real f1 = exp(-lambda*u^2) / (lambda*u*(u^(-2)-2*lambda) - u^(-3));
real f2 = 1.5*(Omega_b + Omega_c)*(1+z)^2 + 2*Omega_r*(1+z)^3;
vector[1] output;
output[1] = f1*f2;
return output;
}
// returns u(z) ≡ 1/E(z) to integrate
real u(real z, real xc, array[] real faketheta, array[] real x_r, array[] int x_i) {
// initial conditions
real z0 = 0;
vector[1] init;
init[1] = 1;
// obtain u(z) using ODE solver
vector[1] vec = ode_rk45(ode, init, z0, {z}, faketheta)[1];
real sol = vec[1];
return sol;
}
}
/* block for datasets used */
data {
real Rexp;
real Rexp_err;
real wbexp;
real wbexp_err;
}
/* block for transformed data */
transformed data {
// create null data values to give to integrate_1d because they are both required arguments
array[0] real x_r;
array[0] int x_i;
}
/* block for model parameters */
parameters {
real<lower=0> h;
real<lower=0> Omega_b;
real<lower=0> Omega_c;
real Omega_r;
}
/* block for transformed parameters */
transformed parameters {
// compute λ
real lambda = 0.5 + lambert_w0( -(Omega_b + Omega_c + Omega_r)/(2*exp(0.5)) );
// parameter array (labeled as fake because it includes λ, which is not a true parameter)
array[5] real faketheta = {h, Omega_b, Omega_c, Omega_r, lambda};
// compute z_star
real g1 = 0.0783*(Omega_b*h^2)^(-0.238)/(1 + 39.5*(Omega_b*h^2)^(-0.763));
real g2 = 0.560/(1 + 21.1*(Omega_b*h^2)^1.81);
real z_star = 1048 * (1 + 0.00124*(Omega_b*h^2)^(-0.738)) * (1 + g1*((Omega_b + Omega_c)*h^2)^g2);
// get integral of u(z) ≡ 1/E(z) from 0 to z_star
real intofu = integrate_1d(u, 0, z_star, faketheta, x_r, x_i);
// wb
real wb = Omega_b*h^2;
// R
real R = sqrt(Omega_b + Omega_c) * intofu;
}
/* block for model definition */
model {
// priors
h ~ normal(0.7, 10);
Omega_b ~ normal(0.05, 10);
Omega_c ~ normal(0.27, 10);
Omega_r ~ normal(0, 10);
// likelihood
Rexp ~ normal(R, Rexp_err);
wbexp ~ normal(wb, wbexp_err);
}
/* block for generated quantites */
generated quantities {
}
The error
Compiling and running the previous model using CmdStanPy (with show_console=True, a single chain and sensible initial values) gives the following output:
15:28:43 - cmdstanpy - INFO - compiling stan file /home/joseferreira/tmp/model.stan to exe file /home/joseferreira/tmp/model
15:29:06 - cmdstanpy - INFO - compiled model executable: /home/joseferreira/tmp/model
15:29:06 - cmdstanpy - INFO - Chain [1] start processing
Chain [1] method = sample (Default)
Chain [1] sample
Chain [1] num_samples = 1000 (Default)
Chain [1] num_warmup = 1000 (Default)
Chain [1] save_warmup = 0 (Default)
Chain [1] thin = 1 (Default)
Chain [1] adapt
Chain [1] engaged = 1 (Default)
Chain [1] gamma = 0.050000000000000003 (Default)
Chain [1] delta = 0.80000000000000004 (Default)
Chain [1] kappa = 0.75 (Default)
Chain [1] t0 = 10 (Default)
Chain [1] init_buffer = 75 (Default)
Chain [1] term_buffer = 50 (Default)
Chain [1] window = 25 (Default)
Chain [1] algorithm = hmc (Default)
Chain [1] hmc
Chain [1] engine = nuts (Default)
Chain [1] nuts
Chain [1] max_depth = 10 (Default)
Chain [1] metric = diag_e (Default)
Chain [1] metric_file = (Default)
Chain [1] stepsize = 1 (Default)
Chain [1] stepsize_jitter = 0 (Default)
Chain [1] num_chains = 1 (Default)
Chain [1] id = 1 (Default)
Chain [1] data
Chain [1] file = /tmp/tmppcx6vxnf/2kwnz5li.json
Chain [1] init = /tmp/tmppcx6vxnf/zyatktdh.json
Chain [1] random
Chain [1] seed = 6437
Chain [1] output
Chain [1] file = /tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv
Chain [1] diagnostic_file = (Default)
Chain [1] refresh = 100 (Default)
Chain [1] sig_figs = -1 (Default)
Chain [1] profile_file = profile.csv (Default)
Chain [1] num_threads = 1 (Default)
Chain [1]
Chain [1]
15:29:06 - cmdstanpy - INFO - Chain [1] done processing
15:29:06 - cmdstanpy - ERROR - Chain [1] error: terminated by signal 11 Unknown error -11
Traceback (most recent call last):
File "/home/joseferreira/tmp/mwe-cmdstanpy.py", line 18, in <module>
fit = model.sample(data=data, show_console=True, chains=1, inits=inits)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/joseferreira/.micromamba/envs/cmdstanpy/lib/python3.11/site-packages/cmdstanpy/model.py", line 1188, in sample
raise RuntimeError(msg)
RuntimeError: Error during sampling:
Command and output files:
RunSet: chains=1, chain_ids=[1], num_processes=1
cmd (chain 1):
['/home/joseferreira/tmp/model', 'id=1', 'random', 'seed=6437', 'data', 'file=/tmp/tmppcx6vxnf/2kwnz5li.json', 'init=/tmp/tmppcx6vxnf/zyatktdh.json', 'output', 'file=/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv', 'method=sample', 'algorithm=hmc', 'adapt', 'engaged=1']
retcodes=[-11]
per-chain output files (showing chain 1 only):
csv_file:
/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv
console_msgs (if any):
/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906_0-stdout.txt
Consider re-running with show_console=True if the above output is unclear!
The important part is that Stan spits out an unknown error (terminated by signal 11 Unknown error -11) and exits, providing me no further information on what caused it.
Additionally, it speaks of two different files (/tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906.csv and /tmp/tmppcx6vxnf/modelobz90mtx/model-20221110152906_0-stdout.txt) which are nowhere to be found.
Although I suspect these two files will not bring any meaningful information
My main guess was that these initial values are in a situation where the denominator goes to zero, but that would trigger a warning, right?
My second guess was that integrating something which is being outputted by an ODE solver might give some sort of issues, but that shouldn’t be the case, right? Worst case scenario, it would simply take ages.
Besides the previous two rather wild guesses, I’m stuck.
Could this be due to Stan itself?
System
Everything took place in a new virtual environment, with packages installed from conda-forge
Python version: 3.11
CmdStanPy version: 1.0.8
CmdStan version: 2.30.1
OS: Ubuntu 18.04 (x64).
Compiler: gcc version 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04) - however, when I installed CmdStanPy, it also installed a newer version of gcc (10.4.0), which does not seem to be available in my PATH. This is the one which shows up when I write gcc -v, but it is the one which is being used by CmdStan?