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

**URL:** <https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776>\
**Category:** Statistics\
**Tags:** question\
**Created:** [March 9, 2020, 10:40pm UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776 "2020-03-09T22:40:43Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [March 9, 2020, 10:40pm UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/1 "2020-03-09T22:40:43Z")

</div>

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:

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

```

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 10, 2020, 6:59am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/2 "2020-03-10T06:59:19Z")

</div>

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:

> **[GitHub - tpapp/DynamicHMCExamples.jl: Examples for Bayesian inference using...](https://github.com/tpapp/DynamicHMCExamples.jl/)**
>
> Examples for Bayesian inference using DynamicHMC.jl and related packages. - GitHub - tpapp/DynamicHMCExamples.jl: Examples for Bayesian inference using DynamicHMC.jl and related packages.

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

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [March 10, 2020, 10:23am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/3 "2020-03-10T10:23:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [March 10, 2020, 10:32am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/4 "2020-03-10T10:32:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [March 10, 2020, 10:33am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/5 "2020-03-10T10:33:35Z")

</div>

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

```nohighlight
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)

```

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [March 10, 2020, 10:35am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/6 "2020-03-10T10:35:07Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [March 10, 2020, 10:48am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/7 "2020-03-10T10:48:43Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [March 10, 2020, 10:49am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/8 "2020-03-10T10:49:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [March 11, 2020, 11:21am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/9 "2020-03-11T11:21:58Z")

</div>

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

```julia
@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):

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

```

Do you know what might be my issue here?

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [March 11, 2020, 11:30am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/10 "2020-03-11T11:30:13Z")

</div>

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.

> <https://github.com/TuringLang/Turing.jl/blob/3fcf7775d595855dff0a342ae13bc490ec5bfd67/src/stdlib/distributions.jl#L40-L45>

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [March 11, 2020, 11:31am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/11 "2020-03-11T11:31:17Z")

</div>

> <https://github.com/TuringLang/Turing.jl/blob/3fcf7775d595855dff0a342ae13bc490ec5bfd67/test/stdlib/distributions.jl#L11-L16>

---

<div class="post-metadata">

**Author:** ![HarrisonWilde](https://avatars.discourse-cdn.com/v4/letter/h/f07891/32.png) [@HarrisonWilde](https://discourse.julialang.org/u/HarrisonWilde)\
**Post date:** [March 11, 2020, 11:55am UTC](https://discourse.julialang.org/t/question-on-bayesian-inference-with-julia-using-turing-and-or-advancedhmc-and-or-distributions/35776/12 "2020-03-11T11:55:25Z")

</div>

> [@trappmartin](#):
>
> A univariate binomial logit distribution.

Working now, thank you for all the help 🙂