Handling multiple (unconstrained) parameter sets in `log_prob`

[tristansnowsill]

(First to briefly say, I love Stan and CmdStan!)

I think there is a row-major / column-major mismatch error in the log_prob command of CmdStan when you want to provide more than one set of (unconstrained) parameters.

I believe it can be fixed by swapping lines 562 and 563 of command_helper.hpp to change the order of filling the results:

cmdstan/src/cmdstan/command_helper.hpp

Lines 562 to 567 in 2ba7253

	for (size_t i = 0; i < u_params_rows; ++i) {
	for (size_t j = 0; j < u_params_cols; ++j) {
	params_r_ind[i][j] = *(u_params_r.data() + idx);
	++idx;
	}
	}

By which I mean, when I swapped them and recompiled CmdStan it all seemed to work as expected for me.

To demonstrate, take a simple model with multiple parameters:

data {
  int<lower=0> J; // number of schools
  array[J] real y; // estimated treatment effect (school j)
  array[J] real sigma; // std err of effect estimate (school j)
}
parameters {
  real mu;
  real<lower=0> tau;
  array[J] real<offset=mu,multiplier=tau> theta;
}
model {
  theta ~ normal(mu, tau);
  y ~ normal(theta, sigma);
}

A data file

{
  "J": 8,
  "y": [28, 8, -3, 7, -1, 1, 18, 12],
  "sigma": [15, 10, 16, 11, 9, 11, 10, 18]
}

And with two different files for unconstrained parameters:

{
  "params_r": [
    0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10
  ]
}

{
  "params_r": [
    [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10],
    [1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.10],
    [2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 2.09, 2.10]
  ]
}

When I run the second parameters file (3 parameter sets), it sets

Parameter	Run 1	Run 2	Run 3
mu	0.01	1.04	2.07
tau	1.01	2.04	0.08
theta[1]	2.01	0.05	1.08
theta[2]	0.02	1.05	2.08
theta[3]	1.02	2.05	0.09
theta[4]	2.02	0.06	1.09
theta[5]	0.03	1.06	2.09
theta[6]	1.03	2.06	0.10
theta[7]	2.03	0.07	1.10
theta[8]	0.04	1.07	2.10

[WardBrian]

Thanks, this is a nice catch and something we should have been testing.

If you're not already aware, you may wish to use https://roualdes.us/bridgestan/latest/ as a better (in the sense that it is both easier to use and faster) way to query the log probability of a model.

[tristansnowsill]

Thanks, this is a nice catch and something we should have been testing.

I will confess I still don't understand why the original code didn't work... But I am not super familiar with modern C++.

If you're not already aware, you may wish to use https://roualdes.us/bridgestan/latest/ as a better (in the sense that it is both easier to use and faster) way to query the log probability of a model.

Oh cool I didn't know about this! I am probably going to be using MPI so might not be relevant this time, but that is super helpful to know about.