Hi all,

I’m trying to solve an example problem from a paper [1] and wish to do it with high-level modelling. I’m looking into Turing.jl but other packages may exist that I do not know of. For a didactic purpose I want to keep the formulation as close to the paper as possible. From what I could gather, I need to implement this hierarchical model and sample from the posterior:

Using Turing.jl I have the following MWE:

```
using Turing, Distributions, LinearAlgebra
#using Plots
n = 40
h = 1/n
t = h/2:h:1-h/2
x_true = 50*(0.75.*((0.1 .< t) .& (t .< 0.25)) + 0.25*((0.3 .< t) .& (t.<0.32)) + ((0.5 .< t) .& (t .< 1.0)).*sin.(2*pi.*t).^4)
#plot(x_true)
# Construct convolution kernel matrix A:
sig = .03 # Kernel width.
kernel(i,j,h,sig) = (h/(sqrt(pi)*sig^2))*exp(-((i-j)h)^2/(2*sig^2))
A = Matrix{Float64}(undef,n,n)
for i in 1:n
for j in 1:n
A[i,j] = kernel(i,j,h,sig)
end
end
Ax = A*x_true # noise-free observations
err_lev = 2
sigma = err_lev/100 * norm(Ax)/sqrt(n)
eta = sigma * randn(n)
b = Ax + eta # add noise to observations
@model problem(x, A, b) = begin
n = length(b)
α_λ = 1
β_λ = 1e-4
α_δ = 1
β_δ = 1e-4
C = diagm(ones(n))
λ ~ Gamma(n/2+α_λ,0.5*norm(A*x-b)^2+β_λ)
δ ~ Gamma(n/2+α_δ,0.5*x'*C*x+β_δ)
x ~ MvNormal(inv(λ*A'A+δ*C)*λ*A'b,Symmetric(inv(λ*A'A+δ*C)))
end
c = sample(problem(zeros(n), A, b), SMC(1000))
```

However, it seems I’m missing some pieces. I can get Turing to sample, but there’s no `x`

parameter in the output, and also, using other samplers, HMC, for instance I get different error messages, so it seems I’m doing something wrong.

I know this problem can be solved faster/better in other ways (that’s what the paper is about) but want to see how close I can get to the mathematical formulation using high-level modelling.

Comments and suggestions will be appreciated,

thanks!

[1] J. M. Bardsley, “An Efficient MCMC Method for Uncertainty Quantification in Inverse Problems,” pp. 1–13, Nov. 2010. http://hs.umt.edu/math/research/technical-reports/documents/2010/28_BayesianUQforInverseProblems.pdf