Understanding the correct syntax for fitting models with Turing.jl

Specific Domains Probabilistic programming
, ,
1

Hi everyone,

I’m trying to learn how to use Turing.jl to fit models to observational data, but I find the syntax somewhat confusing, and the documentation isn’t very beginner-friendly when it comes to real-world examples.

As a test case, I attempted to fit a 2D Sérsic profile to simulated data with Gaussian noise. Below is my current implementation:

julia

using Turing
using Distributions
using Plots
using StatsPlots
using PairPlots
using CairoMakie

# Approximate calculation of b_n
b_n(n) = 2n - 0.324

# 2D Sérsic function
function sersic2d(x, y, I_e, r_e, n, x0, y0)
    r = sqrt.((x .- x0) .^ 2 .+ (y .- y0) .^ 2)
    return I_e .* exp.(-b_n(n) .* ((r ./ r_e) .^ (1 / n) .- 1))
end

# True parameters
I_e_true = 100.0
r_e_true = 5.0
n_true = 2.5
x0_true = 0.0
y0_true = 0.0

# Grid of observational points
x = range(-10, 10, length=100)
y = range(-10, 10, length=100)
X, Y = [repeat(x, 1, length(y));], [repeat(y', length(x), 1);]

# Simulate data with noise
σ_noise = 5.0
intensity_true = sersic2d(X, Y, I_e_true, r_e_true, n_true, x0_true, y0_true)
data_obs = intensity_true .+ rand(Normal(0, σ_noise), size(intensity_true))

# Plot simulated data
Plots.heatmap(x, y, data_obs, title="Simulated Data (Sérsic with Noise)", xlabel="x", ylabel="y")

# Define Turing model
@model function fit_sersic2d(x, y, data_obs)
    # Priors for parameters
    I_e ~ truncated(Normal(50, 50), 0, Inf)
    r_e ~ truncated(Normal(3, 5), 0, Inf)
    n ~ truncated(Normal(2, 2), 0.5, 10)
    x0 ~ Normal(0, 5)
    y0 ~ Normal(0, 5)
    σ_noise ~ truncated(Normal(0, 10), 0, Inf)

    # Sérsic model
    model_intensity = sersic2d(x, y, I_e, r_e, n, x0, y0)

    # Likelihood
    for i in eachindex(data_obs)
        data_obs[i] ~ Normal(model_intensity[i], σ_noise)
    end
end

# Perform sampling
chain = sample(fit_sersic2d(X, Y, data_obs), NUTS(500,0.65), 1000)

# Pair plot
fig = PairPlots.pairplot(chain)

My main questions:

  1. General syntax for fitting a model: If I have an arbitrary model and a dataset, what is the correct way to define the Bayesian inference problem in Turing? Are there best practices for structuring models efficiently?
  2. Handling noise correctly: In my example, I assume Gaussian noise with a prior on its standard deviation. Is this the correct approach, or is there a better way to handle uncertainties in the observations?
  3. Custom log-likelihood functions: Instead of defining likelihoods through the data_obs[i] ~ Normal(...) syntax, is it possible (and advisable) to define a custom log-likelihood function? If so, how can it be integrated into a Turing model?

In general, I’d like help in understanding the proper approach to fitting custom models that I define myself. The documentation doesn’t seem to clearly explain how to do this in a structured way.

I’d really appreciate any insights from those experienced with Turing.jl! Any examples or references to relevant documentation would also be very helpful.

2

Hey! Would you elaborate on what you find confusing? Your code looks to me like you understood the idea. It’s possible to define custom distributions and to add arbitrary log-probabilities. In terms of modelling things, it’s all about your specific data and what kind of assumptions you want to make. Your noise looks quite standard.

3

I think this is a great resource:

It may be of limited value if you don’t have the book, but if you already have a solid understanding of Bayesian analysis/inference, you might be able to get quite a bit out of it.

4

Sorry for the late reply! My question is more about the core concepts behind model definition in Turing. I come from a background where I’ve always used Python packages like emcee or dynesty to fit my models, and I had to explicitly write my own log-likelihood functions for the fitting process.

Turing’s syntax feels quite abstract to me. For example, in the model I wrote, should the observed data be passed as arguments to the function? What log-likelihood does Turing use by default for sampling?

Additionally, if I wanted to define my own log-likelihood, how should I go about it? For instance, if I have data with associated uncertainties, how do I properly incorporate these into the likelihood function?

Or even, how do I start a chain from an initial position in parameter space?

1 Like
5

Right, cool! You’re familiar with the math so it shouldn’t take long to get accustomed to this flavor Bayesian modelling. There may be various differences that I don’t understand and I’m just a fellow learner but I’ll add a few things. By the way, I recently read this new paper by the Turing developers (not yet properly published I guess but readable) and thought it was illuminating.

You can pass observations as arguments (old-style) or use the conditioning style (newer) when the variable is left of ~, that is, use model | (data_obs = data_obs, ) if you have in the model e.g. data_obs ~ MvNormal(predictions, σ^2*I). Everything except these tilde statements are normal Julia. If I understand correctly, here just the multivariate normal probability density function is used to add log-probability to the calculation.

The basic idea is that you write a generative model, i.e., something that looks like a function that you would write to simulate that data. Sorry I’m not sure what you mean by data and associated uncertainties. Maybe more examples would help.

For the initial position things, initial_params keyword argument seems to be often accepted in the sample method as shown here. (Turing doesn’t seem super mature and well-documented yet but definitely useable.)

1 Like