I’m trying to specify a Gaussian process model for several functional observations with shared mean and observation-level random effects. The model I have coded is a bit unwieldy since all the functional observations are stored in a vector. It would be nicer if I could have them in a matrix with each row belonging to one observation, but the shared mean makes this difficult for me to figure out how to do it.

I have several functional observations y_i, i=1,..,N over a common grid of time values t. The model is

y_i(t) = \mu(t) + \theta_i + \epsilon_i(t),

\theta_i \sim N(0, \sigma^2_\theta),

\mu \sim GP(0,K_{\tau_1,l_1}),

\epsilon_i \sim GP(0,K_{\tau_2,l_2} + \sigma^2 \delta),

where K are squared exponential covariance functions and \delta the identity function. Observations have covariance

Cov(y_i(t),y_j(t')) = K_\mu(t,t') + \delta_{ij}K_\epsilon(t,t') + \delta(t,t')\sigma^2.

That is, each observation has a common mean \mu(t), a structured error \epsilon_i(t) that is autocorrelated within observations and independent between observations, and unstructured error common to account for noise in measurement.

The priors are irrelevant to my problem.

My model block is

```
data {
int<lower=1> n_obs; // number of functional observations
int<lower=1> len; // length of each functional observation
int ID[(n_obs*len)]; // ID of each observational unit
vector[(n_obs*len)] y; // Functional observations
real x[(n_obs*len)]; // Time values of each observation
}
transformed data {
int N = n_obs*len; // Total number of data points
}
parameters {
real<lower=0> l1;
real<lower=0> t1;
real<lower=0> l2;
real<lower=0> t2;
real<lower=0> sigma;
real<lower=0> sigma.theta;
vector[n_obs] theta;
}
model {
matrix[N, N] L_K;
matrix[N, N] K = cov_exp_quad(x, t1, l1); // Covariance matrix
real sq_sigma = square(sigma); // Unstructured error variance
// Add unstructured error covariance
for (i in 1:N)
K[i, i] = K[i, i] + sq_sigma;
// Add structured error covariance
for (i in 1:n_obs){
int a = (i-1)*len + 1;
int b = i*len;
K[a:b, a:b] = K[a:b, a:b] + cov_exp_quad(x[1:len], x[1:len], t2, l2);
}
L_K = cholesky_decompose(K);
l1 ~ inv_gamma(5, 5);
t1 ~ std_normal();
l2 ~ inv_gamma(5, 5);
t2 ~ std_normal();
sigma ~ std_normal();
tau ~ std_normal();
for(i in 1:n_obs){
theta[i] ~ normal(0,tau);
}
y ~ multi_normal_cholesky(theta[ID], L_K);
}
```