Problem handing over data to fixed_param run with cmdstanpy

Interfaces
1

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

2

Did you recompile your model after the change?

3

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]
4

So now you are able to call sample with data?

How are the “raw” samples from that fit?

5

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