Known non-positive definite covariance matrix

[KK Sasa]

Dear all,

Because of the natural character of my model, I have to draw samples from multi_normal with known non-positive definite covariance matrix.
Like

x ~ multi_normal(mu1,sigma1);

where mu1 = c(0,0,0),

> sigma1
     [,1] [,2] [,3]
[1,]    2    1   -1
[2,]    1    2    1
[3,]   -1    1    2

However, the Stan requires the covariance matrix to be positive definite. Above method failed.
I know this can be done with the mvtnorm package of R language by issuing rmvnorm(n=1000, mean=c(0,0,0), sigma=sigma1).
I want to do the same thing as in R. But what can I do within Stan code?

Thank you in advance!!

[Michael Betancourt]

I'm not sure what R does, but if the covariance isn't positive-definite then the density (and any samples) are ill-defined.  In particular you get an unnormalizeable measure with most of the mass around infinity.

On Nov 22, 2013, at 9:49 AM, KK Sasa <genw...@gmail.com> wrote:

Dear all,

Because of the natural character of my model, I have to draw samples from multi_normal with known non-positive definite covariance matrix.
Like

randnum ~ multi_normal(mu1,sigma1);

where mu1 = c(0,0,0),

> sigma1
     [,1] [,2] [,3]
[1,]    2    1   -1
[2,]    1    2    1
[3,]   -1    1    2

However, the Stan requires the covariance matrix to be positive definite. Above method failed.
I know this can be done with the mvtnorm package of R language by issuing rmvnorm(n=1000, mean=c(0,0,0), sigma=sigma1).
I want to do the same thing as in R. But what can I do within Stan code?

Thank you in advance!!

--
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.
For more options, visit https://groups.google.com/groups/opt_out.

[KK Sasa]

Thanks for ur kind reply.
but i cannot avoid to do it.
Is it possible to call R function like rmvnorm in Stan?
or is it realistic for me to write a new distribution in Stan to do it?

[Ben Goodrich]

This variance-covariance matrix

     [,1] [,2] [,3]
[1,]    2    1   -1
[2,]    1    2    1
[3,]   -1    1    2

is positive semi-definite with rank 2. To use it with Stan, it depends on whether the multivariate normal distribution pertains to a likelihood or a prior. If it is a likelihood, you can use any two rows and columns of it to model data vectors of length 2. If it is a prior, then the best thing to do would be to calculate its Cholesky factor with pivoting (see help(chol) in R). In your case, L works out to

           [,1]     [,2] [,3]
[1,]  1.4142136 0.000000    0
[2,]  0.7071068 1.224745    0
[3,] -0.7071068 1.224745    0

If you declare z to be a vector of length 3 in the parameters block and put iid standard normal priors on its elements, then

L * z

is distributed multivariate normal with mean vector zero and the above variance-covariance matrix. At that point, you just have to prove (or hope) that the posterior is proper.

Ben

[Bob Carpenter]

That should totally go in the manual!

Could we define a low-rank matrix data type as
a parameter? (I realize you're still waiting for
Cholesky correlation, which is my highest coding
priority if I ever see my way clear to writing new
code rather than maintaining what's accumulating.)

But we could easily build it in for distributions.
That would just involve loosening the test on the
input. I'd be OK with that --- if a user wants to
validate a matrix, they can declare it as a cov_matrix.
(Soon to be renamed spd_matrix for "symmetric positive
definite".)

- Bob

[Ben Goodrich]

On Friday, November 22, 2013 12:21:41 PM UTC-5, Bob Carpenter wrote:

That should totally go in the manual!

It is in the manual in the subsection on multivariate reparameterizations (p.158). You may have mistakenly looked at the other two places in the Stan manual where multivariate reparameterizations are discussed. But that transformation is more general than just the multivariate normal. A similar thing holds for any elliptical distribution and it is one of the things that we could facilitate with a forward-looking transformation function.
 

Could we define a low-rank matrix data type as
a parameter?  

I don't know if it is prevalent enough to be a primitive type. A lot of people fit factor-type models to data or covariance matrices. That is why we needed a Cholesky factor that had more rows than columns so that the rank of  L * L'  is equal to the number of columns of L instead of the number of rows.

Ben

[KK Sasa]

Thanks for your helpful information.
I have tried it successfully in Stan.
But when I expand my model structure like this

> sigma1
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    2    1    1   -1   -1    0
[2,]    1    2    1    1    0   -1
[3,]    1    1    2    0    1    1
[4,]   -1    1    0    2    1   -1
[5,]   -1    0    1    1    2    1
[6,]    0   -1    1   -1    1    2

and run it in the same way in Stan, the result came out as follows

> cov(yy)
            [,1]        [,2]       [,3]        [,4]        [,5]        [,6]
[1,]  1.95608231  1.00459661 1.05309831 -0.95148570 -0.87882833  0.06265375
[2,]  1.00459661  2.05881940 1.08184909  1.05422279  0.06922895 -0.94372370
[3,]  1.05309831  1.08184909 2.06049579  0.02875078  1.00263318  1.01918009
[4,] -0.95148570  1.05422279 0.02875078  2.00570850  0.94805729 -1.00637744
[5,] -0.87882833  0.06922895 1.00263318  0.94805729  6.08575158  3.04545749
[6,]  0.06265375 -0.94372370 1.01918009 -1.00637744  3.04545749  6.81903572

As you see, the fifth and sixth elements are incorrect compared to the true covariance matrix.
Is something wrong with the covariance matrix?
Thank you!!

Ben Goodrich於 2013年11月22日星期五UTC+8下午9時29分13秒寫道：

[Bob Carpenter]

I'm not sure what you're asking here. How are sigma1
and cov(yy) related? Is cov(yy) an estimate of something
else based on a parameter in Stan?

- Bob

On 11/23/13, 10:28 PM, KK Sasa wrote:
> Thanks for your helpful information.

> I have tried it successfully in R.

> But when I expand my model structure like this
>
>> sigma1
> [,1] [,2] [,3] [,4] [,5] [,6]
> [1,] 2 1 1 -1 -1 0
> [2,] 1 2 1 1 0 -1
> [3,] 1 1 2 0 1 1
> [4,] -1 1 0 2 1 -1
> [5,] -1 0 1 1 2 1
> [6,] 0 -1 1 -1 1 2
>
>

> and run it in the same way in R, the result came out as follows

>
>> cov(yy)
> [,1] [,2] [,3] [,4] [,5] [,6]
> [1,] 1.95608231 1.00459661 1.05309831 -0.95148570 -0.87882833 0.06265375
> [2,] 1.00459661 2.05881940 1.08184909 1.05422279 0.06922895 -0.94372370
> [3,] 1.05309831 1.08184909 2.06049579 0.02875078 1.00263318 1.01918009
> [4,] -0.95148570 1.05422279 0.02875078 2.00570850 0.94805729 -1.00637744
> [5,] -0.87882833 0.06922895 1.00263318 0.94805729 6.08575158 3.04545749
> [6,] 0.06265375 -0.94372370 1.01918009 -1.00637744 3.04545749 6.81903572
>
>
> As you see, the fifth and sixth elements are incorrect compared to the true covariance matrix.
> Is something wrong with the covariance matrix?
> Thank you!!
>

> Ben Goodrich嚙踝蕭 2013嚙羯11嚙踝蕭22嚙踝蕭P嚙踝蕭嚙踝蕭UTC+8嚙磊嚙踝蕭9嚙踝蕭29嚙踝蕭13嚙踝蕭g嚙瘩嚙瘦

[KK Sasa]

Sorry for unclear description. yy is distributed from multivariate normal distribution with mean 0 and sigma1.
I used the method as Ben suggested  (p. 158) and run it in Stan.


Bob Carpenter於 2013年11月25日星期一UTC+8上午9時13分01秒寫道：

[Ben Goodrich]

Are there data involved? If so, the posterior distribution of something could easily be different from the prior distribution. That or it is possible that it didn't converge yet.

Ben

[KK Sasa]

In previous example, I have put in the data, but didn't use it.
This time, I simply drew samples from multivariate normal distribution and the R code is like

library(rstan)
model <- '
data{
  matrix[6,6] sigma1;
}
transformed data {
  matrix[6,6] L;
  vector[6] mu1;
  mu1[1] <- 0; mu1[2] <- 0; mu1[3] <- 0;
  mu1[4] <- 0; mu1[5] <- 0; mu1[6] <- 0;
  L <- cholesky_decompose(sigma1);
}
parameters{
  vector[6] randnum;
}
transformed parameters {
  vector[6] y;
  y <- mu1 + L * randnum;
}
model {
  randnum ~ normal(0,1);
}
'
sigma1 <- diag(6) * 2
sigma1[2,1] <- 1; sigma1[1,2] <- sigma1[2,1];
sigma1[3,1] <- 1; sigma1[1,3] <- sigma1[3,1];
sigma1[4,1] <- -1; sigma1[1,4] <- sigma1[4,1];
sigma1[5,1] <- -1; sigma1[1,5] <- sigma1[5,1];
sigma1[6,1] <- 0; sigma1[1,6] <- sigma1[6,1];
sigma1[2,3] <- 1; sigma1[3,2] <- sigma1[2,3];
sigma1[2,4] <- 1; sigma1[4,2] <- sigma1[2,4];
sigma1[2,5] <- 0; sigma1[5,2] <- sigma1[2,5];
sigma1[2,6] <- -1; sigma1[6,2] <- sigma1[2,6];
sigma1[3,4] <- 0; sigma1[4,3] <- sigma1[3,4];
sigma1[3,5] <- 1; sigma1[5,3] <- sigma1[3,5];
sigma1[3,6] <- 1; sigma1[6,3] <- sigma1[3,6];
sigma1[4,5] <- 1; sigma1[5,4] <- sigma1[4,5];
sigma1[4,6] <- -1; sigma1[6,4] <- sigma1[4,6];
sigma1[5,6] <- 1; sigma1[6,5] <- sigma1[5,6];
data <- list(sigma1=sigma1)
fit1 <- stan(model_code=model,data=data,
             iter=100000,chains=1,thin=1)

# extract samples
fit_ss <- extract(fit1, permuted = TRUE)
yy <- fit_ss$y

This output goes like

> cov(yy)
             [,1]          [,2]        [,3]         [,4]          [,5]         [,6]
[1,]  1.775847970  0.8776978664 0.871982398 -0.898150104 -0.9177753396  0.002816607
[2,]  0.877697866  1.7783345347 0.877739161  0.900636668  0.0007111106 -0.910376544
[3,]  0.871982398  0.8777391610 1.743159420  0.005756763  0.8589437393  0.873488914
[4,] -0.898150104  0.9006366683 0.005756763  1.798786772  0.9184864501 -0.913193151
[5,] -0.917775340  0.0007111106 0.858943739  0.918486450  5.3908461961  2.630386314
[6,]  0.002816607 -0.9103765442 0.873488914 -0.913193151  2.6303863139  6.247356763

It seems there is problem in convergence.

Ben Goodrich於 2013年11月25日星期一UTC+8上午11時36分49秒寫道：

[Bob Carpenter]

On 11/24/13, 11:55 PM, KK Sasa wrote:
> In previous example, I have put in the data, but didn't use it.
> This time, I simply drew samples from multivariate normal distribution and the R code is like
>
> library(rstan)
> model <- '
> data{
> matrix[6,6] sigma1;
> }
> transformed data {
> matrix[6,6] L;
> vector[6] mu1;
> mu1[1] <- 0; mu1[2] <- 0; mu1[3] <- 0;
> mu1[4] <- 0; mu1[5] <- 0; mu1[6] <- 0;

You can replace these six assignments with:

mu1 <- rep_vector(0.0, 6);

> L <- cholesky_decompose(sigma1);
> }
> parameters{
> vector[6] randnum;
> }
> transformed parameters {
> vector[6] y;
> y <- mu1 + L * randnum;

If the mu1 are all 0, you could just drop mu1
altogether for this model.

> }
> model {
> randnum ~ normal(0,1);
> }
> '
> sigma1 <- diag(pair_num) * 2

> sigma1[2,1] <- 1; sigma1[1,2] <- sigma1[2,1];
> sigma1[3,1] <- 1; sigma1[1,3] <- sigma1[3,1];
> sigma1[4,1] <- -1; sigma1[1,4] <- sigma1[4,1];
> sigma1[5,1] <- -1; sigma1[1,5] <- sigma1[5,1];
> sigma1[6,1] <- 0; sigma1[1,6] <- sigma1[6,1];
> sigma1[2,3] <- 1; sigma1[3,2] <- sigma1[2,3];
> sigma1[2,4] <- 1; sigma1[4,2] <- sigma1[2,4];
> sigma1[2,5] <- 0; sigma1[5,2] <- sigma1[2,5];
> sigma1[2,6] <- -1; sigma1[6,2] <- sigma1[2,6];
> sigma1[3,4] <- 0; sigma1[4,3] <- sigma1[3,4];
> sigma1[3,5] <- 1; sigma1[5,3] <- sigma1[3,5];
> sigma1[3,6] <- 1; sigma1[6,3] <- sigma1[3,6];
> sigma1[4,5] <- 1; sigma1[5,4] <- sigma1[4,5];
> sigma1[4,6] <- -1; sigma1[6,4] <- sigma1[4,6];
> sigma1[5,6] <- 1; sigma1[6,5] <- sigma1[5,6];

This kind of code is very error prone. I'd try
to keep the LHS ordered so that for sigma1[i,j]
definitions in the first column you have
i > j. Then you can make it symmetric with:

# 1 <= i < j <= 6
for (i in 1:5)
for (j in (i+1):6)
sigma1[i,j] <- sigma1[j,i];

Now that I look at it, this kind of code is very
error prone, too!

> data <- list(sigma1=sigma1)
> fit1 <- stan(model_code=model,data=data,
> iter=100000,chains=1,thin=1)
>
> # extract samples
> fit_ss <- extract(fit1, permuted = TRUE)
> yy <- fit_ss$y
>

What's the effective sample size and R-hat values?
I'd suggest running way fewer iterations and using
more chains.

> This output goes like

>
>> cov(yy)
> [,1] [,2] [,3] [,4] [,5] [,6]

> [1,] 1.775847970 0.8776978664 0.871982398 -0.898150104 -0.9177753396 0.002816607
> [2,] 0.877697866 1.7783345347 0.877739161 0.900636668 0.0007111106 -0.910376544
> [3,] 0.871982398 0.8777391610 1.743159420 0.005756763 0.8589437393 0.873488914
> [4,] -0.898150104 0.9006366683 0.005756763 1.798786772 0.9184864501 -0.913193151
> [5,] -0.917775340 0.0007111106 0.858943739 0.918486450 5.3908461961 2.630386314
> [6,] 0.002816607 -0.9103765442 0.873488914 -0.913193151 2.6303863139 6.247356763
>
>
> It seems there is problem in convergence.

What makes you think that?

- Bob

> Ben Goodrich嚙踝蕭 2013嚙羯11嚙踝蕭25嚙踝蕭P嚙踝蕭嚙瑾UTC+8嚙磕嚙踝蕭11嚙踝蕭36嚙踝蕭49嚙踝蕭g嚙瘩嚙瘦

>
> Are there data involved? If so, the posterior distribution of something could easily be different from the prior
> distribution. That or it is possible that it didn't converge yet.
>
> Ben
>
> On Sunday, November 24, 2013 8:57:09 PM UTC-5, KK Sasa wrote:
>
> Sorry for unclear description. yy is distributed from multivariate normal distribution with mean 0 and sigma1.
> I used the method as Ben suggested (p. 158) and run it in Stan.
>
>

> Bob Carpenter嚙踝蕭 2013嚙羯11嚙踝蕭25嚙踝蕭P嚙踝蕭嚙瑾UTC+8嚙磕嚙踝蕭9嚙踝蕭13嚙踝蕭01嚙踝蕭g嚙瘩嚙瘦

>
> I'm not sure what you're asking here. How are sigma1
> and cov(yy) related? Is cov(yy) an estimate of something
> else based on a parameter in Stan?
>
> - Bob
>
>

[Ben Goodrich]

On Sunday, November 24, 2013 11:55:18 PM UTC-5, KK Sasa wrote:

In previous example, I have put in the data, but didn't use it.
This time, I simply drew samples from multivariate normal distribution and the R code is like

library(rstan)
model <- '
data{
  matrix[6,6] sigma1;
}
transformed data {
  matrix[6,6] L;
  vector[6] mu1;
  mu1[1] <- 0; mu1[2] <- 0; mu1[3] <- 0;
  mu1[4] <- 0; mu1[5] <- 0; mu1[6] <- 0;

  L <- cholesky_decompose(sigma1);
}

The cholesky_decompose function in Stan does not use pivoting and hence that L is calculated very inaccurately with Sigma is singular. That is why I suggested to calculate L in R with pivoting.

Ben

Ben Goodrich於 2013年11月25日星期一UTC+8上午11時36分49秒寫道：

[Bob Carpenter]

On 11/25/13, 12:32 AM, Ben Goodrich wrote:

> The cholesky_decompose function in Stan does not use pivoting and hence that L is calculated very inaccurately with
> Sigma is singular. That is why I suggested to calculate L in R with pivoting.

Is there a way we can make what we're doing more robust
by pivoting. As you know, my linear algebra fu is weak,
at best.

It shouldn't be hard to change the behavior, as it's
a simple function that's already unit tested.

- Bob

[KK Sasa]

Many thanks for Bob's suggestion on coding.

In present stage, I am simply testing the reliability of the the cholesky_decompose function in Stan and make sure it will work if the SIGMA matrix will be parameters to be estimated (with necessary model constraints) in next stage.
If cholesky_decompose function is not stable in this condition, my complete model seems not estimatable in Stan. :(

- cw

Bob Carpenter於 2013年11月25日星期一UTC+8下午1時35分13秒寫道：

[Bob Carpenter]

You don't need to use the Cholesky decomposition --- it's
just an efficiency issue. But then the model may not be
practially estimable if it's not efficient enough.

- Bob

On 11/25/13, 12:46 AM, KK Sasa wrote:
> Many thanks for Bob's suggestion on coding.
>
> In present stage, I am simply testing the reliability of the the cholesky_decompose function in Stan and make sure it

> will works if the SIGMA matrix will be parameters to be estimated (with necessary model constraints) in next stage.

> If cholesky_decompose function is not stable in this condition, my complete model seems not estimatable in Stan. :(
>
> - cw
>

> Bob Carpenter嚙踝蕭 2013嚙羯11嚙踝蕭25嚙踝蕭P嚙踝蕭嚙瑾UTC+8嚙磊嚙踝蕭1嚙踝蕭35嚙踝蕭13嚙踝蕭g嚙瘩嚙瘦

[KK Sasa]

I am a bit confused, Bob. I indeed need to draw samples from such singular SIGMA matrix in estimation.
I have no idea what to do it if not to use the Cholesky decomposition in Stan.
Or I can use eigenvectors and eigenvalues to generate random samples in Stan?

- cw

Bob Carpenter於 2013年11月25日星期一UTC+8下午2時05分26秒寫道：

[Bob Carpenter]

Sorry -- I forgot the context. We don't have
a good way to do that given that we assume
full rank covariance matrices (positive def).

- Bob

On 11/25/13, 1:21 AM, KK Sasa wrote:
> I am a bit comfused, Bob. I indeed need to draw samples from such singular SIGMA matrix in estimation.
> I have no idea what to do it if not to use the Cholesky decomposition in Stan.
>
> - cw
>
> Bob Carpenter嚙踝蕭 2013嚙羯11嚙踝蕭25嚙踝蕭P嚙踝蕭嚙瑾UTC+8嚙磊嚙踝蕭2嚙踝蕭05嚙踝蕭26嚙踝蕭g嚙瘩嚙瘦

>
> You don't need to use the Cholesky decomposition --- it's
> just an efficiency issue. But then the model may not be
> practially estimable if it's not efficient enough.
>
> - Bob
>
>

[Ben Goodrich]

I think you just need to calculate L in R with pivoting and pass that to Stan as data. It is true that at the moment Stan does not have any matrix decomposition algorithms that do pivoting.

Ben

On Monday, November 25, 2013 12:46:40 AM UTC-5, KK Sasa wrote:

Many thanks for Bob's suggestion on coding.

In present stage, I am simply testing the reliability of the the cholesky_decompose function in Stan and make sure it will work if the SIGMA matrix will be parameters to be estimated (with necessary model constraints) in next stage.

If cholesky_decompose function is not stable in this condition, my complete model seems not estimatable in Stan. :(

- cw

Bob Carpenter於 2013年11月25日星期一UTC+8下午1時35分13秒寫道：

[Ben Goodrich]

Well, with pivoting it permutes the rows and columns, so it would lead to a lot of user errors. We do use the 3-factor Cholesky with pivoting internally in some case.

Basically, for fixed covariance matrices, you can do the decomposition externally and pass it as data. Or you could do it in Stan via an eigendecomposition instead of the Cholesky decomposition. In the parameterized case, everyone would be better off working with the Cholesky factor primitive types.

Ben


 

[Bob Carpenter]

I think we're answering too many questions!

KK Sasa wants the matrix to be parameters, so the
decomposition can't be done in R.

I was also trying to ask if we could make the
Cholesky decomposition more accurate in Stan by
using pivoting. We need to do more work on
testing positive-definiteness, too, because we
still have the outstanding issue of high-dimensional
correlation and covariance matrix parameters being
rejected for not being positive definite.

- Bob

On 11/25/13, 2:06 PM, Ben Goodrich wrote:
> I think you just need to calculate L in R with pivoting and pass that to Stan as data. It is true that at the moment
> Stan does not have any matrix decomposition algorithms that do pivoting.
>
> Ben
>
> On Monday, November 25, 2013 12:46:40 AM UTC-5, KK Sasa wrote:
>
> Many thanks for Bob's suggestion on coding.
>
> In present stage, I am simply testing the reliability of the the cholesky_decompose function in Stan and make sure
> it will work if the SIGMA matrix will be parameters to be estimated (with necessary model constraints) in next stage.
> If cholesky_decompose function is not stable in this condition, my complete model seems not estimatable in Stan. :(
>
> - cw
>

> Bob Carpenter嚙踝蕭 2013嚙羯11嚙踝蕭25嚙踝蕭P嚙踝蕭嚙瑾UTC+8嚙磊嚙踝蕭1嚙踝蕭35嚙踝蕭13嚙踝蕭g嚙瘩嚙瘦

>
>
>
> On 11/25/13, 12:32 AM, Ben Goodrich wrote:
>
> > The cholesky_decompose function in Stan does not use pivoting and hence that L is calculated very
> inaccurately with
> > Sigma is singular. That is why I suggested to calculate L in R with pivoting.
>
> Is there a way we can make what we're doing more robust
> by pivoting. As you know, my linear algebra fu is weak,
> at best.
>
> It shouldn't be hard to change the behavior, as it's
> a simple function that's already unit tested.
>
> - Bob
>

[KK Sasa]

I learn that the correlated random samples from the singular matrix Sigma can be obtained via eigendecomposition as well.
So, I tried [V,D] = eig(Sigma) in MATLAB, then obtained the correlated random samples x = A*Z with A = V*sqrt(D); z = randn(6,10000);

sigma =

     2     1     1    -1    -1     0

     1     2     1     1     0    -1

     1     1     2     0     1     1

    -1     1     0     2     1    -1

    -1     0     1     1     2     1

     0    -1     1    -1     1     2

V =
         0       -0.5774   -0.4082   -0.5774         0        -0.4082
   -0.3536    0.5774   -0.2041   -0.5774   -0.3536    0.2041
    0.3536   -0.0000    0.6124   -0.5774    0.3536    0.2041
   -0.3536   -0.5774    0.2041    0.0000   -0.3536    0.6124
    0.3536   -0.0000   -0.6124   -0.0000    0.3536    0.6124
   -0.7071    0             0             0             0.7071     0
D =
         0         0         0         0         0         0
         0         0.00    0         0         0         0
         0         0         0.00    0         0         0
         0         0         0         4.00    0         0
         0         0         0         0         4.00    0
         0         0         0         0         0         4.00
A = V*sqrt(D) =
         0   -0.0000   -0.0000   -1.1547    0            -0.8165
         0    0.0000   -0.0000   -1.1547   -0.7071    0.4082
         0   -0.0000    0.0000   -1.1547    0.7071    0.4082
         0   -0.0000    0.0000    0.0000   -0.7071    1.2247
         0   -0.0000   -0.0000   -0.0000    0.7071    1.2247
         0         0         0             0            1.4142      0

Finally, obtain x = A*randn(6,10000)
So, cov(x') =
    2.0128    0.9696    0.9890   -1.0432   -1.0238    0.0194
    0.9696    1.9679    0.9710    0.9982    0.0014   -0.9968
    0.9890    0.9710    1.9866   -0.0180    0.9976    1.0156
   -1.0432    0.9982   -0.0180    2.0414    1.0252   -1.0162
   -1.0238    0.0014    0.9976    1.0252    2.0214    0.9962
    0.0194   -0.9968    1.0156   -1.0162    0.9962    2.0124
which is approximately equal to the Sigma matrix.

My question is whether I can do the same thing in Stan as in MATLAB when the singular Sigma matrix is parameterized (fixed singular in estimation)?

However, I found that the p. 273 in Stan manual reads

vector eigenvalues_sym(matrix A)
     The vector of eigenvalues of a symmetric matrix A in descending order
matrix eigenvectors_sym(matrix A)
     The matrix with the eigenvectors of symmetric matrix A

The eigenvalues are sorted by default, so are eigenvectors?
In addition, I found the R function 'eigen'  sorts the eigenvalues in decreasing order, too.



Ben Goodrich於 2013年11月26日星期二UTC+8上午3時11分19秒寫道：

[Ben Goodrich]

On Monday, November 25, 2013 9:20:02 PM UTC-5, KK Sasa wrote:

My question is whether I can do the same thing in Stan as in MATLAB when the singular Sigma matrix is parameterized (fixed singular in estimation)?

However, I found that the p. 273 in Stan manual reads

vector eigenvalues_sym(matrix A)
     The vector of eigenvalues of a symmetric matrix A in descending order
matrix eigenvectors_sym(matrix A)
     The matrix with the eigenvectors of symmetric matrix A

The eigenvalues are sorted by default, so does eigenvectors?
In addition, I found the R function 'eigen' also sorts the eigenvalues in decreasing order, too.

Yes, this method works for any matrix M such that M * M' = A. In this case, M is equal to the matrix of eigenvectors post-multiplied by a diagonal matrix with square roots of eigenvalues, so you have to manually intervene to make sure the last eigenvalue(s) are not slightly negative when they are exactly zero theoretically. You can do it inside Stan that way. The ordering doesn't matter if you are just multiplying by iid standard normals.

Ben


Ben

[KK Sasa]

Ben Goodrich於 2013年11月26日星期二UTC+8上午11時21分47秒寫道：

Thanks for your feedback, Ben. But I indeed need the ordering of elements of eigenvectors and eigenvalues to be fixed as I showed previously in MATLAB, because it will be much more difficult to estimate Sigma matrix without knowing the ordering of elements of Sigma.
The ordering of elements of Sigma is corresponding to my data matrix.
Could we modify the eigenvalues_sym and eigenvectors_sym without ordering?

[Ben Goodrich]

On Tuesday, November 26, 2013 12:11:00 AM UTC-5, KK Sasa wrote:

Thanks for your feedback, Ben. But I indeed need the ordering of elements of eigenvectors and eigenvalues to be fixed as I showed previously in MATLAB, because it will be much more difficult to estimate Sigma matrix without knowing the ordering of elements of Sigma.

The ordering of elements of Sigma is corresponding to my data matrix.
Could we modify the eigenvalues_sym and eigenvectors_sym without ordering?

That would require carrying a patch against Eigen, which we are reluctant to do. The order of the eigenvalues and eigenvectors doesn't affect the elements of Sigma. So, you could just change your code so that it assumes the zero eigenvalues come first. Or you could do the eigendecomposition in R beforehand, order things however you want, and pass that to Stan as data.

Ben
 

[KK Sasa]

Thanks for your kind reply. The eigendecomposition in Stan is useful and correct.
Yet, I run into a very simple problem. The normal sampling distribution seems problematic? or I misused it?
the R code like this

model <- '
data{
  int<lower=1> N;
  int<lower=1> I;
}
parameters{
  matrix[N,I] randnum;
}
model {
  for (n in 1:N) {
     for (i in 1:I) {
      randnum[n,i] ~ normal(0,1);
  }}
}
'
Then
fit_ss <- extract(fit1, permuted = TRUE) # extract results

# calculate the covariance matrix of randnum
r <- summary(fit1,pars = "randnum")
r_est <- r$summary[,1] # read estimate 'mean'
dim(r_est) <- c(pair_num,people) # reshape, place it from row
r_est <- t(r_est)
cov(r_est)

I got
> cov(r_est)

              [,1]         [,2]          [,3]         [,4]         [,5]         [,6]

[1,]  1.3485412075 -0.010084119 -0.0008423765 -0.007074001  0.002629817 -0.005378524
[2,] -0.0100841190  1.328471452  0.0346705190 -0.017408832  0.001335260 -0.002542019
[3,] -0.0008423765  0.034670519  1.3281833322 -0.006896448 -0.011058352  0.008618782
[4,] -0.0070740009 -0.017408832 -0.0068964476  1.297968318 -0.012098921  0.006230886
[5,]  0.0026298166  0.001335260 -0.0110583516 -0.012098921  1.345455370 -0.008625754
[6,] -0.0053785239 -0.002542019  0.0086187819  0.006230886 -0.008625754  1.309925881

Why the diagonal variance is not equal to 1?

Ben Goodrich於 2013年11月26日星期二UTC+8下午1時47分45秒寫道：

[Ben Goodrich]

On Thursday, November 28, 2013 9:06:17 AM UTC-5, KK Sasa wrote:

Thanks for your kind reply. The eigendecomposition in Stan is useful and correct.
Yet, I run into a very simple problem. The normal sampling distribution seems problematic? or I misused it?

Try this:

normal <- stan(file = "normal.stan", data = list(N = 100, I = 6))
randnum <- extract(normal, pars = "randnum")[[1]]
dim(randnum) # 4000 100 6
> round(matrix(rowMeans(apply(randnum, 1, cov)), 6, 6), 7)

           [,1]       [,2]       [,3]       [,4]       [,5]       [,6]

[1,]  1.0027809 -0.0049003  0.0018023 -0.0001027  0.0014064 -0.0019516
[2,] -0.0049003  0.9988229  0.0027084 -0.0016156 -0.0042573 -0.0013159
[3,]  0.0018023  0.0027084  1.0038629  0.0016223  0.0006535 -0.0037108
[4,] -0.0001027 -0.0016156  0.0016223  1.0064975 -0.0011702  0.0002508
[5,]  0.0014064 -0.0042573  0.0006535 -0.0011702  1.0062963 -0.0003245
[6,] -0.0019516 -0.0013159 -0.0037108  0.0002508 -0.0003245  1.0082776

Ben

[KK Sasa]

Ben Goodrich於 2013年11月26日星期二UTC+8上午11時21分47秒寫道：

Thank you very much. Your feedback is extremely helpful for me!

Here comes a new question.
If I use the Eigenvector decomposition and want to truncate correlated samples (e.g., y > 0), is it possible to do this in Stan?
I know Stan offer the normal() T[], but it cannot perform when sigma is not positive definite.
So, this is why I resort to Eigenvector decomposition.

 

Ben