Problem handing over data to fixed_param run with cmdstanpy

Hi all,

I want to infer parameters of a system of ODEs. Before working with real data, I want to generate data for many synthetic experiments.
In order to familiarize myself with stan, I’m trying different interfaces and as a first step, I want to solve the ODE from the user guide.

This is the sho.stan file I’m using:

functions {
  vector sho(real t,        // time
             vector y,      // state
             real theta) {  // friction parameter
    vector[2] dydt;
    dydt[1] = y[2];
    dydt[2] = -y[1] - theta * y[2];
    return dydt;
  }
}

transformed data {
  int<lower=1> T = 40;
  vector[2] y0 = [1.0, 0.0]';
  real t0 = 0.0;
  array[T] real ts;
  for (t in 1:T) {
    ts[t] = t;
  }
  real theta = 0.35;
}

model {
}

generated quantities {
  array[T] vector[2] y_sim = ode_rk45(sho, y0, t0, ts, theta);

  // add measurement error
  for (t in 1:T) {
    y_sim[t, 1] += normal_rng(0, 0.1);
    y_sim[t, 2] += normal_rng(0, 0.1);
  }
}

I get expected results, when interfacing that code with cmdstanpy:

from pathlib import Path
from cmdstanpy import CmdStanModel
import matplotlib.pyplot as pt


stan_file = Path('stan/sho.stan')

model = CmdStanModel(stan_file=stan_file)

sim_data = model.sample(data={}, fixed_param=True)

y_sim = sim_data.stan_variable(var='y_sim')
y_sim_mean = np.mean(y_sims, axis=0)

pt.plot(y_sim_mean[:,0], y_sims_mean[:,1])
pt.plot(y_sim[0, :,0], y_sims[0, :,1])
pt.show()

However, I want to get the data from Python. If I set up a dictionary with

data = {
    'T': 40,
    'y0': [1.0, 0.0],
    't0': 0.0,
    'ts': list(range(1, 41)),
    'theta': 0.35,
}

And exchange the tranformed data in the stan code with a data block

data {
  int<lower=1> T;
  vector[2] y0;
  real t0;
  array[T] real ts;
  real theta;
}

I get following error Exception: variable does not exist; processing stage=data initialization; variable name=T; base type=int.

If I try the same with pystan, it works:

from pathlib import Path
import stan
import matplotlib.pyplot as pt


stan_code = Path('stan/sho.stan').read_text()

data = {
    'T': 40,
    'y0': [1.0, 0.0],
    't0': 0.0,
    'ts': list(range(1, 41)),
    'theta': 0.35,
}

posterior = stan.build(stan_code, data=data, random_seed=1)

fit = posterior.fixed_param()

y_sim = fit['y_sim']
y_sim_mean = np.mean(y_sim, axis=2)

pt.plot(y_sim_mean[:,0], y_sim_mean[:,1])
pt.plot(y_sim[:, 0, 0], y_sim[:, 1, 0])
pt.show()

Could someone give me a hint to how I get this working with cmdstanpy? Thank you so much in advance!

Operating System: Linux
Interface Version: pystan 3.4.0, cmdstanpy 1.0.1
Compiler/Toolkit: stan 2.29

Did you recompile your model after the change?

Ah, good point! I set `CmdStanModel(stan_file=stan_file, compile=‘force’) and I arrived at a point, where I think I was alrady yesterday, but I couldn’t remember how I got there.

This is now the result of print(sim_data.summary()):

             Mean  MCSE  StdDev   5%  50%  95%  N_Eff  N_Eff/s  R_hat
name                                                                 
lp__          0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[1,1]    0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[1,2]    0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[2,1]    0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[2,2]    0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
...           ...   ...     ...  ...  ...  ...    ...      ...    ...
y_sim[38,2]   0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[39,1]   0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[39,2]   0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[40,1]   0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN
y_sim[40,2]   0.0   NaN     0.0  0.0  0.0  0.0    NaN      NaN    NaN

[81 rows x 9 columns]

So now you are able to call sample with data?

How are the “raw” samples from that fit?

Oh dear, my ts variable was actually one index too short in the cmdstanpy example.
Now, with the compile='force' and the correct initialisation of ts, everything works just fine.
Thank you so much for the really fast responses!

1 Like