column vector

[Yilong]

Hi:

I suppose someone might have asked this question before, however I could not find it in the mailing list. While setting a multivariate normal distribution, I need to declare a vector column from this distribution. Eg. 

real tmp[M,N]

tmp[,i] ~ multi_normal(mu, Sigma)    ( use R's syntax )

instead of 

tmp[i, ] ~ multi_normal(mu, Sigma) 

Currently, I solved this problem by transform parameters but it of course inefficient.  

Best,

Yilong

[Ben Goodrich]

On Tuesday, June 4, 2013 7:07:52 PM UTC-4, Yilong wrote:

Hi:

I suppose someone might have asked this question before, however I could not find it in the mailing list. While setting a multivariate normal distribution, I need to declare a vector column from this distribution. Eg. 

real tmp[M,N]

tmp[,i] ~ multi_normal(mu, Sigma)    ( use R's syntax )

If you are doing things like multi_normal(), you don't want to use a plain two-dimensional array (e.g. real tmp[M,N]). Often the fastest is to create a N-array of vectors

vector[M] tmp[N]

in which case you can access the i-th "column" with tmp[i].

The most intuitive (but slightly slower) way to do it from an R perspective is to declare a matrix

matrix[M,N] tmp;

in which case you can access (a copy of) the i-th column with col(tmp, i).

Ben

[Bob Carpenter]

On 6/4/13 7:42 PM, Ben Goodrich wrote:
> On Tuesday, June 4, 2013 7:07:52 PM UTC-4, Yilong wrote:
>
> Hi:
>
> I suppose someone might have asked this question before, however I could not find it in the mailing list. While
> setting a multivariate normal distribution, I need to declare a vector column from this distribution. Eg.
>
> real tmp[M,N]
>
> tmp[,i] ~ multi_normal(mu, Sigma) ( use R's syntax )
>
>
> If you are doing things like multi_normal(), you don't want to use a plain two-dimensional array (e.g. real tmp[M,N]).
> Often the fastest is to create a N-array of vectors
>
> |
> vector[M] tmp[N]
> |
>
> in which case you can access the i-th "column" with tmp[i].

If you absolutely don't want to mess with the input format,
you can define this array in the transformed data block.

> The most intuitive (but slightly slower) way to do it from an R perspective is to declare a matrix
>
> |
> matrix[M,N] tmp;
> |
>
> in which case you can access (a copy of) the i-th column with col(tmp, i).

The copy here is just an implementation detail that means it's
a bit slower. It will behave the same way in a Stan model
as using an array of vectors.

As far as efficiency, given that it's showing up in the context
of a multi-normal, it probably won't be too noticeable.

I like the array of vectors because it seems clearer conceptually
than a matrix, though of course, you can think of a matrix as
an array of vectors.

- Bob

[elong zhang]

Thanks Ben and Bob for your kindly reply. I actually simplified my question above. Due to the structure of my model, I both need column and row vectors to define a muti_normal and Dirichlet distribution. (That's the reason why I did not use vector[M] tmp[N] )   Below is my code that works in stan. But, I feel uncomfortable about the "transformed parameters" part. Even through the speed to draw samples from this model is quick, it waste memory to store the transform of matrix alpha. 

I feel col(tmp, i) should works in my situation. However, while I have a three dimension array like

real tmp[M,N,P]; 

how could I get row, column and ranks vector and some sub-matrix like tmp[1, , ] in a same framework ?  

|

data{

  int<lower = 1> M; // kinds of bacterias

  int<lower = 1> N; // Total bacterias

  int<lower = 1> T; // number of obs time

  int<lower = 1> N_case; // number of cases

  int<lower = 1> N_control; // number of control

  real time[T];

  int<lower = 0, upper = N> cases[N_case, T, M];      // count of bacterias in cases

  int<lower = 0, upper = N> control[N_control, T, M]; // count of bacterias in control

}

transformed data {

  vector[T] mu;

  cov_matrix[T] Sigma;

  for (i in 1:T)

    mu[i] <- 0;

  for (i in 1:T)

    for(j in 1:T)

      Sigma[i,j] <- exp(-pow(time[i] - time[j], 2)) +

                    if_else(i == j, 0.1, 0.0);

}

parameters {

  vector[T] alpha[M];

  simplex[M] theta_case[T];

  simplex[M] theta_control[T];

}

transformed parameters{

  vector<lower=0>[M] alpha_1[T];

  for(i in 1:T)

    for(j in 1:M)

      alpha_1[i,j] <- abs( alpha[j,i] );

}

model {

  for(i in 1:M){

    alpha[i] ~ multi_normal(mu, Sigma); 

  }

  for(i in 1:T){

    theta_case[i] ~ dirichlet( alpha_1[i] );

    theta_control[i] ~ dirichlet( alpha_1[i] );

  }

  for(i in 1:N_case)

   for(j in 1:T)

        cases[i,j] ~ multinomial(theta_case[j]);

  

  for(i in 1:N_control)

    for(j in 1:T)

        control[i,j] ~ multinomial(theta_control[j]);

}

|

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

--
Best,

Yilong

[Ben Goodrich]

On Wednesday, June 5, 2013 1:15:34 AM UTC-4, Yilong wrote:

Thanks Ben and Bob for your kindly reply. I actually simplified my question above. Due to the structure of my model, I both need column and row vectors to define a muti_normal and Dirichlet distribution. (That's the reason why I did not use vector[M] tmp[N] )   Below is my code that works in stan. But, I feel uncomfortable about the "transformed parameters" part. Even through the speed to draw samples from this model is quick, it waste memory to store the transform of matrix alpha. 

You can move the construction of alpha to the model block, in which case it won't be stored.
 

I feel col(tmp, i) should works in my situation. However, while I have a three dimension array like

real tmp[M,N,P]; 

how could I get row, column and ranks vector and some sub-matrix like tmp[1, , ] in a same framework ?

I would do an M-array of N by P matrices, in which case you can do row(tmp[i], n) or col(tmp[i], p).

Also,
 

transformed data {

  vector[T] mu;

  cov_matrix[T] Sigma;

  for (i in 1:T)

    mu[i] <- 0;

  for (i in 1:T)

    for(j in 1:T)

      Sigma[i,j] <- exp(-pow(time[i] - time[j], 2)) +

                    if_else(i == j, 0.1, 0.0);

}

You really want to also create Sigma_inverse in transformed data block so that here
 

model {

  for(i in 1:M){

    alpha[i] ~ multi_normal(mu, Sigma); 

  }

you can use multi_normal_prec(mu, Sigma_inverse) here instead of multi_normal(mu, Sigma). Your way internally inverts (basically) the fixed matrix Sigma M times per leapfrog step.

Ben

[Bob Carpenter]

I think you can simplify your model if I understand
it properly. As is, you're going to have horrible
problems with modes because of the abs().

On 6/5/13 1:15 AM, elong zhang wrote:

...

> parameters {
> vector[T] alpha[M];

...

> }
> transformed parameters{
> vector<lower=0>[M] alpha_1[T];
> for(i in 1:T)
> for(j in 1:M)
> alpha_1[i,j] <- abs( alpha[j,i] );
> }
> model {
> for(i in 1:M){
> alpha[i] ~ multi_normal(mu, Sigma);

...
> }
...

I may be missing something here, but I think what you
want is to just put a lower-bound on alpha itself:

parameters {
vector<lower=0> alpha[M];
...
}
model {
for (m in 1:M)
alpha[m] ~ multi_normal(mu, Sigma);
...
}

and then just replace all the alpha_1 uses with alpha.

It eliminates the multimodality and the need to transform.

Have you considered other priors here, such as a multivariate
logistic for the mean and something else for the scale?

// multivariate logistic prior for mean
alpha_mean_logodds[m] ~ multi_normal(mu,Sigma);
alpha_mean[m] <- softmax(alpha_mean_log_odds[m]);

// prior count (could use something stronger)
alpha_count[m] ~ cauchy(0,5);

// scale mean by count to get something to plug into Dirichlet
alpha[m] <- alpha_count[m] * alpha_mean[m];

We'll be speeding up these models soon with vectorization.
You can also speed up the model by using multi_normal precision
instead of covariance.

- Bob

[UD1989]

Sorry but doesnt

vector[M] tmp[N]

translate to matrix[N,M] ( and not matrix[M,N]) and wont tmp[i] access rows?

[Jonah]

If tmp is declared as a matrix then tmp[i] will access the ith row. But if tmp is declared as an array of column vectors then tmp[i] will access the ith column vector. 

[UD1989]

@Jonah : Thank you! That was a major point of confusion....thanks for clearing it up !

[Jonah]

Glad to help. I can see why that would be confusing. The way to keep it straight intuitively is to think of arrays and matrices/vectors as fundamentally different ideas and types of objects. An array is like a container. It can contain scalers, vectors, matrices, etc. If tmp is a container of vectors then it makes sense that tmp[i] is the ith vector rather than the ith row of the matrix you could conceivably make from the vectors if you happened to want to.

Jonah