Vectorization of real valued parameter

Modeling
1

I am struggling with this rather simple statement in my STAN code. Is there a way to vectorize a real-valued parameter? Normal(mu,sd) returns a real-valued array. However, as you can see I have defined x to be a vector[n_obs]. What is puzzling me is the fact that in the statement

X=append_col(ones,x);

it does treat x to be a vector and hence allows me to use append_col. However, in the statement,

P=(x’)*x

it throws me an error of type mismatch, where the base type is a matrix and the right-hand side is a real. Can you guys help me try to understand how the same value is treated differently in two subsequent statements and any suggestions on how to avoid this?

The STAN code is attached here and I realize many of the statements in the transformed data block are redundant and be clubbed together, but I have broken down the block as such since many several people are working on this one project and have varying knowledge of STAN. Thank you in advance.

> transformed data{
  vector[n_obs] x;
  int r=0;
  matrix [n_obs,r] phi;
  vector[r] zeros;
  vector[n_obs] ones;
  matrix [n_obs,n_cov] X;
  matrix [n_obs,n_obs] G;
  matrix [n_obs,n_obs] P;
  matrix [n_obs,n_obs] M;
  matrix [n_obs,n_obs] Qt;
  matrix [r,r] C;
  matrix [r,r] D;
  matrix [r,r] Z;
  matrix [r,r] A;
  matrix [n_obs,1] x1;
  matrix[r, r] L = cholesky_decompose(A);
  vector[n_obs] values;
  vector [r] V1;
  vector[r] delta;
  zeros = rep_vector(0, r);
  delta = rep_vector(0, r);
  ones = rep_vector(0, n_obs);
  X=append_col(ones,x);
  P=(x')*x;
  G=(I-x*inverse(P)*(x'))*W*(I-x*inverse(P)*(x'));
  values = eigenvalues_sym(G);
  M=eigenvectors_sym(G);
  for(i in 1:n_obs)
  {
    if(values[i]>0)
    {r=r+1;}
  }
  phi=block(M,1,1,n_obs,r);
  Qt=I-0.1571172*W;
  C=quad_form(Qt,phi);
  V1=eigenvalues_sym(C);
  for(i in 1:r)
  {
    if(V1[i]>0) delta[i]=V1[i];
    else delta[i]=0;
  }
  D=diag_matrix(delta);
  Z=eigenvectors_sym(C);
  A=quad_form(D,Z);
}
parameters {
  vector[n_cov] beta;
  vector[r] v;
}
transformed parameters{
  vector[n_obs] eta;
  vector[n_obs] F;
  F=X*beta;
  eta = phi*v;
}
model
{
  x ~ normal(mu_X,Sd_X);
  y ~ poisson_log(F+eta);
  beta~ normal(0,10);
  v ~ multi_normal_cholesky(zeros, L);
}
2

I tried to simplify the code to just the parts that return the unexpected result. As I see it, if you want P to be an n_obs x n_obs matrix you should generate it with

P = x * x'

because variables of the type vector are column vectors.

Here is my code and the result of printing out the variables X, P and P2 (a variable I invented)

MODEL <- "
data {
    int<lower=0> N;
    vector[N] x;
}
transformed data {
    vector[N] ones;
    matrix [N,2] X;
    matrix [N,N] P;
    real P2;
    ones = rep_vector(0, N);
    X = append_col(ones,x);
    P = x * x';
    P2 = x' * x;
    print(X);
    print(P);
    print(P2);
}
"
vec <- c(1,2,3,4)
dataList = list(N = length(vec), x = vec)
Fit <- stan(model_code = MODEL, data = dataList, algorithm = "Fixed_param", chains = 1, iter = 1)

[[0,1],[0,2],[0,3],[0,4]]  <- Result of print(X)
[[1,2,3,4],[2,4,6,8],[3,6,9,12],[4,8,12,16]]  <- Result of print(P)
30  <- Result of print(P2)

I am sorry if I misunderstood your question. I am fairly new to Stan myself.

1 Like
3

Thank you very much. I feel utterly embarrassed for having worked on this for so long and not checking the dimensions. Thank you again.

4

The matrix x * (x') is a rank 1 and therefore not invertible. I think what you want here is

P = x*(x')/dot_self(x);
G = quad_form(W, diag_matrix(rep_vector(1, n_obs)) - P);
5

@nhuurre x is a column vector of size n_obs. Then, shouldn’t x*(x') gives us a n_obs x n_obs matrix? All these two statements are trying to do is implement the Moran’s I operator, which is basically,

 G=(I-X(X'X)^{-1}X')*W*(I-X(X'X)^{-1}X')
6

Yes, XX' is an n_obs by n_obs matrix, but it’s not invertible. Moran’s operator inverts X'X which is a scalar (and equal to dot_self(X)).

7

@nhuurre I got your point. But when I incorporate that and run the code, it throws the following error.

Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) : 
  Exception: eigenvalues_sym: m is not symmetric. m[1,2] = nan, but m[2,1] = nan  (in 'model3bec6b921450_count_basis_new' at line 39)
8

Weird. Does it work the other way? Maybe print intermediate values to see were it goes NaN.

Actually, I just noticed you have vector x and matrix X. Which one are you supposed to use to build Moran’s?

9

@nhuurre vector x is what we are supposed to use. matrix X is for a different purpose. As for whether it works, we have previously worked on this as a mixed effect model, where the vector x was considered fixed and fed into the STAN modeling block as a data and it has worked just fine. However, now we are trying to change it to a random effect model and hence have included vector x as a modeling parameter. This entire basis function set up has worked in the previous case, so we know the theory is there.