I tried to implement an algorithm which generates the target values with a bit of noise in R. (based on an example in this paper " Monte Carlo Implementation of Gaussian Process Models for Bayesian Regression and Classification",Neil,1997.

Target Distribution with adjusted deviations

f.x = function(x){

s = sample(c(1, 2), size = 1, replace = TRUE, prob = c(0.95, 0.05))

ifelse(s==1,0.3+0.4 *x+0.5* sin(2.7 *x)+1.1/(1+x^2)+rnorm(1,0,0.1),*

0.3+0.4 x+0.5 *sin(2.7* x)+1.1/(1+x^2)+rnorm(1,0,1))

}

f.x=Vectorize(f.x)

#Generate train and test set

#set.seed(123)

x = seq(-3,3,length=200)

targets = f.x(x)

dataset = data.frame(cbind(x,targets))

train = (data.frame(cbind(train=dataset$x[1:100],targets=dataset$targets[1:100])))

test = (data.frame(cbind(test=dataset$x[101:200],targets=dataset$targets[101:200])))

I assumed Student distributed noise and my model looks like this:

functions{

// radial basis, or square exponential, function

matrix rbf(int Nx, vector x, int Ny, vector y, real eta, real rho){

matrix[Nx,Ny] K;

for(i in 1:Nx){

for(j in 1:Ny){

K[i,j] = eta*exp(-rho*pow(x[i]-y[j], 2));

}

}

return K;

}

```
vector post_pred_rng(real dof,real s,real eta, real rho, real sn, int No, vector xo, int Np, vector xp, vector yobs){
matrix[No,No] Ko;
matrix[Np,Np] Kp;
matrix[No,Np] Kop;
matrix[Np,No] Ko_inv_t;
vector[Np] mu_p;
matrix[Np,Np] Tau;
matrix[Np,Np] L2;
vector[Np] f_predic;
//vector[Np] Yp;
// use of matrix multiplication for the sn noise
Ko = rbf(No, xo, No, xo, eta, rho) + diag_matrix(rep_vector(1, No))*sn ;
Kp = rbf(Np, xp, Np, xp, eta, rho) ;
Kop = rbf(No, xo, Np, xp, eta, rho) ;
Ko_inv_t = Kop' / Ko;
mu_p = Ko_inv_t * yobs;
Tau = Kp - Ko_inv_t * Kop;
L2 = cholesky_decompose(Tau);
f_predic = mu_p + L2*rep_vector(normal_rng(0,1), Np);
return f_predic;
}
```

}

data{

int<lower=1> N1;

int<lower=1> N2;

vector[N1] X;

vector[N1] Y;

vector[N2] Xp ;

int<lower=1, upper=N2> observed_idx[N1];

}

transformed data{

vector[N1] mu;

for(n in 1:N1) mu[n] = 0;

}

parameters{

real<lower=0> eta_sq;

real<lower=0> inv_rho_sq;

real<lower=5> dof; // Degrees of freedom for t

real<lower=0> s; // scale parametr for T

vector[N2] f_tilde;

}

transformed parameters{

real<lower=0> rho_sq;

rho_sq = inv(inv_rho_sq);

}

model{

matrix[N1,N1] Sigma;

matrix[N1,N1] L_S;

vector[N1] f_predict;

// chol dec

f_tilde ~ normal(0, 1);

Sigma = rbf(N1, X, N1, X, eta_sq, rho_sq) + diag_matrix(rep_vector(1, N1))*s;

L_S = cholesky_decompose(Sigma);

f_predict = L_S * f_tilde;

s ~ student_t(4,0,1);

dof ~ normal(4.5,0.1);

Y ~ student_t(dof,f_predict, s);

eta_sq ~ normal(0,1);

inv_rho_sq ~ cauchy(0,5);

}

generated quantities{

vector[N2] f_pred;

vector[N2] Ypred;

f_pred = post_pred_rng(dof,s,eta_sq, rho_sq, s,N1, X, N2, Xp, Y);

for (n in 1:N2)

Ypred[n] = student_t_rng(dof, f_pred[n], s);

}

this is the resulted plot for student noise