# Gaussian Process Model with Turing

**URL:** <https://discourse.julialang.org/t/gaussian-process-model-with-turing/42453>\
**Category:** Probabilistic Programming\
**Tags:** regression, turing, gaussian-process\
**Created:** [July 2, 2020, 9:10pm UTC](https://discourse.julialang.org/t/gaussian-process-model-with-turing/42453 "2020-07-02T21:10:36Z")\
**Posts on this page:** 1\
**Showing post:** 23

<div class="post-metadata">

**Author:** ![luiarthur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luiarthur/32/16405_2.png) [@luiarthur](https://discourse.julialang.org/u/luiarthur)\
**Post date:** [August 31, 2020, 10:24pm UTC](https://discourse.julialang.org/t/gaussian-process-model-with-turing/42453/23 "2020-08-31T22:24:42Z")

</div>

This runs, but I haven’t looked at the results closely.

Note that I used `KernelFunctions.jl`, which implements the squared exponential covariance function previously written as `sqexp_cov_fn`.

Also cleaned up a few things. You’ll probably want to think about the priors for the application.

```julia
using CSV
using DataFrames
using Turing
using Distributions
using LinearAlgebra
using KernelFunctions

df = DataFrame(
    subject = repeat(1:5, inner=3),
    obese = repeat(rand(Bool, 5), inner=3),
    timepoint = [1,2,3,1,3,4,1,2,5,1,4,5,1,3,5],
    bug = rand(Beta(0.9, 5), 15),
    nutrient = rand(Beta(0.9,5), 15)
)

sekernel(alpha, rho) = 
  alpha^2 * KernelFunctions.transform(SEKernel(), sqrt(0.5)/rho)

@model function myGP(y, Z, X, jitter=1e-6)
    # Dimensions of GP predictors. Note that if X is a single column 
    # it needs to have dimension Nx1 (a matrix), not an array of 
    # length N, so that it can be processed properly by the distance function.
    N, P = size(X)
    
    # Dimensions of linear model predictors
    J = size(Z, 2) # Z should be N x J
    
    # Priors.
    mu ~ Normal(0, 1)
    sig2 ~ LogNormal(0, 1)
    alpha ~ LogNormal(0, 0.1)
    rho ~ LogNormal(0, 1)

    # Prior for linear model coefficients
    beta ~ filldist(Normal(0, 1), J)
    
    # GP Covariance matrix
    kernel = sekernel(alpha, rho) # covariance function
    K = kernelmatrix(kernel, X, obsdim=1) # cov matrix
    K += LinearAlgebra.I * (sig2 + jitter)

    # Sampling Distribution.
    y ~ MvNormal(mu .+ Z * beta, K)
end

gp = myGP(df.bug,
          Matrix(df[!,[:subject, :obese, :nutrient]]),
          Matrix{Float64}(df[!,[:timepoint]]))

@time chain = sample(gp, HMC(0.01, 100), 200)

```

---

_[View the full topic](https://discourse.julialang.org/t/gaussian-process-model-with-turing/42453)._