# Bayesian logistic regression with Turing.jl

**URL:** <https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105>\
**Category:** Probabilistic Programming\
**Tags:** turing, monte-carlo\
**Created:** [April 27, 2021, 12:55pm UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105 "2021-04-27T12:55:42Z")\
**Posts on this page:** 10\
**Page:** 2

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [May 9, 2021, 4:21pm UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/21 "2021-05-09T16:21:21Z")

</div>

For example, yes:) Tullio.jl that I mentioned above often ends up using LoopVectorization.jl under the hood. Though using `@addlogprob!` _within_ the for-loop probably won’t buy you much, i.e. store the results in some other variable and the sum/reduce this in a single `@addlogprob!` instead. With Tullio.jl you can also perform efficient reductions + it will automatically derive adjoint-rules for Zygote.jl and Tracker.jl.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [May 10, 2021, 10:14am UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/22 "2021-05-10T10:14:42Z")

</div>

I’m seeing a slow down using Tullio for summing log probabilities. Probably I’m just doing it wrong.

```julia
using BenchmarkTools, LoopVectorization, Tullio

X = rand(960*950)
Y = rand(960*950)
sig = 1

julia> @benchmark log_lik1 = sum(map(StatsFuns.normlogpdf, X, sig, Y))
BenchmarkTools.Trial: 
  memory estimate: 288 bytes
  allocs estimate: 6
  --------------
  minimum time: 696.503 ns (0.00% GC)
  median time: 903.497 ns (0.00% GC)
  mean time: 1.075 μs (3.42% GC)
  maximum time: 203.784 μs (99.31% GC)
  --------------
  samples: 10000
  evals/sample: 143

julia> @benchmark @tullio log_lik2 := StatsFuns.normlogpdf(X[i],sig,Y[i])
BenchmarkTools.Trial: 
  memory estimate: 6.80 KiB
  allocs estimate: 158
  --------------
  minimum time: 1.191 ms (0.00% GC)
  median time: 1.502 ms (0.00% GC)
  mean time: 1.693 ms (0.00% GC)
  maximum time: 4.586 ms (0.00% GC)
  --------------
  samples: 2914
  evals/sample: 1

```

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [May 10, 2021, 10:53am UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/23 "2021-05-10T10:53:19Z")

</div>

A couple of notes here:

1. Your `map` is actually only computing a _single_ element:) That’s why it’s significantly faster. `map` and `broadcast` have different behavior when the arguments are of different shapes: `map` I believe falls back to the smallest common size, e.g. in this case it’s 1, so it will only call the function _once_, while `broadcast` will do what you want.
2. When using BenchmarkTools.jl, you need to use interpolation to ensure that you’re not seeing performance decrease because of scoping issues and whatnot.
3. If-statements within the function you’re broadcasting over really hampers performance. In the case of `StatsFuns.normlogdf` there are two versions of the method: 1) takes mean, std, and value, from which it computes the z-value and then finally the logpdf, 2) takes the z-value directly and computes logpdf. In (1) there are some if-statements while in (2) there are none. As a result, (2) will be muuuuch faster when broadcasted/vectorized.
4. Tullio.jl can’t derive the gradients automatically when the RHS include “complex” function calls, e.g. `StatsFuns.normlogpdf`. So in this particular case you probably shouldn’t use it.

All in all, your example should be something like:

```julia
using Tullio, CUDA

X = rand(960*950);
Y = rand(960*950);
sig = 1.0

f(X, Y, sig) = sum(StatsFuns.normlogpdf.(X, sig, Y))
g(X, Y, sig) = @tullio lp = StatsFuns.normlogpdf(X[i], sig, Y[i])

@benchmark $f($X, $Y, $sig)
# BenchmarkTools.Trial: 
# memory estimate: 6.96 MiB
# allocs estimate: 2
# --------------
# minimum time: 6.359 ms (0.00% GC)
# median time: 6.662 ms (0.00% GC)
# mean time: 6.847 ms (1.98% GC)
# maximum time: 11.060 ms (23.53% GC)
# --------------
# samples: 730
# evals/sample: 1

@benchmark $g($X, $Y, $sig)
# BenchmarkTools.Trial: 
# memory estimate: 128 bytes
# allocs estimate: 5
# --------------
# minimum time: 6.571 ms (0.00% GC)
# median time: 6.799 ms (0.00% GC)
# mean time: 6.821 ms (0.00% GC)
# maximum time: 8.835 ms (0.00% GC)
# --------------
# samples: 732
# evals/sample: 1

# Circumvent if-statements used in `StatsFuns.normlogpdf(::Number, ::Number, ::Number)`.
h(X, Y, sig) = sum(StatsFuns.normlogpdf.(Y - X / sig))
i(X, Y, sig) = @tullio (+) lp = StatsFuns.normlogpdf(Y[i] - X[i] / sig)

@benchmark $h($X, $Y, $sig)
# markTools.Trial: 
# memory estimate: 20.87 MiB
# allocs estimate: 6
# --------------
# minimum time: 2.947 ms (0.00% GC)
# median time: 3.497 ms (0.00% GC)
# mean time: 4.436 ms (12.28% GC)
# maximum time: 13.297 ms (25.54% GC)
# --------------
# samples: 1123
# evals/sample: 1

@benchmark $i($X, $Y, $sig)
# BenchmarkTools.Trial: 
# memory estimate: 128 bytes
# allocs estimate: 5
# --------------
# minimum time: 476.989 μs (0.00% GC)
# median time: 556.951 μs (0.00% GC)
# mean time: 562.737 μs (0.00% GC)
# maximum time: 999.113 μs (0.00% GC)
# --------------
# samples: 8765
# evals/sample: 1

```

But as mentioned above, you want gradients through it, so your best bet here is probably just to use `h` in the above example. EDIT: Actually, it might work very nicely with ForwardDiff.jl? I just know that for reverse-rules then you can’t use “complex” function calls on the RHS.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [May 10, 2021, 11:07am UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/24 "2021-05-10T11:07:58Z")

</div>

Thanks a lot! I knew I must have been doing something silly.  
I am using ForwardDiff, so I’ll give ‘i’ a shot.

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [May 11, 2021, 3:22pm UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/25 "2021-05-11T15:22:32Z")

</div>

Looking at the StatsFuns docs, the log likelihood function should be `normlogpdf((X - Y) / sig) - log(sig)`. Just incase anyone is copy pasting into their own model.

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [May 11, 2021, 3:44pm UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/26 "2021-05-11T15:44:45Z")

</div>

Ah yes, nice catch! Forgot about the logabsdet jacobian term

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [May 18, 2021, 12:26am UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/27 "2021-05-18T00:26:12Z")

</div>

> [@torfjelde](#):
>
> We do have the `.~` which should cover these cases, i.e.
> 
> ```julia
> I .~ BernoulliLogit.(z)
> 
> ```
> 
> but it seems like this only supports a `Vector` of distributions, not arbitrary arrays of distributions. This is surprising to me though, so 'll try getting this fixed!

Btw, this will be fixed as soon as [[Merged by Bors] - Fix use of arrays of Distributions by devmotion · Pull Request #245 · TuringLang/DynamicPPL.jl · GitHub](https://github.com/TuringLang/DynamicPPL.jl/pull/245) is merged, which should be veeery soon:)

---

<div class="post-metadata">

**Author:** ![EvoArt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evoart/32/25357_2.png) [@EvoArt](https://discourse.julialang.org/u/EvoArt)\
**Post date:** [May 18, 2021, 7:59am UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/28 "2021-05-18T07:59:35Z")

</div>

Using Tullio in the log likelihood gave me a big speed up, using forwarddiff. But rhat values were higher and the posterior visibly messy. So I’ve stuck with broadcasting for now. Obviously I could increase the iterations due to the speed up. But seems something is off, so I don’t quite trust it.

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [May 18, 2021, 4:33pm UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/29 "2021-05-18T16:33:52Z")

</div>

Could you share the code using Tullio.jl? Sort of surpirsed here. If it’s indeed just doing the same computation, it shouldn’t be an issue.

Btw, the new version of DPPL now allows you to do the `I .~ BernoulliLogit.(z)` which should also be pretty fast:)

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [May 18, 2021, 4:55pm UTC](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/30 "2021-05-18T16:55:27Z")

</div>

I want to try this. What packages are you using ? Is there any other prep not in the MWE?

[Previous page](https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105.md?page=1)