For efficiency take attention to the different way of parameterization:

d = fabs(x[i] - x[j])

sigma = variance

https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function

Also check correctness before use.

```
real Matern12(real d, real rho, real sigma) {
real d_times_rho = d * rho;
return sigma * exp(- d_times_rho);
}
real Matern32(real d, real rho, real sigma) {
real d_times_rho = 1.73205080756887729352744634 * d * rho;
return sigma * exp(- d_times_rho) * (1 + d_times_rho);
}
real sqrt5() {
return 2.2360679774997896964;
}
real Matern52(real d, real rho, real sigma) {
real d_times_rho = d * rho;
return sigma * exp(- sqrt5() * d_times_rho) * (1 + sqrt5() * d_times_rho + 5.0 / 3.0 * square(d_times_rho));
}
```

(It was too late last night, my mistake.) I meant sensor noise

```
vector[N] noise;
real<lower=0> sigma_noise;
noise ~ normal(0, 1); // or student_t
sigma_noise ~ cauchy(0, 1);
y ~ bernoulli_logit(fhat + noise * sigma_noise);
```

Yuedong’s stuff is splines. The answer to what’s better? I don’t know.

Thought again: If using splines the “knots” are fixed a-prior, but if you sample in approximation GP’s

each iteration you’d normally sample **new** knots (anchor points), in MCMC style, eg. in Gibbs sampling.

Ideally the number of “points” M in first first post should vary, to have a trade-off between the frequencies.

As far as I remember correctly, Michael have a paper out, that is not “working well” in HMC.