Same model, different notation, different results

[Luca Campanelli]

Dear Stan users,
I modified a Stan model using matrix notation both for efficiency (as discussed in the manual) and flexibility (I can change the number of predictors without changing the Stan model), but I get different results and I'm not sure why. I tried the two types of notation on a simple regression example and it worked, but when used them on my actual model it doesn't. I’m using the same data, priors, and init values for testing.

I'm pasting below some parts of the two models: 1) works fine (meaningful results and no convergence issues or warnings; diagnostics fine), 2) it runs but nothing is good (meaningless results, bad diagnostics, take a lot of time).

Can you identify anything in model 2 that is not correct or can help explaining the different results?

Thank you for any suggestion you may have.
Luca

-----------------------------------------------------------------------------
1) NO MATRIX

// Model 1: NO MATRIX
data {
    int<lower=0> N;    // number of observations
    int<lower=0> J;    // number of subjects
    int subj[N];       // subjects
    vector[N] rt;      // predictor vector
    vector[N] dp;      // outcome vector
}
parameters {
    vector[J] incp;          // intercept
    vector[J] rate;          // rate
    vector[J] asym;          // asymptote

    real<lower=0> sigma;       // L1 error
   
    real mu_ua;
    real mu_ur;
    real mu_ui;

    real<lower=0> sigma_ua;    // L2 error asy
    real<lower=0> sigma_ur;    // L2 error rate
    real<lower=0> sigma_ui;    // L2 error int
}
transformed parameters {
    vector[N] mu;
   
    vector[N] ast;
    vector[N] rat;
    vector[N] itt;

    for (n in 1:N) {
        ast[n] = asym[subj[n]];
        rat[n] = rate[subj[n]];
        itt[n] = incp[subj[n]];
       
        mu[n] = fmax(0, ast[n] * (1.0 - exp( -rat[n] * (rt[n] - itt[n] ) )) );
    }
}
model {
    asym ~ normal(mu_ua, sigma_ua);
    incp ~ normal(mu_ui, sigma_ui);
    rate ~ normal(mu_ur, sigma_ur);

    mu_ua ~ normal(0,100);
    mu_ui ~ normal(0,100);
    mu_ur ~ normal(0,100);
   
    sigma ~ cauchy(0,5);  
    sigma_ua ~ cauchy(0,5);
    sigma_ui ~ cauchy(0,5);
    sigma_ur ~ cauchy(0,5);

    dp ~ normal(mu, sigma);
}

-----------------------------------------------------------------------------
2) MATRIX

// Model 2: MATRIX
data {
    int<lower=0> N;    // number of observations
    int<lower=0> J;    // number of subjects
    int subj[N];       // subjects indicator
    vector[N] rt;      // predictor vector
    vector[N] dp;      // outcome vector
   
    // these matrices have only one column with the value of 1 (for intercept)
    matrix[N,1] int_x;         // interc: fixed effects model
    matrix[N,1] rat_x;         // rate: fixed effects model
    matrix[N,1] asy_x;         // asym: fixed effects model
}
parameters {
    matrix[J,1] inte;       // interc: fixed effects coefficients
    matrix[J,1] rate;       // rate: fixed effects coefficients
    matrix[J,1] asym;       // asym: fixed effects coefficients
   
    real<lower=0> sigma;      // L1 error
   
    real<lower=0> sigma_i;      // subj inte ranef
    real<lower=0> sigma_r;      // subj rate ranef
    real<lower=0> sigma_a;      // subj asym ranef

    real mu_i;
    real mu_r;
    real mu_a;    
}
transformed parameters{
    vector[N] mu;
   
    vector[N] ast;
    vector[N] rat;
    vector[N] itt;
   
    itt = rows_dot_product(int_x, inte[subj] );
    rat = rows_dot_product(rat_x, rate[subj] );
    ast = rows_dot_product(asy_x, asym[subj] );
   
    for (n in 1:N) {
        mu[n] = fmax(0, ast[n] * (1.0 - exp( -rat[n] * (rt[n] - itt[n] ) )) );
    }
}
model {  
    sigma ~ cauchy(0,5);
    sigma_a ~ cauchy(0,5);
    sigma_i ~ cauchy(0,5);
    sigma_r ~ cauchy(0,5);
   
    mu_i ~ normal(0,100);
    mu_r ~ normal(0,100);
    mu_a ~ normal(0,100);
   
    inte[,1] ~ normal(mu_i, sigma_i);
    rate[,1] ~ normal(mu_r, sigma_r);
    asym[,1] ~ normal(mu_a, sigma_a);

    dp ~ normal(mu, sigma);
}

[Bob Carpenter]

I'd suggest testing that you're getting the same results with
the matrix versions. You can either print out the results and
compare by eye, or you can write both in one program and test
with Stan code, or you can write a Stan function and test through R.
They can't be calculating the same function if you get different
results. Stan isn't sensitive to how things are computed other than
in terms of speed.

fmax is really bad news as it breaks derivatives.

You want to use non-centered parameterizations for asym,
rate, etc. if you don't have a lot of data.

You can replace this:

>
> vector[N] ast;
> vector[N] rat;
> vector[N] itt;
>
> for (n in 1:N) {
> ast[n] = asym[subj[n]];
> rat[n] = rate[subj[n]];
> itt[n] = incp[subj[n]];

with

vector[N] ast = asym[subj];

and so on.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Luca Campanelli]

Thank you for your answer, Bob.
I finally found the problem. It was that my init values were correctly passed in model 1, but not in model 2 - and given the nature of the model without good init values things don’t work.
I was passing vectors while the parameters were matrices, and given the absence of warnings or errors messages I thought things were working. Then I looked at get_inits() and saw the problem.

For fmax, I could code things differently, for example:

      mu[n] = my_equation ;
      if (mu[n] > 0)
          dp[n] ~ normal(mu[n], sigma);
      else
          dp[n] ~ normal(0, sigma);

or

    mu[n] = ast[n] * (1.0 - exp( -rat[n] * (rt[n] - itt[n] ) )) ;
    mu[n] = mu[n] > 0 ? mu[n] : 0 ;
    dp[n] ~ normal(mu[n], sigma);

but I don’t think that would address the issue you mentioned. Is there any better way to do it?

Re your suggestion to recode asym etc, that would be the same model, we are just moving them outside of the loop. It’s only an efficiency issue, correct?

Thank you.
Luca

[Jonah Gabry]

On Friday, March 3, 2017 at 4:03:50 PM UTC-5, Luca Campanelli wrote:

Thank you for your answer, Bob.
I finally found the problem. It was that my init values were correctly passed in model 1, but not in model 2 - and given the nature of the model without good init values things don’t work.
I was passing vectors while the parameters were matrices, and given the absence of warnings or errors messages I thought things were working. Then I looked at get_inits() and saw the problem.

Hi Luca, how were you specifying the initial values? It's possible you specified them in a way that is valid but that didn't correspond to what you were intending, in which case that wouldn't result in any warnings/errors. Or maybe it wasn't a valid specification and we need to fix something so an error or warning is thrown. If rstan is silently ignoring incorrectly specified inits then that's a bug. 

 

Re your suggestion to recode asym etc, that would be the same model, we are just moving them outside of the loop. It’s only an efficiency issue, correct? 

Yeah that doesn't change the statistical model at all. What Bob's version gives you is cleaner code. It's probably more efficient too, but the loop you have isn't super inefficient. In general loops are pretty fast in Stan (it's ultimately C++). However, where you really do lose a lot efficiency by looping is when calls to the density functions are involved. When those can be vectorized you can get a big boost in efficiency.

Jonah

[Luca Campanelli]

Hi Jonah,
for this parameter:
    parameters {
        matrix[J,n_cond_int] inte;

this is the init that didn’t work and there were no warnings:
    "inte[,1]" = c(rep(2, sdata$J))

and this is what worked:
    inte = matrix(c(rep(2, sdata$J)), nrow = sdata$J, ncol = 1, byrow = FALSE)

It seems that RStan correctly returns error messages when the name of the parameter in the init list is correct but the values are not (number of values or data type are incorrect).

However, there are no error messages or warnings if the name of the parameter in the init list does not match any parameter in the Stan model; it could be something that one has tried (like in my case) or a typo and the init values are not passed and there are no warnings.
I think that in those cases it could be useful to have warnings saying that there is no matching parameter in the model for init “initname”.

Thank you,
Luca

[Bob Carpenter]

I think a warning is a good idea. Right now, the I/O is handled on a per-interface
basis, so it'd have to be separate code for RStan, PyStan, and CmdStan.

We don't want to throw an error because a common use case is to set up all
the data and then build models that use various subsets of predictors, etc.

- Bob