Re: [stan-users] Estimation for a Markov Transitional Matrix

[Bob Carpenter]

I don't know what your n is, but if you want to compute
A^n, the best way to do it is to knock out powers of
2 by self-multiplying. But it's going to be a hugely
expensive operation with derivatives no matter how you do it.

I'm afraid I don't understand the rest of the model.
But it looks like your loops don't have the
right bounds unless D = K - 1.

> row_vector[D] gamma[(K-1)];
...
> for(i in 1:(K-1)){
> ...
> for (j in 1:(K-1)){
> gamma[i,j] ~ normal(0, 100);
> }


- Bob

> On Jul 17, 2016, at 8:29 PM, tomonn <tomoh...@gmail.com> wrote:
>
> Hello.
> I need your help.
> I tried to estimate a Markov chain matrix from a repeated cross-section data (i.e. a set of cross-section data at different time periods).
> Individual responses are three options (y_n = 1,2, or 3).
> Suppose individuals are identical, after controlling individual heterogeneities.
> The responses then changes according to a Markov chain.
> The Markov transitional matrix is given by
> M=
> [phi_{11} phi_{12} phi_{13}
> phi_{21} phi_{22} phi_{23}
> phi_{31} phi_{32} phi_{33}]
> My target is to estimate this matrix and the effects of individual heterogeneities on the matrix parameters.
>
>
> The following is my stan code with sessionInfo.
> I can compile this, but the computation does not go ahead.
> I don't know what causes the problem.
> Could you give me any suggestions?
> I also attach my simulation data, R code and stan code for your reference.
>
> Your kind assistance is highly appreciated.
>
> Best regards,
> Tom
>
>
> //// Stan Code "Gmail.stan"////
>
> functions{
> matrix matrix_pow(matrix a, int n);
> matrix matrix_pow(matrix a, int n) {
> if (n == 0)
> return diag_matrix(rep_vector(1, rows(a)));
> else
> return a * matrix_pow(a, n - 1);
> }
> }
> //
>
> data {
> int<lower=0> N; // Total # Observations
> int<lower=2> K; // Total # Responses K=3
> int<lower=1> D; // # of Individual Features
> int<lower=1, upper=3> y[N]; // Responses (Discrete Number/ 1,2,3)
> int<lower=0> period[N]; // Time Period Index
> vector[D] x[N]; // Explanatory Variables (Individual Features)
> simplex[K] Ps_init[N]; // Initial Distribution
> }
>
> parameters {
> simplex[K] phi[N,K]; // Matrix Elements
> real beta[(K-1)];
> row_vector[D] gamma[(K-1)];
> }
>
> transformed parameters {
> vector<lower=0>[K] alpha[N,K]; // Dirichlet Weight Parameters
> matrix<lower=0, upper=1>[K,K] Mmatrix[N]; // Markov Transitional Matrix
> simplex[K] Ps[N]; //
>
> for(n in 1:N){
> for(i in 1:K) {
> alpha[n, i, 2] <- 1;
> for(j in 1:(K-1)) {
> alpha[n,i, (2*j-1)] <- exp(beta[j] + gamma[j] * x[n]);
> }
> }
> }
>
>
> for (n in 1:N){
> for (i in 1:3){
> for (j in 1:3)
> Mmatrix[n,j,i] <- phi[n,i,j];
> }
> Ps[n] <- matrix_pow(Mmatrix[n], (period[n]-1)) * Ps_init[n];
> }
>
>
> }
>
> model {
> // Priors
>
> for(i in 1:(K-1)){
> beta[i] ~ normal(0,100);
> for (j in 1:(K-1)){
> gamma[i,j] ~ normal(0, 100);
> }
> }
> for(n in 1:N) {
> for (i in 1:K){
> phi[n,i] ~ dirichlet(alpha[n,i]);
> }
> }
>
> // Likelihood
> for (n in 1:N){
> y[n] ~ categorical(Ps[n]);
> }
> }
>
> ///// sessionInfo() /////
> R version 3.3.0 (2016-05-03)
> Platform: x86_64-apple-darwin13.4.0 (64-bit)
> Running under: OS X 10.11.5 (El Capitan)
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> loaded via a namespace (and not attached):
> [1] tools_3.3.0
>
>
>
> --
> 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.
> <GMailingList.zip>

[Bob Carpenter]

If you really need to evaluate matrix powers, you want to
knock of powers of two to the extent possible to reduce
the number of operations. That reduces a problem involving
N multiplications to one involving only log N multiplications.
And you want separate base cases for n == 0 and n == 1 to
avoid that needless constant unit matrix multiplication.

Also, you want to compute that matrix power only once. I'm
not even sure you really need to do it at all. You can just
estimate Markov models directly.

You also want to vectorize wherever possible.


Also, you are putting a lot of effort into estimating hierarchical
Dirichlet priors, but they tend not to be very powerful as priors
because there's no covariance information.

- Bob