# Bayesian inference in the presence of intractable integral

**URL:** <https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061>\
**Category:** Probabilistic Programming\
**Tags:** bayesian-inference, monte-carlo\
**Created:** [May 24, 2020, 9:40am UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061 "2020-05-24T09:40:07Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 24, 2020, 9:40am UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/1 "2020-05-24T09:40:07Z")

</div>

Hi guys, recently i’ve learned about the powerful packages for bayesian inference in julia  
i would like some tips on how to attack this problem, given this probabiliy distribution over the angles \theta:

P(\vec \theta\_1,\vec \theta\_2|J,\vec h\_1,\vec h\_2) = \frac{W(\vec \theta\_1,\vec \theta\_2|J,\vec h\_1,\vec h\_2) }{Z(J,\vec h\_1,\vec h\_2)} 

where W and Z are defined as

W(\vec \theta\_1,\vec \theta\_2|J,\vec h\_1,\vec h\_2) = \exp\big(\sum\_{i,j=1}^{2} J\_{i,j} \Sigma(\vec\theta\_i) \cdot \Sigma(\vec\theta\_j)+\sum\_{i=1}^2 h\_i \cdot \Sigma(\vec\theta\_i)\big)

Z(J,\vec h\_1,\vec h\_2)=\int\_{[0,2\pi]^{N\_{\theta\_1}+N\_{\theta\_2}}} d\theta\_1 \, d\theta\_2\,W(\vec \theta\_1,\vec \theta\_2|J,\vec h\_1,\vec h\_2) 

finally the \Sigma\_i are defined as

\Sigma(\theta)=\frac{1}{N\_\theta}\sum\_{j=1}^{N\_\theta} \big(\cos(\vec \theta \cdot \hat e\_j),\sin(\vec \theta \cdot \hat e\_j)\big)

I have some observations of these sets of angles \vec \theta\_1,\vec \theta\_2 and i want to infer the posterior distribution of the matrix J and the vectors h\_{1/2}. but i have no clue on how to leverage the marvelous tools in the julia ecosystem like Turing.jl etc… since the integral in Z will unavoidably appear in the `logpdf`.

if this question is considered off-topic i will delete it. Thank you in advance

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 24, 2020, 3:59pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/2 "2020-05-24T15:59:03Z")

</div>

What is N\_\theta?

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 24, 2020, 4:00pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/3 "2020-05-24T16:00:42Z")

</div>

The length of the \theta vector, Sorry for the bad notation

---

<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:** [May 24, 2020, 4:10pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/4 "2020-05-24T16:10:24Z")

</div>

IIUC, you can define your own distribution for `theta = [theta1; theta2]` parameterized by `J`, `h1` and `h2`. Then you can define a Turing model as follows:

```julia
@model mymodel(theta) = begin
    J ~ J's prior
    h1 ~ h1's prior
    h2 ~ h2's prior
    theta ~ MyDist(J, h1, h2)
end

```

Here is how to define your own distribution and make it Turing-compatible [Advanced Usage](https://turing.ml/dev/docs/using-turing/advanced#1-define-the-distribution-type).

The main part in your case is defining the `logpdf` function because the pdf functions seems involved.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 24, 2020, 4:12pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/5 "2020-05-24T16:12:30Z")

</div>

That Is very nice, but…  
Logpdf involves Z, which Is the intractable integral…

---

<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:** [May 24, 2020, 5:10pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/6 "2020-05-24T17:10:32Z")

</div>

How large are N\_{\theta\_1} and N\_{\theta\_2}?

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 24, 2020, 5:12pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/7 "2020-05-24T17:12:15Z")

</div>

Each one roughly \sim 500

---

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [May 24, 2020, 6:30pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/8 "2020-05-24T18:30:29Z")

</div>

I work a lot with models that have intractable normalizing constants in the likelihood (specifically, [ERGMs](https://en.wikipedia.org/wiki/Exponential_random_graph_models)). For that kind of problem, I don’t know that you have a choice except to write your own algorithm using e.g., ABC, auxiliary variable MCMC, things like the [DMH sampler](https://www.tandfonline.com/doi/abs/10.1080/00949650902882162) algorithm, etc. I don’t think there’s much anywhere (Julia or otherwise) that does “off-the-shelf” for such models. Well, [Knet](https://denizyuret.github.io/Knet.jl/latest/nce/) _may_ have NCE (noise-contrastive estimation; it’s not clear to me if they actually have it, or if they had it a long time ago, and the docs are outdated). NCE works a lot like MCMC-MLE, and is good for problems with intractable normalizing constants.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 24, 2020, 8:53pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/9 "2020-05-24T20:53:56Z")

</div>

How big are the eigenvalues of the J matrix? Maybe you can do a steepest descent approximation

---

<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:** [May 24, 2020, 9:43pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/10 "2020-05-24T21:43:09Z")

</div>

Could you please give a reference paper or book for this method?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 24, 2020, 10:54pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/11 "2020-05-24T22:54:21Z")

</div>

Wikipedia has some references: [Method of steepest descent - Wikipedia](https://en.m.wikipedia.org/wiki/Method_of_steepest_descent)

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 25, 2020, 5:22am UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/12 "2020-05-25T05:22:54Z")

</div>

That Is the core of the inference problem, they can be 0 or very big, thanks anyway

---

<div class="post-metadata">

**Author:** ![robsmith11](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robsmith11/32/29641_2.png) [@robsmith11](https://discourse.julialang.org/u/robsmith11)\
**Post date:** [May 25, 2020, 6:18am UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/13 "2020-05-25T06:18:38Z")

</div>

I’m still a beginner at Bayesian inference, but have had good results using Approximate Bayesian Computation (ABC) for parameter inference for some mixture models with Markovian time dependencies that make the likelihood intractable.

I’ve been using the parallel implementation of ABC-SMC from [GitHub - marcjwilliams1/ApproxBayes.jl: Approximate Bayesian Computation (ABC) algorithms for likelihood free inference in julia](https://github.com/marcjwilliams1/ApproxBayes.jl)

I’ve been meaning to try out [GitHub - tanhevg/GpABC.jl](https://github.com/tanhevg/GpABC.jl) but so far ApproxBayes.jl has been good enough.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 25, 2020, 2:43pm UTC](https://discourse.julialang.org/t/bayesian-inference-in-the-presence-of-intractable-integral/40061/14 "2020-05-25T14:43:53Z")

</div>

@opera_malenky  
thank you for the pointers, the DMH algorithm is very interesting, i will try to implement it, although it is very sad that such an algorithm is not available in the main general purpose packages that deal with bayesian computation.

@robsmith11  
thank you for the suggestion on ApproxBayes.jl, that could be a winning approach if DHM turns out too slow.