Thank you both for the quick replies and very helpful resources.
I tried to amend the model I had:
matrix[N, N] cov = cov_exp_quad(x, alpha, rho) + diag_matrix(rep_vector(sigma, N));
matrix[N, N] L_cov = cholesky_decompose(cov);
vector[N] f = L_cov * y_tilde;
y ~ neg_binomial_2_log(f, phi);
As suggested by @bbbales2 I’ll start with a GP on the mean and estimate sigma, which I think is all I need. The problem is that the model is not yet correct (my little experience with Stan doesn’t help).
I also would like to fit this model to two different y’s and find where these two GPs differ significantly (I’ll also need to contain a specific offset to each y)
This sounds like a generated quantities thing. Fit two GPs to data, and then in generated quantities draw from them and compare the numbers you get out? Is this a model comparison thing? Or is it figuring out the differences in two separate but comparable processes in your application? Doing model comparison stuff is always trickier than it seems is why I’m asking (so I don’t point you in the wrong direction).
A description of my bioinformatics problem might help to clarify. Briefly, I have a phenotypic response, y var, to multiple individual mutations located across the chromosome, x var. Now this was acquired for multiple cell lines and a control condition. My idea is to fit a GP in the chromosome for each condition and compare them to the GP fitted for the control. Basically, I would like to find the x’s for which the difference between the conditions and control is greatest.
I presume I don’t necessarily need to do a complete model comparison, but would be rather drawing from the GPs and compare.