INLA.jl implementation help

bparbhu · May 16, 2023, 10:18pm

So I started trying to reverse engineer this comparison done in R that performs a very close implementation of INLA from scratch. Though, I’m having difficulty reproducing the results of the “from-scratch” implementation due to differences in how Julia can perform computations dealing with infinite values versus R.

Here’s the link to my example.

github.com

bparbhu/INLA.jl/blob/master/src/inla_example.jl

using DataFrames, DataFramesMeta, Gadfly, Random, Distributions, StatsBase, SparseArrays, Optim, RCall, Compose, LinearAlgebra, CategoricalArrays


Random.seed!(1234)
n = 100
alpha_true = 0.6
sigma_true = 0.1
beta_true = 10

# simulate x_t as AR1, p_t as logit(x_t), and y_t as Bernoulli(p_t)
function simulate_ar1(n, alpha, sigma)
    x = Vector{Float64}(undef, n)
    x[1] = randn() * sigma
    for t in 2:n
        x[t] = alpha * x[t-1] + randn() * sigma
    end
    return x
end

x_true = simulate_ar1(n, alpha_true, sigma_true)

This file has been truncated. show original

What I’m trying to reproduce is essentially this example done in R

https://stefansiegert.net/inla-project/inla-from-scratch

Where I am trying to just bring in the res object into Julia and then use it for comparison against the julia implementation.

I could definitely use a second set of eyes with this but if we can reproduce the example done in R. I think we can have a good standing point for what to implement for an INLA.jl package. We just need to get the Julia version scratch going and then we can build from there.

Also happy to have some collaborators on this.

Thanks again,

-Brian

bparbhu · May 16, 2023, 10:48pm

The comparison I have is really bad and I’m trying to mimic the distribution that is done in the tutorial

juliohm · May 17, 2023, 10:53am

Very interesting tutorial on INLA @bparbhu , thanks for sharing. I didn’t have time to read it carefully, but INLA is something that we are pursuing in our GeoStats.jl stack for spatial processes. We have implemented the Laplacian operators for arbitrary meshes already and the underlying Gaussian solvers for high-dimensional simulation. Any help is appreciated

bparbhu · May 17, 2023, 6:42pm

Hey @juliohm ,

Thanks again for the feedback I appreciate it! I didn’t know that JuliaEarth was working on an INLA solution and its very cool that you guys were able to replicate the Laplacian operators, that’s awesome. If that’s the case. Would you want to move this under JuliaEarth or do you have a branch in GeoStats.jl that we can work this out. Also, curious just about identifying differences in implementation through this kind of exercise as well. So if you have an implementation this can also serve as a test to see how close it is to R-INLA. So happy to help, just want to finish this up for my own sake.

Thanks again,

-Brian

juliohm · May 17, 2023, 7:23pm

Check the SPDEGS solver. It implements Gaussian simulation via Laplace operators (the latent process used in INLA):

https://juliaearth.github.io/GeoStats.jl/stable/solvers/builtin.html#GeoStatsSolvers.SPDEGS

Here is an example simulation on the sphere:

using GeoStats
using GeoStatsViz
import GLMakie as Mke

mesh = simplexify(Sphere((0,0,0), 1))

# ask for 3 realizations of the process over the sphere
problem = SimulationProblem(mesh, :Z=>Float64, 3)

# latent process of INLA
solver = SPDEGS()

# perform simulation
ensemble = solve(problem, solver)

# view first realization
viewer(ensemble[1])