Constraints on the factor loading matrix

[maxwell]

Dear Stan users,

I am currently working with a factor model. It is know that for identifiability of the factor model, it is imperative to impose a constraint on the factor loading matrix. One way to achieve that is to impose that the factor loading matrix is a full rank matrix by choosing a lower triangular matrix where the elements of the diagonal are forced to be strictly positive.  My questions are:

1) Can such model be implemented in Rstan? If yes, can you please point me to an example or reference(s).

2)  Is there a lower triangular matrix variable type predefined in Rstan? 

Thank you.  

[Bob Carpenter]

1) Stan manual, section 22.4 on matrix and vector types.

2) If you define:

cholesky_factor_cov[K] L;

then L is a lower-triangular matrix with the diagonals constrained
to be positive.

matrix[K,K] Sigma;

Sigma <- multiply_lower_tri_self_transpose(L);

is the efficient way to set Sigma <- L * L’; See the specialized
products subsection of section 32.2, Matrix Arithmetic Operations.

It’s nice to be able to just answer “yes” once in a while!

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

[Marcus Brubaker]

Bob,

In these models it's often the case that you want non-square factor loading matrices.  Is that implemented or must L be square?  I don't remember any more...

Maxwell,

If Stan doesn't currently support non-square versions of that type you can also do this by hand with:

parameters {

  vector[D] L_nondiag; // D is the number of lower diagonal entries (i don't know the formula off the top of my head...)

  vector<lower=0>[N] L_diag; // N is the number of factors

}

transformed parameters {

  matrix[M,N] L; // M is the number of rows and N is the number of factors

  int index = 1;

  for (i in 1:N) {

    L[i,i] <- L_diag[i];

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

      L[j,i] <- L_nondiag[index];

      L[i,j] <- 0;

      index <- index + 1;

    }

  }

}

Cheers,

Marcus

[Bob Carpenter]

Nope — the matrix doesn’t have to be square. I’ll clarify that in
the doc for the next release.

It was your (Marcus’s) idea to generalize — I still count matrix ops on
my fingers.

And there are specialized efficient auto-diffs for it.

- Bob

[Roger Zoh]

Indeed, it works for me! The matrix does not have to be square.

Thanks much guys!

On Thu, Feb 20, 2014 at 6:51 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Nope -- the matrix doesn't have to be square.  I'll clarify that in

the doc for the next release.

It was your (Marcus's) idea to generalize -- I still count matrix ops on

--

Roger Zoh, Ph.D.

Department of Statistics

and Nutrition

Blocker 525B and Kleberg 318

Texas A&M University

[Joachim Vandekerckhove]

Could you post an example of that working?  I'm trying to get Stan to do exploratory factor analysis but I keep running into convergence issues.  Even a simple model seems to fail.

For example, in this (unconstrained) model, I expected the matrix M to converge easily (and then I'd go on to apply penalty terms to L to add constraint), but it doesn't.

data { 

  int<lower=1> A;  # Number of latent factors

  int<lower=1> P;  # Number of participants

  int<lower=1> T;  # Number of manifest variables

  matrix[T,P] X;  # Raw data

}

parameters {

  matrix[A,P] F;  # Factor scores

  matrix[T,A] L;  # Loading matrix

}

transformed parameters {

  matrix[T,P] M;  # Prediction

  M <- L*F;

}

model { 

  for (t in 1:T)

    X[t] ~ normal(M[t], 1);  # Likelihood

[小杉考司]

Dear Bob and Marcus,

My name is Koji. 

I am a beginner of rstan, but have become to love this gradually.

I write Marcus's code, but error message showed that  "int " can not be declared in transformed parameters block.

What can I do for this problem?

P.S. Using choleskey_factor_cov is very elegant way. But the restriction of positive diagonal is not good for my data. For EFA, the restriction is need for only lower triangle.

2014年2月21日金曜日 6時04分06秒 UTC+9 Marcus:

[Bob Carpenter]

It would help if you sent the model itself.

You can define an integer local variable in the model itself,
or you can define an integer generated quantity.

I believe you need to have the positive diagonal on the Cholesky
decomposition in order for the product to be positive definite.
Do you want zeroes so that it's not full rank and only positive
semi-definite?

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

> For more options, visit https://groups.google.com/d/optout.

[小杉考司]

Dear Bob,

Thank you for your kindly message.
I attach my model. I just want to make EFA for Kendall data.

When I use cholesky_factor_cov vectors, the model goes fine.
But the diagonal of factor loadings matrix L(M by N) could have a
negative value.
I think that its not unusual way to use choeskey-vec, so I want to
make diag as zero.

> 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/HmkqJkT8i3o/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.

--
*************************************************
Don't Feel, Think!!
　小杉 考司 ( Dr.K.E.K )
e-mail : kos...@yamaguchi-u.ac.jp (Office)
　kosu...@gmail.com (Private)
WebSite: https://sites.google.com/site/kosugitti/
*************************************************

[小杉考司]

Bob,

I had many mistakes on my code about the size of matrix, but I got
solution for my problem.

My friend told me to make brackets, I can declare local values as int.

Now I been struggling on another issue, it would stop R at compile
time, but I think cause it probably points to elements outside of the
matrix.

I in going to work towards resolution.
Thank you very much!