Gaussian process regression with heteroscedasticity - Initialisation Issues

[Benjamin Frot]

Hi Stan users,

I have been using (R)STAN  for a few weeks now and I find it great. 

Recently, I have been trying to implement a simple model of Bayesian Linear regression that corrects for heteroscedasticity.

The idea is to put a Gaussian Process prior on the variance like so:

 y_i =  a + bX_i + e_i with e_i ~ N(0, r(x) = e^g(x))  and g(x) ~ GP(mu0,k(x,x') ).

k(x,x') is the standard squareq-exponential covariance function.

Attached is my attempt at implementing this model in Stan. Here, g should be a constant function. 

But it keeps failing at initialisation, even when I provide my own initial parameters.

I think my problem stems from a misunderstanding of the way STAN works... or my code is simply not implementing what I think it is.

I would be happy if anyone could tell me whether something is fundamentally wrong with my STAN code.

Thanks!

Ben