Expression ill formed error in generated quantities block

Modeling
1

I have this simple model and the generated quantities is returning and error: expression is ill formed around this line: vector[N] x;. what is happening here?

x= seq(1:100)
y = 3+.7*x + rnorm(100,2,10)
group = rep(c("A","B"),each=50)
dat = data.frame(x = x, y = y, group = group)
dat
plot(x,y)
summary(lm(y~x))


library(rstan)
standat <- list(
  N = nrow(dat),
  y = dat$y,
  x = dat$x
)

stanmodelcode = '
data {
int<lower=1> N; 
real y[N];                                 
real x[N];  
}

transformed data{
}
parameters {
real intercept;  
real beta;                   
real<lower=0> sigma;  
}
transformed parameters {  // none needed
}

model {
real mu;

// priors
intercept~ normal(0,10);
beta ~ normal(0,10);
sigma ~ normal(0,10);

// likelihood
for(i in 1:N){
mu = intercept +  beta*x[i];
y[i] ~ normal(mu, sigma);

}

}

generated quantities{
  vector[N] x;
  vector[N] yhat;                // linear predictor
  real<lower=0> rss;             // residual sum of squares
  real<lower=0> totalss;         // total SS              
  real Rsq;                      // Rsq
  
  yhat = x * beta;
  rss = dot_self(y-yhat);
  totalss = dot_self(y-mean(y));
  Rsq = 1 - rss/totalss;


}
'

stan(model_code=stanmodelcode, data=standat, iter=10000, warmup=2000, chains = 1)
2

You cannot multiply a real array by a scalar. It will work if you declare

vector[N] x;
3

I just inserted that line and now that line returns the error. please see edit.

4

You have to declare x as a vector in the data block, rather than declaring it as a real array. Then you can use it in the model and generated quantities block.

5

Thank you. Just wondering why won’t real array work?

6

The multiplication operator is not defined for it.

7

The deeper reason why it doesn’t work is that in general, linear algebra operations are ambiguous for arrays. Stan separates row vectors and column vectors, so it can now that row vector times column vector produces a scalar result whereas vector times row vector produces a matrix result.

I have to admit it’s been confusing to our users familiar with R. Yet it’s exactly the ambiguity in R I wanted to avoid:

> c(1, 2) %*% c(10, 100)
     [,1]
[1,]  210

> t(c(1, 2)) %*% c(10, 100)
     [,1]
[1,]  210

Why do I get the same answer here? And more surprisingly, why do I get a different answer if I transpose the right-hand side?

> c(1, 2) %*% t(c(10, 100))
     [,1] [,2]
[1,]   10  100
[2,]   20  200

> t(c(1, 2)) %*% t(c(10, 100))
Error in t(c(1, 2)) %*% t(c(10, 100)) : non-conformable arguments
>

Even the basic typing is very confusing:

> is.array(c(1, 2))
[1] FALSE

> is.vector(c(1, 2))
[1] TRUE

> is.matrix(c(1, 2))
[1] FALSE

> is.array(t(c(1, 2)))
[1] TRUE

> is.vector(t(c(1, 2)))
[1] FALSE

> is.matrix(t(c(1, 2)))
[1] TRUE

Why when I transpose a vector do I get a matrix? Because there’s no row vector type, so it’s the only option. I based Stan’s design more on MATLAB, where row vectors and column vectors are two different things and the langauge doesn’t try to guess what you mean like in R when you multiply two column vectors together and it somehow decides to return a dot product rather than an error.