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

Specific Domains Statistics
1

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
2

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.

3

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.

2 Likes
4

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.

5

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)
    end
end

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

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)
3 Likes
7

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?

8

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

1 Like
9

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)
    end

    @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))
end

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?

10

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
11
12

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

1 Like