# Applied Gaussian Processes in Stan, Part 1. A Case Study

**URL:** <https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442>\
**Category:** Modeling\
**Created:** [August 31, 2018, 11:59am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442 "2018-08-31T11:59:35Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![anon79882417](https://avatars.discourse-cdn.com/v4/letter/a/f05b48/32.png) [@anon79882417](https://discourse.mc-stan.org/u/anon79882417)\
**Post date:** [August 31, 2018, 11:59am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/1 "2018-08-31T11:59:35Z")

</div>

Hey all - More specifically, hey GP people -

Writing up some stuff I’ve done over the summer at Aalto (although a bit late…). Aiming for an applied, pedagogical intro to applied Gaussian Processes in Stan, using a variety of open source data sets.

This is Part 1: Intro to GP notation, motivation behind GPs, some properties of kernels, GP regression and ARD.

Some things I could use help with:

1. Hyperparameter priors: I could at least use a length scale that uses the scale of the data, as in Betancourt’s case study.

2. Out of sample: Using Loo package instead of a naive RMSE.

3. I was trying to show the smoothness changes in different cases of Matern kernel, but the sample size to high/vizualizations too small to show the “wiggliness” of the kernels.

Part 2 is coming. It has some different models including GP classificaiton, GP time series (modeling directly and using a GP prior in GARCH), and then a survival model. I split it up because the case study was becoming too long.

[applied\_gaussian\_processes\_in\_stan.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/291cd8cbd6568dd32512b336de7d7f0f4482e615.pdf) (419.7 KB)

---

<div class="post-metadata">

**Author:** ![jjramsey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jjramsey/32/3319_2.png) [@jjramsey](https://discourse.mc-stan.org/u/jjramsey)\
**Post date:** [August 31, 2018, 12:51pm UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/2 "2018-08-31T12:51:23Z")

</div>

I noticed that you are using some functions, such as `add_diag` and `gp_exp_quad_cov`, that do not appear to be defined in the `functions` block and also do not appear in the documentation for Stan. Are these undocumented features in the latest CmdStan release, or are they from a development version?

---

<div class="post-metadata">

**Author:** ![anon79882417](https://avatars.discourse-cdn.com/v4/letter/a/f05b48/32.png) [@anon79882417](https://discourse.mc-stan.org/u/anon79882417)\
**Post date:** [September 1, 2018, 10:42pm UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/3 "2018-09-01T22:42:33Z")

</div>

good question.

I think `add_daig` is in 2.18 math library, but it hasn’t been added to the language yet. There are covariance functions such as the matern32 and matern52 that I’m using locally, and still pushing through PRs to get them merged into the Stan library. Good call on the docs, I need to update the stan docs to show that at least `add_diag` has been included. Besides these technical issues, do you have any comments about modeling issues or the pedagogy of the doc in general? I know it’s basic, but I’m trying to show users how to build up from nothing!

---

<div class="post-metadata">

**Author:** ![jjramsey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jjramsey/32/3319_2.png) [@jjramsey](https://discourse.mc-stan.org/u/jjramsey)\
**Post date:** [September 2, 2018, 12:51pm UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/4 "2018-09-02T12:51:26Z")

</div>

I am curious about how identifiable the length scales in ARD are. I’ve had issues with that myself in dealing with GPs, and I’m curious as to how much of that is at my end and how much is just inherent to GP regression.

---

<div class="post-metadata">

**Author:** ![anon79882417](https://avatars.discourse-cdn.com/v4/letter/a/f05b48/32.png) [@anon79882417](https://discourse.mc-stan.org/u/anon79882417)\
**Post date:** [September 2, 2018, 4:48pm UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/5 "2018-09-02T16:48:25Z")

</div>

yeah, that warrants another case study. I could trivially make the `real<lower=0>length_scale[D]` an ordered vector, but then we would be assuming that the length\_scale1 \< length\_scale2 \< length\_scale3 \< … and I don’t think this is a good assumption. I could run a model a bunch of times w/ different random seeds and see where they tend to fall, and re-order the dataset in that way but it’s a hack… I’ll put some actually though into it, thank you.

---

<div class="post-metadata">

**Author:** ![increasechief](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/increasechief/32/2039_2.png) [@increasechief](https://discourse.mc-stan.org/u/increasechief)\
**Post date:** [September 2, 2018, 6:03pm UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/6 "2018-09-02T18:03:46Z")

</div>

Really fantastic work @anon79882417 I took in a lot of lessons from your methodology and presentation.

---

<div class="post-metadata">

**Author:** ![jgcolman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jgcolman/32/836_2.png) [@jgcolman](https://discourse.mc-stan.org/u/jgcolman)\
**Post date:** [December 11, 2018, 10:53am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/7 "2018-12-11T10:53:31Z")

</div>

That is a good question. I find there is often an interplay between the length scales and the variance, so that it may not make sense to treat the two as independent hyper-parameters.

---

<div class="post-metadata">

**Author:** ![NTorabi](https://avatars.discourse-cdn.com/v4/letter/n/b487fb/32.png) [@NTorabi](https://discourse.mc-stan.org/u/NTorabi)\
**Post date:** [July 15, 2019, 11:48pm UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/8 "2019-07-15T23:48:41Z")

</div>

Hey @anon79882417, thanks for putting all these together. By any chance, have you prepared the part related to classification? I am trying to understand the function gp\_pred\_rng(…) and I thought asking your help/write-up. Thanks

---

<div class="post-metadata">

**Author:** ![NTorabi](https://avatars.discourse-cdn.com/v4/letter/n/b487fb/32.png) [@NTorabi](https://discourse.mc-stan.org/u/NTorabi)\
**Post date:** [July 16, 2019, 3:06am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/10 "2019-07-16T03:06:31Z")

</div>

@anon79882417 Yes I am a real person with a real question. Considering that the likelihood is not Guassian in the GP classification, I am trying to understand how did you come up with the posterior predictive for the test data in function gp\_pred\_rng(…)? I am basically referring to [this post](https://github.com/stan-dev/docs/issues/16) that you have for classification.

---

<div class="post-metadata">

**Author:** ![anon79882417](https://avatars.discourse-cdn.com/v4/letter/a/f05b48/32.png) [@anon79882417](https://discourse.mc-stan.org/u/anon79882417)\
**Post date:** [July 16, 2019, 8:44am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/11 "2019-07-16T08:44:10Z")

</div>

The integral were evaluation for a non Gaussian likelihood is this:

\int \pi(y|f, X, \theta) \pi(f | X, \theta) \pi(x|\theta)

In Gaussian likelihood case the \pi(f | •) gets marginalized out, but here it does not. Instead you generate latent f from:

1. Computer covariance function K
2. Compute cholesky decomp K=LL’
3. Then use L to sample from MV Gaussian and f is a realization from that MV Gaussian so: f = \rho \* L, \rho \sim N(0, 1)
4. And then use that f in your likelihood which I depends on the model you specify. f could be a prior for a parameter in a linear model, you could use it directly in the logistic or multi nominal likelihood, etc.

A good reference is GPML, and I think most of this is in the Classification chapter except step 4.

The equations for ‘gp\_pred\_rng’ are in the document, I believe on the first page. Keywords are “posterior predictive”

---

<div class="post-metadata">

**Author:** ![anon79882417](https://avatars.discourse-cdn.com/v4/letter/a/f05b48/32.png) [@anon79882417](https://discourse.mc-stan.org/u/anon79882417)\
**Post date:** [July 16, 2019, 8:47am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/12 "2019-07-16T08:47:34Z")

</div>

> [@NTorabi](#):
>
> Yes I am a real person with a real question.

Yes, apologies. Please let me know if you have any other questions, and feel free to open another thread with code if you have questions on implementation details. There are better references than that document.

Andre

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [July 16, 2019, 9:02am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/13 "2019-07-16T09:02:44Z")

</div>

See Section 3.3 in [https://betanalpha.github.io/assets/case\_studies/gp\_part1/part1.html#33\_simulating\_from\_a\_gaussian\_process\_conditional\_on\_non-gaussian\_observations](https://betanalpha.github.io/assets/case_studies/gp_part1/part1.html#33_simulating_from_a_gaussian_process_conditional_on_non-gaussian_observations) for an example of how to implement a non-Gaussian observational model. If you replace the Poisson distributions (for integer-valued observations) with Bernoulli distributions (for binary-valued observations) then you get a Gaussian process model amenable to classification. In other words you end up modeling the logit probability with the Gaussian process.

---

<div class="post-metadata">

**Author:** ![Max\_Mantei](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/max_mantei/32/7912_2.png) [@Max\_Mantei](https://discourse.mc-stan.org/u/Max_Mantei)\
**Post date:** [November 21, 2019, 10:57am UTC](https://discourse.mc-stan.org/t/applied-gaussian-processes-in-stan-part-1-a-case-study/5442/14 "2019-11-21T10:57:30Z")

</div>

Hi! Cool stuff, thanks! Did you ever find the time to write that second part? :)