Example using Turing + Stheno

[oizin]

I was trying to put together a simple example using Stheno and Turing. Based on other examples I thought the below would work but I am getting an error about a missing method. Any ideas what is going wrong? I'm using Julia 1.5.2, Stheno v0.6.15 and Turing v0.14.10.

using Turing
using Distributions
using Random
using Stheno

# make some fake data
l = 0.4
σ² = 1.3
σ²_n = 0.05 # noise
x_ = collect(range(-4.0, 4.0; length=100))
k = σ² * stretch(Matern52(), 1 / l)
f = GP(k, GPC())
fx = f(x_,σ²_n) 
y_ = rand(fx)

# prior model
@model gp0(y,x) = begin
    # priors
    σ² ~ LogNormal(0, 1)
    l ~ LogNormal(0, 1)
    σ²_n ~ LogNormal(0, 1)  
    k =  σ² * stretch(Matern52(), 1 ./ l)
    gp = GP(k, GPC())
    # y ~ normal(GP,σ²_n)
    y ~ gp(x,σ²_n + 1e-3)
end

# sample from posterior
m = gp0(y_, x_)
@time chain = sample(m, HMC(0.01, 100), 10)

Resulting error message:

ERROR: MethodError: no method matching logdetcov(::Stheno.FiniteGP{GP{Stheno.ZeroMean{Float64},Stheno.Scaled{Array{Float64,1},Stheno.Stretched{Array{Float64,1},Matern52,typeof(identity)},typeof(identity)}},Array{Float64,1},LinearAlgebra.Diagonal{Float64,FillArrays.Fill{Float64,1,Tuple{Base.OneTo{Int64}}}}})

Thanks for any help - really like the package.

[willtebbutt]

Hmmm looks like I need to update how some of Stheno's internals are implemented. Will open a PR and add this example to our tests.

[andreasko]

Is there any workaround at the moment? Maybe somewhat similar to your example with Soss.jl as described in the docs. No pressure though :)

[willtebbutt]

Hey, sorry for the delay with this. I'm a bit swamped, so I can't promise I'll be able to work on it any time soon :( IIRC there were more things preventing me making this work that I had initially anticipated.

[andreaskoher]

With Stheno 0.7 this issue seems to be solved. For completeness, here is modified example from above:

using Turing
using Distributions
using Random
using Stheno
using AbstractGPs

# make some fake data
l = 0.4
σ² = 1.3
σ²_n = 0.05 # noise
x_ = collect(range(-4.0, 4.0; length=100))
k = σ²*transform(Matern52Kernel(), 1/l)
f = GP(k)
fx = f(x_,σ²_n) 
y_ = rand(fx)

# prior model
@model gp0(y,x) = begin
    # priors
    σ² ~ LogNormal(0, 1)
    l ~ LogNormal(0, 1)
    σ²_n ~ LogNormal(0, 1)  
    k = σ²*transform(Matern52Kernel(), 1/l)
    gp = GP(k)
    y ~ gp(x,σ²_n + 1e-3)
end

# sample from posterior
m = gp0(y_, x_)
chain = sample(m, HMC(0.01, 100), 1000)

using StatsPlots
chain = chain[200:end, :, :]
plot(chain)