Hi,

I’m working on a problem involving multiple, independent gaussian processes that I would like to fit jointly together with a common mean function. Specifically, I would like to fit a model of the form

mu ~ GP(0, C_mu)

f_i ~ GP(0, C_i), independent for i = 1, …, n

Y_ij | mu, f_i ~ N(beta0 + mu(t_ij) + f_i(t_ij), sigma^2)

but with the constraint that sum_{i=1}^n f_i(t) = 0 so that each GP f_i represents a subject-specific deviations from the common GP mean function mu.

I thought that I could do something like

`y ~ gp(t, k = 10, c = 5/4) + gp(t, k = 10, c = 5/4, by = id)`

since the gp-terms in brms are zero mean and the intercept is included by default, but that fails miserably. Sometimes it won’t even go into warmup (failing after a lot of “Gradient evaluated at the initial value is not finite”), and other times it starts but with very bad sampling behavior.

Is there something I’m missing in the parametrization or something specific related to the default priors?

Thanks in advance!

Here’s an example

```
library(brms)
library(mvtnorm)
seCov <- function(s, t, alpha, rho) {
alpha^2 * exp(-(s-t)^2 / (2 * rho^2))
}
mu <- function(t) {
sin(t*10) + t*10
}
simGP <- function(t, alpha = 1, rho = 0.1, sigma = 0.1) {
cMat <- outer(t, t, seCov, alpha = alpha, rho = rho)
f <- as.numeric(rmvnorm(1, sigma = cMat))
y <- mu(t) + f + rnorm(length(t), 0, sigma)
list(f = f, y = y)
}
n <- 25 #25 GPs
t <- seq(0, 1, length.out = 20) #20 obs for each
set.seed(12345)
d <- do.call("rbind", lapply(1:n, function(id) {
gp <- simGP(t)
data.frame(id = id, t = t, f = gp$f, y = gp$y)
}))
d$id <- as.factor(d$id)
m <- brm(y ~ gp(t, k = 10, c = 5/4) + gp(t, k = 10, c = 5/4, by = id), data = d)
m
```