How would one declare a covariance matrix with different bounds on different entries of the matrix?

[Qian Zhang]

J is an integer read in as data.

I'd like to use a J x J covariance matrix Sigma in the model block that has diagonal entries bounded between [0,1], but off-diagonal entries bounded in terms of diagonal entries like so:

0 < Sigma[i,j] < sqrt(Sigma[i,i] * Sigma[j,j])

How would I pull that off in STAN?

This is what I attempted. The difficulty with this approach is described in the comment between /* and */.

parameters{

       // Declare diagonal entries of cs matrix (declared later)

       vector<lower=0, upper=1>[J] c1s;

   

       /* Declare off-diagonal entries of cs matrix (declared later).

       

         CONTINUE HERE: This is wrong because it defines a vector c12s for eacher (j, jprime),

         whereas I want a single c12s vector with each entry bounded differently,

         according to the entries of c1s. I don't know how to do that in STAN.

       

         If I knew J, e.g. J = 3, then instead of having a vector c12s,

         I could have a number c12s_1, c12s_2, ..., c12_3 for each upper diagonal entry 

         of the cs matrix, and have 1 statement per c12s_i,

         e.g. real<lower = 0, upper = sqrt(c1s[1]*c1s[2])>c12s_1,

         but I don't know J. J is variable. */

 

       for (j in 1:J){

              for (jprime in (j+1):J){

                     real<lower=0, upper=sqrt(c1s[j]*c1s[jprime])>c12s[(J^2 - J)/2];

              }

       }

}

transformed parameters {

        // Fill in diagonal entries with the parameters defined with constraints earlier

        for(j in 1:J){

                cs[j,j] <- c1s[j];

        }

  

        // Declare off-diagonal entries

        for(j in 1:J){

                for(jprime in (j+1):J){

                        cs[j,jprime] <- c12s[(jprime - j) + (j-1)*(J-1)];

                        cs[jprime,j] <- cs[j,jprime];

                }

        }

}

I saw something similar here: https://groups.google.com/forum/#!topic/stan-users/ibEPXA7apAU, for bounding the first entry of a vector, leaving the other entries unbounded, but the bound in that example was done on a single entry, not a vector, like so:

parameters { 
   real<lower=0,upper=1> x1; 
   real x2_N[N-1]; 
} 
transformed parameters { 
   real x[N]; 
   x[1] <- x1; 
   for (n in 1:(N-1)) 
     x[n+1] <- x2_N[n]; 
} 

[Ben Goodrich]

On Tuesday, January 26, 2016 at 11:32:38 PM UTC-5, Qian Zhang wrote:

J is an integer read in as data.

I'd like to use a J x J covariance matrix Sigma in the model block that has diagonal entries bounded between [0,1], but off-diagonal entries bounded in terms of diagonal entries like so:

0 < Sigma[i,j] < sqrt(Sigma[i,i] * Sigma[j,j])

How would I pull that off in STAN?

I would do it by declaring a Cholesky factor of a correlation matrix

cholesky_factor_corr[J] L_corr;

and a vector of standard deviations that are bounded

vector<lower=0,upper=1>[J] sigma;

as parameters, and then

diag_pre_multiply(sigma, L_corr);

is a Cholesky factor of a covariance matrix whose diagonal elements would be between 0 and 1.

Ben

[Qian Zhang]

Thanks Ben!

So it turns out I have informative priors on entries of the Sigma matrix. The problem with using these sigma and Cholesky matrices you mention is I don't know how to derive priors on the sigma and Cholesky matrices from the priors on entries of Sigma, and it seems not-worth-it to go down that route.

I tried 2 approaches:

(1) Declare Sigma as a cov_matrix in the parameters block, put priors on its entries in the model block, and rely on rejection of proposed Sigma values that are out of bounds (Code Excerpt 1 below)

(2) Define Sigma as a cov_matrix in the transformed parameters block, define it in terms of sigma and Cholesky matrices, put priors on Sigma entries but not on sigma or Cholesky matrices (Code Excerpt 2 below)

Both approaches gave me draws. Approach 1 was a lot faster (~ a few minutes) than Approach 2 (~ 30 min), so I'm going with Approach 1. (I guess you put bounds in Approach 2 to encourage movement to a new position in the parameter space at the next draw, while in Approach 1 you tend to stay put at the last parameter draw because you propose something out of bounds and then reject it, so the only difference is mixing time, both real human time and time in terms of number of draws.)

Question 0: Both approaches work, i.e. I'm not coding up something totally non-sensical either way? They should both implement MCMC procedures to draw from the same posterior.

Question 1: Since Sigma is a cov_matrix, in the model block I only need to put priors on the upper diagonal entries of Sigma, and the lower diagonal entries will automatically be filled in with the upper diagonal draws?

I got many warnings in Code Excerpt 1 of the following 3 types.

Type 1:

                                                                                                                                                                count

validate transformed params: thetaB is -0.00107406, but must be greater than or equal to 0                 4

 

Type 2:

Exception thrown at line 97: beta_log: Random variable is 1.00308, but must be less than or equal to  2

Type 3:

Warning messages:

1: There were 132367 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. 

2: Examine the pairs() plot to diagnose sampling problems

Question 3: Is the following understanding correct? My understanding is Types 1 and 2 come from drawing Sigma entries that are out of bounds, or result in out-of-bounds thetaB, in which case such draws are rejected. Also, I never need to worry about this divergent transitions warning message, as you say here: “There are almost always some divergent transitions during the warmup phase but we don't really worry about them. It is trying to find a stepsize during the warmup phase that achieves the target acceptance probability without diverging and without hitting the max_treedepth.” (https://groups.google.com/forum/#!topic/stan-users/hjYWqJZfSgw) But then why do you have this warning message at all?

Code Excerpt 1:

parameters {

  //  Covariance matrix. I don't define constraints on it and rely on rejection in the model block.

  cov_matrix[J] Sigma;

  

}

transformed parameters {

  

  real<lower=0,upper=1> thetaB;

  

  thetaB <- 0;

  

  for(j in 1:J){

        for(jprime in 1:J){

            if(j != jprime){

                 thetaB <- thetaB + Sigma[j,jprime];

            }

        }

  }

  

  thetaB <- thetaB / (J^2 - J);

}

model{

  

  // Try putting priors on the entries of Sigma

  

  for(j in 1:J){

      Sigma[j,j] ~ beta(1,3);

  }

  

  for(j in 1:J){

      for(jprime in (j+1):J){

          

            // CONTINUE HERE: Hope that I only need to specify one off-diagonal

            // and that the other will get filled in by the fact that Sigma is a covariance matrix

            Sigma[j, jprime] ~ uniform(0.0, sqrt(Sigma[j,j] * Sigma[jprime,jprime]));    

      }    

  }

}

Code Excerpt 2: 

parameters {

  // Cholesky factor of a correlation matrix

  cholesky_factor_corr[J] L_corr;

  // Vector of standard deviations that are bounded

  vector<lower=0,upper=1>[J] sigma;

}

transformed parameters {

  

    cov_matrix[J] Sigma;

   # p.74 STAN manual shows how to multiply matrices

   Sigma <- diag_pre_multiply(sigma, L_corr) *  diag_pre_multiply(sigma, L_corr)';

  

}

model{

  

  // Try putting priors on the entries of Sigma

  

  for(j in 1:J){

      Sigma[j,j] ~ beta(1,3);

  }

  

  for(j in 1:J){

      for(jprime in (j+1):J){

          

            // CONTINUE HERE: Hope that I only need to specify one off-diagonal

            // and that the other will get filled in by the fact that Sigma is a covariance matrix

            Sigma[j, jprime] ~ uniform(0.0, sqrt(Sigma[j,j] * Sigma[jprime,jprime]));    

      }    

  }

}

[Bob Carpenter]

> On Jan 28, 2016, at 1:12 PM, Qian Zhang <qian.s...@gmail.com> wrote:
>
> Thanks Ben!
>
> So it turns out I have informative priors on entries of the Sigma matrix. The problem with using these sigma and Cholesky matrices you mention is I don't know how to derive priors on the sigma and Cholesky matrices from the priors on entries of Sigma, and it seems not-worth-it to go down that route.

It'll be worth it computationally, but the Jacobians
are painful. They're spelled out in the manual.

> I tried 2 approaches:
>
> (1) Declare Sigma as a cov_matrix in the parameters block, put priors on its entries in the model block, and rely on rejection of proposed Sigma values that are out of bounds (Code Excerpt 1 below)

This is bad for Stan. It causes its nice Hamiltonian flow to
devolved to a random walk. And if rejection's bad, it'll just lock.

> (2) Define Sigma as a cov_matrix in the transformed parameters block, define it in terms of sigma and Cholesky matrices, put priors on Sigma entries but not on sigma or Cholesky matrices (Code Excerpt 2 below)

You need to define it as a cov_matrix in the parameters block and put
priors on it. Otherwise, you'll have to deal with the Jacobian in
the prior on Sigma entries.

> Both approaches gave me draws. Approach 1 was a lot faster (~ a few minutes) than Approach 2 (~ 30 min), so I'm going with Approach 1.

You need to (a) measure convergence and time per effective sample size,
and you need to (b) make sure with simulated data you recover the right
results. Fast is meaningless if it's wrong!

> (I guess you put bounds in Approach 2 to encourage movement to a new position in the parameter space at the next draw, while in Approach 1 you tend to stay put at the last parameter draw because you propose something out of bounds and then reject it, so the only difference is mixing time, both real human time and time in terms of number of draws.)
>
> Question 0: Both approaches work, i.e. I'm not coding up something totally non-sensical either way? They should both implement MCMC procedures to draw from the same posterior.

See above.

> Question 1: Since Sigma is a cov_matrix, in the model block I only need to put priors on the upper diagonal entries of Sigma, and the lower diagonal entries will automatically be filled in with the upper diagonal draws?

It depends what you want the priors to do!
You are putting priors on the covariance matrix at this point.
We've already done the Jacobian for you.

What's really going on is that the unconstrained parameters
are the Cholesky factor L_Sigma and then it gets multipled
to give you Sigma = L_Sigma * L_Sigma'. So you need the
Jacobian if you're going to be putting priors on Sigma and
the declared parameter is L_Sigma.

>
> I got many warnings in Code Excerpt 1 of the following 3 types.
>
> Type 1:
>
> count
> validate transformed params: thetaB is -0.00107406, but must be greater than or equal to 0 4
>
>
> Type 2:
>
>
>
> Exception thrown at line 97: beta_log: Random variable is 1.00308, but must be less than or equal to 2

I find these truncated messages hard to read --- the constraint is that it must
be less than or equal to 1, and I believe the 2 means it showed up twice.

>
> Type 3:
>
> Warning messages:
> 1: There were 132367 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
> 2: Examine the pairs() plot to diagnose sampling problems

This means it's not working! You get numerical issues.

> Question 3: Is the following understanding correct? My understanding is Types 1 and 2 come from drawing Sigma entries that are out of bounds, or result in out-of-bounds thetaB, in which case such draws are rejected.

And 3. Whenever you get a numerical exception.

> Also, I never need to worry about this divergent transitions warning message, as you say here: “There are almost always some divergent transitions during the warmup phase but we don't really worry about them.

You should still worry about them. It means you'll bias
warmup. Starting with lower stepsize can mitigate in
some cases, but won't here because the models are broken.

> It is trying to find a stepsize during the warmup phase that achieves the target acceptance probability without diverging and without hitting the max_treedepth.” (https://groups.google.com/forum/#!topic/stan-users/hjYWqJZfSgw) But then why do you have this warning message at all?

The message indicates that you're going to get biases in your
estimates. We don't want to set our default parameters too far out
or it'll give the impression of hanging. We'd rather hit max
depth and then diagnose.

- Bob

> On Tuesday, January 26, 2016 at 9:23:52 PM UTC-8, Ben Goodrich wrote:
> On Tuesday, January 26, 2016 at 11:32:38 PM UTC-5, Qian Zhang wrote:
> J is an integer read in as data.
>
> I'd like to use a J x J covariance matrix Sigma in the model block that has diagonal entries bounded between [0,1], but off-diagonal entries bounded in terms of diagonal entries like so:
>
> 0 < Sigma[i,j] < sqrt(Sigma[i,i] * Sigma[j,j])
>
> How would I pull that off in STAN?
>

> I would do it by declaring a Cholesky factor of a correlation matrix
>
> cholesky_factor_corr[J] L_corr;
>
> and a vector of standard deviations that are bounded
>
> vector<lower=0,upper=1>[J] sigma;
>

> as parameters, and then
>
> diag_pre_multiply(sigma, L_corr);
>
> is a Cholesky factor of a covariance matrix whose diagonal elements would be between 0 and 1.
>
> 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+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.