functions {

  matrix exp_cov_matrix (matrix x, real sill, real range, real nugget){

    int S;
    matrix[dims(x)[1],dims(x)[1]] cov;

    S <- dims(x)[1];


    // off-diagonal elements
    for (i in 1:(S-1)) {
      for (j in (i+1):S) {
        # beta0 = partial sill = scale^2
        # beta1 = range = phi
        cov[i,j] <- sill * exp(-1.0/range * x[i,j]);
        cov[j,i] <- cov[i,j];
      }
    }

    // diagonal elements
    for (k in 1:S)
      cov[k,k] <- sill + nugget;

    return cov;

  }

}

data {
  int<lower=1> S; // number of locations
  int<lower=1> P; // number of space predictors (including intercept)
  vector[S] y; // observations
  matrix[S,S] dist; //distance matrix stations 
  matrix[S,P] covars;
}

transformed data{
  vector[S] zero_vector;
  for(s in 1:S)
    zero_vector[s] <- 0.0;
}

parameters {

  real<lower=0> sill;
  real<lower=0> range;
  real<lower=0> nugget;

  // regression coefficients
  vector[P] betas;

  real<lower=0> sigma;
} 

model {    

  // compute covariance matricies on the fly
  // i.e. dont estimate posteriors for the whole matrix
  matrix[S,S] cov;
  vector[S] resid;
  vector[S] mu;

  // covariance and inverse covariance based on knot locations 
  cov <- exp_cov_matrix(dist, sill, range, nugget);

  # exp link function
  mu <- covars * betas;
  resid <- y - mu;

  y ~ normal(mu,sigma);

  // The Y(s) = XB + w(s)
  // µ = Xß + w(s) ; w(s) ~ MVN(0,∑)
  resid ~ multi_normal(zero_vector, cov);

  // priors on space parameters
  sill ~ normal(0,100);
  range ~ normal(0,100);
  nugget ~ normal(0,100);

  // regression coefficients (weakly informative priors)
  betas ~ normal(0, 100);
}