Modeling gaussian process with Turing.jl

Hi, I was trying to build a gaussian process model with Turing.jl

I’ve read Gaussianprocesses.jl and Stheno.jl. They must be good.
But for the purpose of learning, I’d like to go through step by step in frame of Turing.jl.

I found this one in search.

Also I’ve looked up this great reference.

So Followings are what I wrote and got an error. :face_exhaling:

begin
	using Turing, Distributions
	using Random, LinearAlgebra
	using Plots, StatsPlots
	using AbstractGPs, KernelFunctions
	using CSV, HTTP, DataFrames
end
url = "https://raw.githubusercontent.com/delphinH/imagetoy/main/dataset/Fish.csv"
fish = CSV.read(HTTP.get(url).body, DataFrame)
begin
	xs = fish[:, :Length1]
	ys = fish[:, :Weight]
end
#Exponentiated quadratic kernel

function EQkernel(x1, x2, l, σ)
    t₁ = sum(x1.^2, dims=2)
    t₂ = sum(x2.^2, dims=2)
    t₃ = 2*x1*x2'
    t = (t₁ .+ t₂') - t₃
    return σ^2 * exp.(-0.5/l^2 * t)
end
EQkernel(xs, xs, 1, 1) #this work!
@model gptest(x, y) = begin
	#Prior
	σ ~ LogNormal(0.0, 0.1)
	l ~ LogNormal(0.0, 1.0)
	σy ~ LogNormal(0.0, 1.0)
	
	#covariance function
	kernel = EQkernel(x, x, σ, l)
    K = kernelmatrix(kernel, x[:, :], obsdims=1)

	#add noise
    y ~ MvNormal(K + LinearAlgebra.I * σy^2)
end
model = gptest(xs, ys)
chain = sample(model, NUTS(), 2000)  # Got an MethodError 

MethodError: no method matching kernelmatrix(::Matrix{Float64}, ::Matrix{Float64}; obsdims=1)

How can I fix this error?
Are there better ways?
I don’t know what kernelmatrix function is.

Thanks for reading and your attention! :pray:

1 Like

Since you said you want to go step by step you’d have to define that function yourself.
This function should return an NxN matrix that contains kernel(X[i], X[j]) at position (i, j), where X is your NxM matrix of independent variables (M = number of independent variables).

Note also that instead of:

you want kernel to be a function that is called inside kernelmatrix for every (i, j), i.e. you’d probably want to instead define:

kernel(x1, x2) = EQkernel(x1, x2, σ, l)
1 Like

Thx for replying!
EQkernel that I made gives me N by N matrix. so I already have one.
:+1:

Based on your advice, I delete K = kernelmatrix(kernel, x[:, :], obsdims=1)
this line.

and now I have another error
MethodError: no method matching +(::Type{KernelFunctions.Kernel}, ::LinearAlgebra.UniformScaling{Float64})

Endless error :joy:

Try to understand the error (it is actually quite informative).
It tells you that in this line:

you’re trying to add something that is of type Type{KernelFunctions.Kernel} to something of type LinearAlgebra.UniformScaling{Float64}.
Is this what you expect? No. You want the first summand to be of type Matrix{Float64}.
In other words: Your K is not a matrix (you said you deleted the line, but you didn’t say what you replaced it with).

My general advice is: Try to understand the errors and go step by step from there (e.g. make sure you understand every line and the output of every line before continuing iterating on the full script/model)

1 Like

I’ve tried to find out that.
It should not be the Type Type{KernelFunctions.Kernel} .
And it wasn’t.
it was a capital letter problem…I should have typed kernel not Kernel

Really sorry for bothering :pleading_face:

No worries at all! :slight_smile:
glad you found it!

2 Likes