# 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?

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.

2 Likes

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

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

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

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)), Diagonal(repeat([σ], size(X_real))))

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 * coefs)), ins), +, zip(eachrow(X_real), y_real))
@logpdf() += w * mapreduce(ins -> logpdf(Bernoulli(logistic(ins * coefs)), ins), +, 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?

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 1 Like