Spatial Model Compiling Errors

[Colin]

Dear Stan Group,

I am trying to fit a GP spatial model over T time points (years) where there is spatial dependence within each year, but across years independence is assumed and the spatial parameters do not vary.  Thus, the covariance matrix is block diagonal (D \otimes Identity matrix).  The spatial structure and parameters are the same across years, so I am trying to just have to take the inverse of the covariance matrix once for all T years.  

I have fit the model for individual years in Stan, but when I try and estimate the parameters over all T, I'm getting long source code/compiler errors. I initially tried to separate the T time points using an array for the model matrix, but ran into trouble so I created an ID variable instead for each T.  

Below is my Rstan code, and I can attach the data file if that would be helpful.  The data set is not terribly large (180 observations).  I am wondering if the way I am separating the years is causing the issue, but I am not sure.  Any insight would be much appreciated. 

data{

  int<lower=1> N; //total obs (L*T)

  int<lower=1> L; // total number of locations per year

  int<lower=1> T; //number of years, assume independence across years

  int<lower=1> p; //number of covariates

  int <lower=1> id[N];  // ID By T (year)

  vector [N] Y; // response vector

  matrix[N,p] X;// Model Matrix

  matrix [L,L] d; //distance matrix, same across all years so only need to define once

}

parameters {

  vector[p] beta; # coefficient

  real<lower=0> sigma; # error sd

  real<lower=0> tau; # sd spatial process

  real<lower=0> phi; # range/decay parameter

}

transformed parameters {

  real<lower=0> sigma2;

  real<lower=0> tau2;

  sigma2 = sigma^2;

  tau2 = tau^2;

  

}

model {

  vector[N] mu; 

  matrix[L,L] H; // Same Spatial Parameters for Each Year

  matrix[L,L] Hinv;

  

  for (i in 1:(L-1)) {

    for (j in (i+1):L) {

        H[i,j] = tau2 * exp(-d[i,j]/ phi);

        H[j,i] = H[i,j]; // do not calculate twice, utilize symmetry

      }

    }

  for (k in 1:L) {

    H[k, k] = sigma2 + tau2; //diagonal of covariance matrix

  }

  

  Hinv = inverse_spd(H); //inverse of covariance matrix, should just have to do this once 

  

  

  target += -.5*log_determinant(H); // this step only has to be done once across all years

  

  //calculate the mean for each year using index variable

  for (n in 1:N){

  mu[n] = Y[id[n]] - X[id[n]]*beta;

  }

  

  // increment the likelihood for each T (year)

  for (t in 1:T){

  target += -.5*(mu[id[t]])'*Hinv*(mu[id[t]]);  //MVN is block diagonal so likelihood is incremented by year

  }

  //priors

  //inproper prios on betas.

  sigma ~ cauchy(0,1);

  tau ~ cauchy(0,1);

  phi ~ gamma(0.5,0.5);

}

[Ben B]

Real quick, is there a reason you're incrementing the log probability manually like this and not using statements like:

y ~ multi_normal(mu, H);

Cause if you can do stuff like that, then you can just do:

L = cholesky_decompose(H);

Followed by a bunch of (with the ys and mus changed to whatever they need to be):

y ~ multi_normal_cholesky(mu, L);

As described in the manual (section 17.4). Computing the Cholesky just once per step should be the speedup you're looking for. Maybe if you can write your model like this it'll be easier to figure out what's going wrong in the model?

Ben

[Colin Lewis-Beck]

Hi Ben,

Thanks for the response.  It seemed more difficult to make a big block diagonal covariance matrix for all observations, and since I am assuming each time point (T) is independent, factoring the likelihood into individual contributions for each year would be a computationally better choice as you are inverting a much smaller dimensional matrix.

Colin

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben B]

Howdy Colin

Yeah I agree with you. I meant that the 'H' variable would be one of these small blocks.

So if I understand your model correctly, each set of data you model as a Gaussian process with different means but the same covariance.

So:

for (d in 1:D)
    y[d] ~ multi_normal(mu[d], H);

Where 'd' are your days, and H is the small covariance matrix (the same H as in your model, not the covariance matrix for all D days lumped together).

Can't you just write your model like this? There is the issue that at each day you need to have data measured at the same points for them to share the covariance matrix, but if you have that, then you can get the speedup you're looking for by doing:

L = cholesky_decomposition(H)
for (d in 1:D)
    y[d] ~ multi_normal_cholesky(mu[d], L);

It's a Cholesky instead of an inversion cause computing matrix inverses by themselves is a Bad Idea. Solving Ax = b is strongly preferred over trying to compute A^-1 and then multiplying it by b. The Cholesky factorization makes it way easier to do the solves, and so will have the effect you're looking for.

Ben

[Colin Lewis-Beck]

Hi Ben,

Thanks so much for your reply!  It is very helpful and exactly what I am trying to do.  I have one more quick question about how to subset the response vector y and the model matrix X, by the time variable.  Right now y is a vector and X is a matrix, each of length N, the full data set.  It seems like to be able to loop over time, say T (or as you wrote, D), you need to index the vector y and matrix X into separate groups.  I guess the alternative would be to pass them individually in the data step, but that seems more complicated.  I've seen indexing done over all N in the Stan manual and other examples, but not when it comes to looping over a subgroup.

Below is the set up I am currently working with per your suggestion.  From looking at the Stan manual you can subset vectors via indexing y[1:30], but since I do have an id vector in the data block, it seems like there should be an easier way.  Rather than subsetting X, I could also calculate mu=X*beta and then try and subset mu based on the id variable?  Sorry if these are silly questions, but I'm still trying to get familiar with Stan's indexing.  Any help you could provide would be much appreciated.

Thanks again,

Colin

data{

  int<lower=1> N; //total obs (L*T)

  int<lower=1> L; // total number of locations per year

  int<lower=1> T; //number of years, assume independence across years

  int<lower=1> [N] id; // vector of length N to id Year

  int<lower=1> p; //number of covariates

  vector [N] Y; // response vector

  matrix[N,p] X;// Model Matrix

  matrix [L,L] d; //distance matrix, same across all years

}

parameters {

  vector[p] beta; # coefficient

  real<lower=0> sigma; # error sd

  real<lower=0> tau; # sd spatial process

  real<lower=0> phi; # range/decay parameter

}

transformed parameters {

  real<lower=0> sigma2;

  real<lower=0> tau2;

  sigma2 = sigma^2;

  tau2 = tau^2;

  

}

model {

  vector[T] mu[L];  //save T vectors of means

  matrix[L,L] H; // Same Spatial Paramters for Each Year

  

  for (i in 1:(L-1)) {

    for (j in (i+1):L) {

      H[i,j] = tau2 * exp(-d[i,j]/ phi);

      H[j,i] = H[i,j]; // do not calculate twice, utilize symmetry

    }

  }

  for (k in 1:L) {

    H[k, k] = sigma2 + tau2; //diagonal of covariance matrix

  }

  

  L = cholesky_decomposition(H)

  

  for (t in 1:T){      

    mu[t] = X[t]*beta  //how to break up model Matrix X into Years?

  }

  

  for (t in 1:T)

    y[t] ~ multi_normal_cholesky(mu[t], L);  

  

  //priors

  sigma ~ cauchy(0,1);

  tau ~ cauchy(0,1);

  phi ~ gamma(0.5,0.5);

}

--

[Ben B]

The way I'd do it is arrange y like:

vector[L] Y[T];

So that means you have T vectors of length L (like you have arranged mu in your model block). Then you can just use:

for (t in 1:T)
    y[t] ~ multi_normal_cholesky(mu[t], L);

Ben

[Bob Carpenter]

You can also just vectorize multi_normal if y is an array
of vectors.

- 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.

[Bob Carpenter]

I should've responded to this one. If you declare

vector[L] y[T];

then you can just write

y ~ multi_normal(mu, Sigma);

rather than Cholesky decomposing and the multi_normal will do
it for you and everything gets vectorized internally.

- 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.

[Colin Lewis-Beck]

Thanks for the follow up Bob.  That will definitely speed up the code.  

One more quick question about subsetting the y and mu vector.   As you can see below, I am sub-setting them by indexing manually.  I know Rstan has a segment function to subset vectors.  However, is there a way to append or add an id variable for each year so I can separate y and mu by year using a loop?  

Thanks again,

Colin

 data{

  int<lower=1> N; //total obs (L*T)

  int<lower=1> L; // total number of locations per year

  int<lower=1> T; //number of years, assume independence across years

  int id [N]; // vector of length N to id Year

  int<lower=1> p; //number of covariates

  vector [N] Y; // response vector

  matrix[N,p] X;// Model Matrix

  matrix [L,L] d; //distance matrix, same across all years

}

transformed data{

  vector [L] Ysub [T];  //break up Y by years T with L obs each Year

  Ysub[1] = Y[1:30];

  Ysub[2] = Y[31:60];

  Ysub[3] = Y[61:90];

  Ysub[4] = Y[91:120];

  Ysub[5] = Y[121:150];

  Ysub[6] = Y[151:180];

}

parameters {

  vector[p] beta; # coefficient

  real<lower=0> sigma; # error sd

  real<lower=0> tau; # sd spatial process

  real<lower=0> phi; # range/decay parameter

}

transformed parameters {

  real<lower=0> sigma2;

  real<lower=0> tau2;

  sigma2 = sigma^2;

  tau2 = tau^2;

  

}

model {

   vector [N] mu  

   vector [L] mu_yr[T];  //for Group Means

   matrix[L,L] H; // Same Spatial Paramters for Each Year

   

   for (i in 1:(L-1)) {

     for (j in (i+1):L) {

       H[i,j] = tau2 * exp(-d[i,j]/ phi);

       H[j,i] = H[i,j]; // do not calculate twice, utilize symmetry

     }

   }

   for (k in 1:L) {

     H[k, k] = sigma2 + tau2; //diagonal of covariance matrix

   }

   

  mu = X*beta; //overall mean vector

   

   //seperate by year

  mu_yr[1] = mu[1:30];

  mu_yr[2] = mu[31:60];

  mu_yr[3] = mu[61:90];

  mu_yr[4] = mu[91:120];

  mu_yr[5] = mu[121:150];

  mu_yr[6] = mu[151:180];

   

   

  Ysub ~ multi_normal(mu_yr, Sigma);

   

   //priors

   sigma ~ cauchy(0,1);

   tau ~ cauchy(0,1);

   phi ~ gamma(0.5,0.5);

}

On Thu, May 4, 2017 at 8:48 AM, Bob Carpenter <ca...@alias-i.com> wrote:

I should've responded to this one.  If you declare

  vector[L] y[T];

then you can just write

  y ~ multi_normal(mu, Sigma);

rather than Cholesky decomposing and the multi_normal will do
it for you and everything gets vectorized internally.

- Bob

> On May 3, 2017, at 11:52 AM, Ben B <bbba...@gmail.com> wrote:
>
> The way I'd do it is arrange y like:
>
> vector[L] Y[T];
>
> So that means you have T vectors of length L (like you have arranged mu in your model block). Then you can just use:
>
> for (t in 1:T)
>     y[t] ~ multi_normal_cholesky(mu[t], L);
>
> Ben
>
> --
> 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+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

> On May 4, 2017, at 11:14 AM, Colin Lewis-Beck <clewi...@gmail.com> wrote:
>
> Thanks for the follow up Bob. That will definitely speed up the code.
>
> One more quick question about subsetting the y and mu vector. As you can see below, I am sub-setting them by indexing manually. I know Rstan has a segment function to subset vectors. However, is there a way to append or add an id variable for each year so I can separate y and mu by year using a loop?

I'm not sure what you're asking in the way of modeling. Obviously
you can add an ID variable,

int ID[Y];

but what do you want to do with it? Do you want a varying intercept
or something?

- Bob

> > 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.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --

> 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.

[Colin Lewis-Beck]

Hi Bob,

Thanks.  Yes, I would like a varying slope and intercept for each time point, but with a constant covariance matrix.  I think I figured out the way to subset the data correctly and model seems to be running fine now.  See below if anyone else is interested.  

data{

  int<lower=1> L; // total number of locations per year

  int<lower=1> T; //number of years, assume independence across years

  matrix [L,T] Y; // response vector

  matrix [L,T] X;

  matrix [L,L] d; //distance matrix, same across all years

}

parameters {

  real alpha [T]; //intercept

  real beta [T]; //slope

  real<lower=0> sigma; # error sd

  real<lower=0> tau; # sd spatial process

  real<lower=0> phi; # range/decay parameter

}

transformed parameters {

  real<lower=0> sigma2;

  real<lower=0> tau2;

  sigma2 = sigma^2;

  tau2 = tau^2;

  

}

model {

   matrix[L,L] H; // Same Spatial Paramters for Each Year

   

   for (i in 1:(L-1)) {

     for (j in (i+1):L) {

       H[i,j] = tau2 * exp(-d[i,j]/ phi);

       H[j,i] = H[i,j]; // do not calculate twice, utilize symmetry

     }

   }

   for (k in 1:L) {

     H[k, k] = sigma2 + tau2; //diagonal of covariance matrix

   }

  

   for (i in 1:T){

   Y[,i] ~ multi_normal(alpha[i] + X[,i]*beta[i], H);

}

   

   //priors 

   //inproper priors on coefficients

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> 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+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

Thanks much. We really appreciate it when people report back. The new discourse
interface lets us mark issues as closed/answered, which should be helpful.

- Bob

> > > 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.
> >
> > --
> > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.

> > To unsubscribe from this group and all its topics, 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.
> >
> >
> > --
> > 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.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/iEbPDMqPZN4/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --
> 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.