# Different samples distribution after execution

**URL:** <https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077>\
**Category:** Probabilistic Programming\
**Tags:** question, turing\
**Created:** [September 22, 2020, 3:09pm UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077 "2020-09-22T15:09:43Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![mcql](https://avatars.discourse-cdn.com/v4/letter/m/59ef9b/32.png) [@mcql](https://discourse.julialang.org/u/mcql)\
**Post date:** [September 22, 2020, 3:09pm UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077/1 "2020-09-22T15:09:43Z")

</div>

Hello,  
I wrote a this simple Turing model:

```julia
@model function λ(y)
    lambda1 ~ Exponential(alpha)
    lambda2 ~ Exponential(alpha)
    τ ~ Uniform(1, n_count_data)
    for i in 1:length(y)
        y[i] ~ Poisson(τ > i ? lambda1 : lambda2)
    end
    # return τ
end

e = sample(λ(count_data), HMC(0.05,10), 10000)
savefig(plot(e), "sample.png")

```

Every time I run the program (on the same data) the plotted sample distributions looks very different.  
I have upload an example (but as new user I can embedded only one file):  
W

 ![2](https://global.discourse-cdn.com/julialang/original/3X/7/4/74b1f481941ed2b8a5c201b722beca0d383ba674.png)  
I have the same problem using different samplers.  
I don’t know how this could happen.

Thank you.

---

<div class="post-metadata">

**Author:** ![mcql](https://avatars.discourse-cdn.com/v4/letter/m/59ef9b/32.png) [@mcql](https://discourse.julialang.org/u/mcql)\
**Post date:** [September 22, 2020, 3:10pm UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077/2 "2020-09-22T15:10:47Z")

</div>

This is the second plot:

 ![1](https://global.discourse-cdn.com/julialang/original/3X/0/f/0f6d972f4df4cb5180e703e13f88574e3f52971a.png)

---

<div class="post-metadata">

**Author:** ![yiyuezhuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiyuezhuo/32/27101_2.png) [@yiyuezhuo](https://discourse.julialang.org/u/yiyuezhuo)\
**Post date:** [September 22, 2020, 6:48pm UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077/3 "2020-09-22T18:48:22Z")

</div>

Looks like that step size is too large, try to use `HMC(0.01, 10)` or something like that, or use `NUTS()` to let it tune step size automatically for you. A MWE with data is better.

---

<div class="post-metadata">

**Author:** ![mcql](https://avatars.discourse-cdn.com/v4/letter/m/59ef9b/32.png) [@mcql](https://discourse.julialang.org/u/mcql)\
**Post date:** [September 23, 2020, 2:47pm UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077/4 "2020-09-23T14:47:52Z")

</div>

Thank you @yiyuezhuo for the advice.  
Now the sample distributions looks more stable (but often their density converges in different points with HMC()) but the result is different from what I expect.  
I’m trying to wrote in Julia this PyMC example:

> **[Jupyter Notebook Viewer](https://nbviewer.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_PyMC3.ipynb)**
>
> Check out this Jupyter notebook!

And my posterior distributions looks very different:  
(with NUTS() sampler)

 ![fnuts](https://global.discourse-cdn.com/julialang/original/3X/4/a/4afb4ced49a673695da94c67ae84571ac9fedb09.png)  
Of course I’m running the program with the same data.

Thank you!

---

<div class="post-metadata">

**Author:** ![yiyuezhuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiyuezhuo/32/27101_2.png) [@yiyuezhuo](https://discourse.julialang.org/u/yiyuezhuo)\
**Post date:** [September 24, 2020, 7:43am UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077/5 "2020-09-24T07:43:22Z")

</div>

There are two problems in your code:

1. `pm.Exponential` uses a different parameterization compared to `Distributions` or `Turing`, so you should specify `alpha = mean(count_data)` instead of `alpha = 1 / mean(count_data)`.
2. `tau` is in fact a discrete value, though you can use a continuous approximation in some sense, the `HMC` or gradient will not work for `tau`. You should specify a so-called `CompoundStep` in terms of `pymc`.

For example, if you replace the `pm.Metropolis` with `pm.HamiltonianMC`, you will find that parameters use different sampler:

```python
### Mysterious code to be explained in Chapter 3.
with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000,step=step)
"""
Sequential sampling (2 chains in 1 job)
CompoundStep
>Metropolis: [tau]
>Metropolis: [lambda_2]
>Metropolis: [lambda_1]
100%|██████████| 15000/15000 [00:06<00:00, 2194.08it/s]
100%|██████████| 15000/15000 [00:06<00:00, 2469.70it/s]
The number of effective samples is smaller than 25% for some parameters.
"""

### Mysterious code to be explained in Chapter 3.
with model:
    step = pm.HamiltonianMC()
    trace = pm.sample(10000, tune=5000,step=step)
"""
Sequential sampling (2 chains in 1 job)
CompoundStep
>HamiltonianMC: [lambda_2, lambda_1]
>Metropolis: [tau]
100%|██████████| 15000/15000 [00:07<00:00, 1899.67it/s]
100%|██████████| 15000/15000 [00:07<00:00, 2078.96it/s]
The estimated number of effective samples is smaller than 200 for some parameters.
"""

```

My solution, while I don’t know much how `Gibbs` work and `MH` is not exactly equivalent to `pymc.Metropolis`:

```julia
# alpha = 1 / mean(count_data)
alpha = mean(count_data)

@model function λ(y)
    lambda1 ~ Exponential(alpha)
    lambda2 ~ Exponential(alpha)
    τ ~ Uniform(1, n_count_data)
    for i in 1:length(y)
        y[i] ~ Poisson(τ > i ? lambda1 : lambda2)
    end
    # return τ
end

sampler = Gibbs(MH(:τ), HMC(0.01, 5, :lambda1, :lambda2))
e = sample(λ(count_data), sampler, 10000)
plot(e)

```

 ![download](https://global.discourse-cdn.com/julialang/original/3X/7/6/76d22351e2fd5b73089edbbfda293833701c1aef.png)

EDIT: Ok, it seems that just using HMC does not matter much. Random momentum is valid proposal though:

```julia
sampler = HMC(0.01, 10)
e = sample(λ(count_data), sampler, 10000)

```

 ![download (1)](https://global.discourse-cdn.com/julialang/original/3X/1/1/112c8e2502932a12618e15ffc8d30c83dcefc564.png)

---

<div class="post-metadata">

**Author:** ![mcql](https://avatars.discourse-cdn.com/v4/letter/m/59ef9b/32.png) [@mcql](https://discourse.julialang.org/u/mcql)\
**Post date:** [September 24, 2020, 3:14pm UTC](https://discourse.julialang.org/t/different-samples-distribution-after-execution/47077/6 "2020-09-24T15:14:57Z")

</div>

With alpha=mean(count\_data) result are better.  
With Gibbs() sampler I obtain always valid samples.  
Using HMC(0.01, 5) sometimes I obtain something like this:

 ![hu](https://global.discourse-cdn.com/julialang/original/3X/a/c/ac57f12a475b1b4e894d234683199fc8d2c0d911.png)