INLA.jl implementation help

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.

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

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,


1 Like

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

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 :slight_smile:

1 Like

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,


1 Like

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

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

1 Like