# Trouble translating a simple PyMC3 model to Turing

**URL:** <https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856>\
**Category:** New to Julia\
**Tags:** turing\
**Created:** [November 21, 2021, 10:27am UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856 "2021-11-21T10:27:12Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 7:44pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/22 "2021-11-21T19:44:51Z")

</div>

> [@dodoplus](#):
>
> (would be great if Turing gave similar commentary…)

This is more a role for the summary function on the Chains type. But yeah, agreed.

I keep getting a segfault:

```julia

signal (11): Segmentation fault
in expression starting at REPL[6]:1
getproperty at ./task.jl:170 [inlined]
produce at /home/dlakelan/.julia/packages/Libtask/8YyQv/src/ctask.jl:154
#1 at /home/dlakelan/.julia/packages/AdvancedPS/31ABw/src/container.jl:12 [inlined]
#8 at /home/dlakelan/.julia/packages/Libtask/8YyQv/src/ctask.jl:85
unknown function (ip: 0x7ffb9cc85dbf)
_jl_invoke at /buildworker/worker/package_linux64/build/src/gf.c:2245 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:2427
jl_apply at /buildworker/worker/package_linux64/build/src/julia.h:1787 [inlined]
start_task at /buildworker/worker/package_linux64/build/src/task.c:878
Allocations: 120671065 (Pool: 120653739; Big: 17326); GC: 188
Segmentation fault

```

running some version of 1.7 I should maybe update it to a more recent one.

Still getting a segfault with 1.7.0-rc3 and latest Turing etc, seems to be a libtask issue, ugh!

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2021, 9:27pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/23 "2021-11-21T21:27:46Z")

</div>

Turing doesn’t work with 1.7.0 yet, unfortunately.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 9:29pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/24 "2021-11-21T21:29:40Z")

</div>

Well, “doesn’t work” seems a bit overstating. I’ve been using it successfully for months, but yeah seems like some functionality isn’t 100%

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2021, 9:42pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/25 "2021-11-21T21:42:49Z")

</div>

It’s not just that – Turing will behave erratically and force-downgrade to version 0.15 IIRC. (Unfortunately, 0.15.0 did not set an upper bound on the newest possible version of Julia, and this can’t really be fixed.)

In reality, I’m kind of shocked it’s working at all, and this might explain some of your bugs. 1.6.4 is the newest non-beta version of Julia, and the newest version supported by Turing. (We hope to fix this in the future by removing our dependence on Libtask, which is tied to specific Julia versions.)

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 9:46pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/26 "2021-11-21T21:46:26Z")

</div>

Gotcha. Thanks for that warning. Unfortunately 1.6 is unusable on my computer due to the old way of doing Registry handling. I pinned Libtask to #master and Turing to 0.19.0 and that allowed the install, it runs and gives meaningful inference in many cases… but I’ll be sure not to do anything important and final until a later version is available.

What exactly does Libtask do for Turing? and what does it take to remove that dependency?

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2021, 9:56pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/27 "2021-11-21T21:56:23Z")

</div>

Libtask is necessary for task-copying, which we use in our SMC/Particle Gibbs implementations; this does explain why PG+SMC are broken. If Turing is pinned to 0.19.0, any non-particle samplers should work.

If SMC/PG are broken for you, you can try Metropolis-Hastings, which should work better.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 10:03pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/28 "2021-11-21T22:03:10Z")

</div>

That explains it then, this is the first time I’ve used PG or SMC. most of what I’ve done is NUTS and a little MH thrown in.

How does MH handle discrete parameters?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 10:35pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/29 "2021-11-21T22:35:34Z")

</div>

Ok on the basis of this discussion and some searching on the Turing manual pages… It seems like this should work:

```julia
trace9 = sample(model9(k01, k10, k11), Gibbs(NUTS(1000,.8,:p0,:p1),MH(:N => RandomWalkProposal(DiscreteUniform(-2,+2)))),1000; init_params=[.6,.2,134],thinning=100)

```

and yet it’s stuck and doesn’t move at all. am I doing something wrong? @ParadaCarleton

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2021, 11:00pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/30 "2021-11-21T23:00:01Z")

</div>

First, I’d disable thinning; thinning only helps if you’re constrained by memory requirements, and I doubt you have a chain that’s long enough that you can’t store it in your computer’s memory. This is also some _very_ aggressive thinning; taking only 1 out of every 100 samples is pretty crazy. This should only leave you with 10 samples at the end.

Also, you should probably be using `RWMH` (Random-walk metropolis-hastings) and not MH.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 11:09pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/31 "2021-11-21T23:09:38Z")

</div>

That doesn’t seem to be the way thinning works, it runs 100000 samples and takes every 100’th one to give me a chain with 1000 samples. The goal was to have it spend a lot of time trying to actually do something useful, but not have to see 100k samples of nothing. In the end, I got 1000 samples of nothing. it never moved

On the other hand, this does work:

```julia

@model function test2(Nobs)
    p ~ Beta(10,10)
    N ~ DiscreteUniform(Nobs,200)
    Nobs ~ Binomial(N,p)
end

ch = sample(test2(22),Gibbs(NUTS(1000,.8,:p),MH(:N => RandomWalkProposal(DiscreteUniform(-1,1)))),1000;init_params=[.4,21],thinning=40)

plot(ch)

```

If I use `RWMH(:N)` it fails with:

```julia
julia> ch = sample(test2(22),Gibbs(NUTS(1000,.8,:p),RWMH(:N)),1000;init_params=[.4,21],thinning=40)
ERROR: MethodError: no method matching getspace(::MetropolisHastings{RandomWalkProposal{false, Symbol}})
Closest candidates are:
  getspace(::Union{DynamicPPL.SampleFromPrior, DynamicPPL.SampleFromUniform}) at ~/.julia/packages/DynamicPPL/RcfQU/src/sampler.jl:9
  getspace(::HMCDA{<:Any, space}) where space at ~/.julia/packages/Turing/nfMhU/src/inference/Inference.jl:440
  getspace(::MH{space}) where space at ~/.julia/packages/Turing/nfMhU/src/inference/Inference.jl:437
  ...

```

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 21, 2021, 11:12pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/32 "2021-11-21T23:12:41Z")

</div>

Oh I see, sorry about that; I must have confused Turing with some other PPL’s thinning.

Was the problem with the initialization, then? It looks like you passed 2 parameters in this case, instead of 3.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 21, 2021, 11:13pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/33 "2021-11-21T23:13:43Z")

</div>

This is a different simpler model, with one continuous and one discrete parameter. I don’t know about that other proposed model above. I don’t know why it won’t move, it could be legitimately nearly impossible to sample from, or it could be a software issue.

---

<div class="post-metadata">

**Author:** ![bdeonovic](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bdeonovic/32/3928_2.png) [@bdeonovic](https://discourse.julialang.org/u/bdeonovic)\
**Post date:** [November 22, 2021, 4:03pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/34 "2021-11-22T16:03:02Z")

</div>

Just a sanity check:  
are we sure that the parameterization from PyMC3 ` pm.HyperGeometric('y', N=N, k=k, n=n, observed=x)` is equivalent to Turing’s `Hypergeometric(k, N - k, n)`?

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [November 22, 2021, 5:28pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/35 "2021-11-22T17:28:06Z")

</div>

> [@bdeonovic](#):
>
> are we sure that the parameterization from PyMC3 ` pm.HyperGeometric('y', N=N, k=k, n=n, observed=x)` is equivalent to Turing’s `Hypergeometric(k, N - k, n)` ?

It _should_ be.  
The PyMC3 documentation [says](https://docs.pymc.io/en/v3/api/distributions/discrete.html):

> The probability of `x` successes in a sequence of `n` bernoulli trials taken without replacement from a population of `N` objects, containing `k` good (or successful or Type I) objects.

while Julia’s Distribution docs [say](https://juliastats.org/Distributions.jl/v0.14/univariate.html#Distributions.Hypergeometric):

> A _Hypergeometric distribution_ describes the number of successes in `n` draws without replacement from a finite population containing `s` successes and `f` failures

In other words, `N=s+f`. In Julia we pass `s` and `f` and in Python we pass `N` and `s`.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 22, 2021, 8:33pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/37 "2021-11-22T20:33:25Z")

</div>

Ok In this multinomial model I think I’ve discovered one problem:

```julia
@model function model9(k01, k10, k11)
    p0 ~ Beta(1, 1)
    p1 ~ Beta(1, 1)
    num_seen = k01 + k10 + k11
    N ~ DiscreteUniform(num_seen, 350)
    
    q0 = 1-p0
    q1 = 1-p1
    ps = [q0*q1, q0*p1, p0*q1, p0*p1]
    
    y = [N - num_seen, k01, k10, k11]
    y ~ Multinomial(N, ps)
end

```

by looking at the chain output, I see that y is NOT treated like data, instead y is **overwritten** with samples from the Multinomial. This is because y doesn’t appear in the arguments to the model function (and it can’t obviously).

So I modified it as follows:

```julia
@model function model9(k01, k10, k11)
    p0 ~ Beta(1, 1)
    p1 ~ Beta(1, 1)
    num_seen = k01 + k10 + k11
    N ~ DiscreteUniform(num_seen, 350)
    
    q0 = 1-p0
    q1 = 1-p1
    ps = [q0*q1, q0*p1, p0*q1, p0*p1]
    
    y = [N - num_seen, k01, k10, k11]
    Turing.@addlogprob!(logpdf(Multinomial(N, ps),y))
end

```

and now using:

```julia
trace9 = sample(model9(k01, k10, k11), Gibbs(MH(:N => RandomWalkProposal(DiscreteUniform(-1,1))),NUTS(1000,.8,:p0,:p1)),1000; init_params=[.09,.04,152],thinning=10)

```

it samples as follows:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/2/82d512e643f563dbf3340e50356f41f25b944d08.png)

Or with a bit more samples and thinning = 100

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/d/6d98146932bc300b90e238e580dd278afbf91bc4.png)

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 22, 2021, 10:49pm UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/38 "2021-11-22T22:49:25Z")

</div>

Also for the original model, the following seems to work:

```julia

k = 23
n = 19
x = 4

@model function model6(y,k,n)
    N ~ DiscreteUniform(50, 500)
    y ~ Hypergeometric(k, N - k, n)
end

trace6 = sample(model6(x,k,n), MH(:N => AdvancedMH.RandomWalkProposal(DiscreteUniform(-1,1))), 1000; thinning=1000);

plot(trace6)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/4/04ae056fccc7418b41478e297d538260187bcbe0.png)

```julia
Chains MCMC chain (1000×2×1 Array{Float64, 3}):

Iterations = 1:1000:999001
Number of chains = 1
Samples per chain = 1000
Wall duration = 2.79 seconds
Compute duration = 2.79 seconds
parameters = N
internals = lp

Summary Statistics
  parameters mean std naive_se mcse ess rhat ess_per_sec 
      Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64 

           N 193.0430 94.2292 2.9798 14.4653 24.1587 1.0606 8.6466

Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5% 
      Symbol Float64 Float64 Float64 Float64 Float64 

           N 70.9750 125.0000 170.0000 235.2500 426.0500

```

Or using a wider proposal, for much more effective sampling:

```julia

trace6 = sample(model6(x,k,n), MH(:N => AdvancedMH.RandomWalkProposal(DiscreteUniform(-10,10))), 1000; thinning=100);

plot(trace6)
trace6

Chains MCMC chain (1000×2×1 Array{Float64, 3}):

Iterations = 1:100:99901
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.29 seconds
Compute duration = 0.29 seconds
parameters = N
internals = lp

Summary Statistics
  parameters mean std naive_se mcse ess rhat ess_per_sec 
      Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64 

           N 165.4380 80.1469 2.5345 7.3142 114.7284 1.0002 391.5646

Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5% 
      Symbol Float64 Float64 Float64 Float64 Float64 

           N 66.9750 109.0000 144.0000 200.2500 386.0500

```

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 23, 2021, 2:44am UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/39 "2021-11-23T02:44:13Z")

</div>

I’m glad this works! That’s definitely some very weird and unexpected behavior.

@yebai do you know what’s up with this? Looks like a bug with how `@model` handles `~`.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 23, 2021, 2:53am UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/40 "2021-11-23T02:53:39Z")

</div>

I’m not 100% sure but I think this is a feature or a intentional design decision. In general the left hand side of ~ is considered a parameter unless it appears in the arguments of the model function.

I’m not clear on what the intended way to handle your example problem is. However the `@addlogprob!` is a general purpose way to handle stuff like this.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [November 23, 2021, 3:19am UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/41 "2021-11-23T03:19:27Z")

</div>

`~` should work even if it’s not in the arguments – actually, we’re starting to work on moving data out of model arguments. If `y` is already defined, it seems very weird to me for `y` to be overwritten by a sampled variable.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 23, 2021, 4:23am UTC](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856/42 "2021-11-23T04:23:22Z")

</div>

> [@ParadaCarleton](#):
>
> moving data out of model arguments

Where would the data go instead?

[Previous page](https://discourse.julialang.org/t/trouble-translating-a-simple-pymc3-model-to-turing/71856.md?page=1)