GP model

[mvsoom]

Hi everyone,

I am unable to get a simple Gaussian process (GP) model to work and need help as I can't seem to figure it out. The model consists of observations of a GP gp at times t with known observation noise σ².

using RxInfer, AbstractGPs

gp = GP(with_lengthscale(Matern32Kernel(), 5.))
t = 1.:100.

@model function gp_model(t, gp, σ²::ConstVariable)
  n = length(t)
  u = randomvar(n)
  y = datavar(Float64, n)
  
  gpt = gp(t) # Condition the GP and convert it to a MvNormal
  
  u .~ MvNormal(μ = mean(gpt), Σ = cov(gpt))
  y .~ NormalMeanVariance(u, σ²)
end

result = infer(
  model = gp_model(t, gp, noise_variance),
  data  = (y = observations, ),
  #constraints = FullFactorisation(),
  #initmarginals = (u = vague(NormalMeanPrecision),),
  #initmessages = (u = vague(NormalMeanPrecision),),
  free_energy   = true,
  iterations    = 10, 
  showprogress  = true
)

This yields:

ProductOf object cannot be used as a functional form in inference backend. Use form constraints to restrict the functional form of marginal posteriors.

I tried all form constraints I could find in the docs, but I suspect I am overlooking something basic here. Any help would be much appreciated.

---

The code above is a summary, the Pluto notebook can be found here: gp.jl.txt. I am on latest RxInfer with Julia version 1.10.1.

[albertpod]

Hi @mvsoom!

I am not sure you can use broadcasting in RxInfer.jl this way.
What do you mean by u .~ MvNormal()?

Can't we say something like this instead:

using RxInfer, AbstractGPs

gp = GP(with_lengthscale(Matern32Kernel(), 5.))
t = 1.:100.
@model function gp_model(t, gp, σ²::ConstVariable)
    y = datavar(Vector{Float64})
    
    gpt = gp(t) # Condition the GP and convert it to a MvNormal

    u ~ MvNormal(μ = collect(mean(gpt)), Σ = cov(gpt))
    y ~ MvNormal(μ=u, Σ=σ²*diageye(length(t)))
end

result = infer(
  model = gp_model(t, gp, 1.0),
  data  = (y = randn(length(t)), ),
  free_energy   = true,
  showprogress  = true
)

So now, we will use a single variable, u, as a prior for our observations.
As for the observations, I changed it to multivariate distribution with diagonal covariance.

Note that I needed to use collect when calling mean(gpt) as it doesn't return a vector.

[mvsoom]

Excellent, thank you! I got confused with the difference between a single randomvar() that is vector-valued vs randomvar(n). I have another question but I will post that in a different discussion.