Question on Bayesian Inference with Julia, using Turing and/or AdvancedHMC and/or Distributions

I am looking to run some experiments in Julia after growing tired of jumping in and out of Stan and suffering from some performance issues that I think Julia could help eliminate. I am looking to do MCMC to learn parameters or form posterior predictives based on transformations of a standard Bayesian update. I.e. in Stan where one would usually write out:

mu ~ ...
sigma ~ ...
target += normal_lpdf(y | mu, sigma)

I would like to update using some transformation of the lpdf, e.g.

target += w * normal_lpdf(y | mu, sigma) - exp(normal_lpdf(y | mu / 10, sigma ^ 1/2))

Which is not a real example but just to illustrate the kind of thing I am doing. Is this possible in Julia, so far I am having difficulty wrapping my head around how I might achieve this using any of the packages I mentioned in the title. I am much more proficient in Stan, R and Python which is probably the reason for this so I apologise if this question is trivial. But so far all I can see is that in AdvancedHMC I could set my target to be equal to something imported from Distributions.jl and update with that as the target. Is there some way to easily define a custom acceptable target through transformations of log pdfs from Distributions, or perhaps using straight mathematical formulae for the distributions?

Thanks and please let me know if I can provide any more information.

1 Like

In DynamicHMC, you are essentially coding a function that evaluates the log posterior on \mathbb{R}^n so you are free to do any kind of transformation β€” instead of incrementing a variable, I would recommend just returning a sum (AD has an easier time with that usually).

See examples here:

If you need help with a concrete model, please specify it fully, ideally with an MWE that simulates example data.

If by target you mean the log joint probability which is equal to the log posterior up to an offset, then you can write a Julia function for that and use DynamicHMC or AdvancedHMC. Or if you want to use the ~ notation for most of your model, you can use Turing to define your model and use @logpdf() += ... in the model to add a custom quantity to the log joint probability accumulated replacing ... with your custom calculation. Using Turing also means you can use non-HMC and Gibbs inference algorithms and have discrete variables in your model.


Thank you these examples are useful, DynamicHMC is looking promising, will try and potentially get back to you with a more explicit question if needed.

If you want to use Turing, as @mohamed82008 mentioned you can write the model as follows:

using Turing

@model mymodel(Y, w) = begin
    mu ~ Normal()
    sigma ~ truncated(Cauchy(0, 5), 0, Inf)
    for y in Y
        @logpdf() += w * logpdf(Normal(mu, sigma), y) - pdf(Normal(mu/10, sigma^(1/2)), y)

result = sample(mymodel(data, 1.0), NUTS(0.7), 1000)

You can also write the for loop more efficiently using mapreduce, i.e.

@logpdf() += mapreduce(y -> w * logpdf(Normal(mu, sigma), y) - pdf(Normal(mu/10, sigma^(1/2)), y), +, Y)

This sounds very promising, I don’t suppose you have a link to some examples where @logpdf() is defined in a custom way as you describe?

Ah this looks great, will give it a try now, thank you!

1 Like

So I have had a go with this in a couple of ways for a logistic regression type problem:

@model logistic_regression(data, params) = begin
	@unpack y_real, X_real, y_synth, X_synth = data
	@unpack w, Οƒ = params
    coefs ~ MvNormal(zeros(size(X_real)[2]), Diagonal(repeat([Οƒ], size(X_real)[2])))

    for (x, y) in zip(eachrow(X_real), y_real)
    	@logpdf() += logpdf(Bernoulli(logistic(dot(x, coefs))), y)

    @logpdf() += mapreduce(ins -> logpdf(Bernoulli(logistic(ins[1] * coefs)), ins[2]), +, zip(eachrow(X_real), y_real))
    @logpdf() += w * mapreduce(ins -> logpdf(Bernoulli(logistic(ins[1] * coefs)), ins[2]), +, zip(eachrow(X_real), y_real))
    @logpdf() += sum(logpdf.(Bernoulli.(logistic.(X_real * coefs)), y_real))
    @logpdf() += w * sum(logpdf.(Bernoulli.(logistic.(X_synth * coefs)), y_synth))

But run into an error no matter how I do it (I believe the last method is most efficient):

LoadError: ArgumentError: Bernoulli: the condition zero(p) <= p <= one(p) is not satisfied.

Do you know what might be my issue here?

As far as I remember, this arises due to numerical instabilities in the logistic function. We have the following Distribution for this purpose in Turing.

1 Like

Working now, thank you for all the help :slight_smile:

1 Like