Parameters declared in the model block

Hi,

I am building a model Y=X\beta + \epsilon, where both X and Y are matrices. I want to model the correlation between the Y.

As you can see, I have some pretty big parameter matrices, so I want to declare them in the model block so that they can be discarded. But I got error, saying “std_normal_lpdf: Random variable[1] is nan, but must not be nan!”

However, if I move the declaration to the transformed parameters blocks, then it runs fine, although very slow.

Can you advise where it went wrong?

Also beside the error message, I feel that there are many places that the code can be improved, can you also advise on that?

Thanks a lot!

data {

  int<lower=0> nSamples; //number of samples 
  int<lower=0> nBio; //number of y
  int<lower=0> nX; // number of x

  matrix[nSamples, nBio] y; // the multivariate outcome matrix
  matrix[nSamples, nX] X; // predictor matrix
  
}

parameters {

  matrix[nX, nBio] beta; // betas from N(0,1)
 
  vector<lower=0>[nBio] sigma_eps; // sd of the residual
  
  vector<lower=0>[nBio] tau;     // prior scale
  
  cholesky_factor_corr[nBio] L;

}

transformed parameters {

}

model {
  matrix[nSamples, nBio] z;
  matrix[nSamples, nBio] mu;
  
  mu = X * beta + z * (diag_pre_multiply(tau,L))';

  sigma_eps ~ exponential(1);
 
  tau ~ cauchy(0, 2.5);

  L ~ lkj_corr_cholesky(1);
 
  for(i in 1:nX){
    to_vector(beta[i]) ~ student_t(10, 0, 10);
  }

  for(i in 1:nSamples){
    to_vector(z[i]) ~ std_normal();
  }
  
  for(i in 1:nSamples){
    for(j in 1:(nBio)){
      y[i][j] ~ normal(mu[i][j], sigma_eps[j]);
    }
  }
}

z seems to be a true parameter proper, so you can’t omit it from the parameters section.

Also, this:

  for(i in 1:nX){
    to_vector(beta[i]) ~ student_t(10, 0, 10);
  }

is way more efficient as:

to_vector(beta) ~ student_t(10, 0, 10);

Ditto:

  for(i in 1:nSamples){
    to_vector(z[i]) ~ std_normal();
  }

should be:

    to_vector(z) ~ std_normal();

Cool! Thanks so much!

Oh, I also think something is fundamentally wrong with your model; it doesn’t make sense that mu and z have the same dimensions.

Oh, nevermind, that’s fine. I’d been misreading what mu was.

Yup, since we don’t care about z in here. Anyway to throw it away after each sampling iteration? carrying them along eats up too much memory.

Unfortunately no, though they aren’t kept in memory but instead written to file. Are you observing gradually increasing ram during sampling?

I see. I didn’t specifically test if Z increase my memory, but keeping mu definitely increased my fitted object size. And they significantly slowed down the program.

Btw, see the “hwg” model at the repo here for a demo of a possibly-faster version of what you’re doing.

I worked on some other project with many Ys. In that case, keeping all the intermediate values gave me hundred of Gb of fitted object. Getting rid of those slim the object to a couple of Mb

will do! Thanks a bunch!

1 Like