data {
  int<lower=0> N;
  real y[N];
  matrix[N,N] A;
  vector[N] mu_u_prior;
}


transformed data {

matrix[N,N] A_inv;

A_inv <- inverse(A);

}


parameters{
 real beta0;
 vector[N] u;	
 real<lower=0> sigmasq_e;
 real<lower=0> sigmasq_u;	
} 


model {


  ##-----------------------------------------------
  ## likelihood	
  for(i in 1:N){
       y[i] ~ normal(beta0 + u[i], sqrt(sigmasq_e)); 
     }  


  ##---------------------------------
  ## priors
  
 lp__ <- lp__ - 0.5 * N * log(sigmasq_u); // non-constant part of log(det(sigmasq_u * A)^-0.5)
 lp__ <- lp__ - (u' * A_inv * u) / ( 2 * sigmasq_u); // log of kernel of multivariate normal

  beta0 ~ normal(0,1000);

  sigmasq_e ~ scaled_inv_chi_square(2,0.25);
  sigmasq_u ~ scaled_inv_chi_square(2,0.25);	

}