# Help converting pymc3 switchpoint analysis code to Turing.jl

**URL:** <https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880>\
**Category:** New to Julia\
**Tags:** question, turing, machine-learning\
**Created:** [May 21, 2020, 11:50am UTC](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880 "2020-05-21T11:50:04Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![3210jr](https://avatars.discourse-cdn.com/v4/letter/3/779978/32.png) [@3210jr](https://discourse.julialang.org/u/3210jr)\
**Post date:** [May 21, 2020, 11:50am UTC](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880/1 "2020-05-21T11:50:05Z")

</div>

Hey all,

I’ve been trying to practice turing.jl and Julia as a whole and its been a great experience so far!

I am having some trouble coming from pymc3 background, any help and direction is much appreciated!

Given the code from the (amazing) book on Bayesian methods for Hackers:

```julia
import pymc3 as pm
import theano.tensor as tt

with pm.Model() as model:
    alpha = 1.0/count_data.mean() # Recall count_data is the
                                   # variable that holds our txt counts
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)
    
    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)

    idx = np.arange(n_count_data) # Index
    lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)

    observation = pm.Poisson("obs", lambda_, observed=count_data)

    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000,step=step)

```

my translation is:

```julia
@model text_model(y) = begin
    α = 1/mean(y)
    λ₁ ~ Exponential(α)
    λ₂ ~ Exponential(α)
    
    τ ~ DiscreteUniform(1, length(y))
    
    function equality(idx)
        if τ > idx
            return λ₁
        end
        return λ₂
    end
    
    times = 1:length(y)
    λ = map(equality, times)

    for idx in times
        y[idx] ~ Poisson(λ[idx])
    end
    return τ
end

model = model(data.count)
trace = sample(model, Gibbs(HMC(0.01, 10, :λ₁, :λ₂), PG(20, :τ)), 1000)

```

There are no errors (other than warning of rejected proposals for numerical reasons once in a while) and the sampling goes on without a problem.

The problem is that I get the wrong results 😑 .

Question:  
Do you think the code translation makes sense? not super sure about the switchpoint method I am using.

Thanks in advance Julia and Turing Communities!

---

<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:** [May 21, 2020, 1:12pm UTC](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880/2 "2020-05-21T13:12:26Z")

</div>

What exactly is wrong in the results? Also, in the python script the upper bound of the discrete uniform is n\_count - 1 while it is n\_count in the Turing script. Maybe that is a problem? Hard to examine without data. Is the data available somewhere?

---

<div class="post-metadata">

**Author:** ![3210jr](https://avatars.discourse-cdn.com/v4/letter/3/779978/32.png) [@3210jr](https://discourse.julialang.org/u/3210jr)\
**Post date:** [May 21, 2020, 1:16pm UTC](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880/3 "2020-05-21T13:16:24Z")

</div>

Thanks for reaching out @trappmartin , the reason for python to be n\_count -1 is since its 0 indexed - therefore the length of a list (data) is different from the index of the last item.

And the code is talked about here: [Jupyter Notebook Viewer](https://nbviewer.jupyter.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_PyMC3.ipynb)

And does the switch statement make sense to you? I this how you would implement `pm.switch` in turing?

Thanks again!

---

<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 21, 2020, 1:56pm UTC](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880/4 "2020-05-21T13:56:58Z")

</div>

This seems to work fine:

```julia
@model text_model(y) = begin
    α = 1/mean(y)
    λ₁ ~ Exponential(α)
    λ₂ ~ Exponential(α)
    τ ~ DiscreteUniform(1, length(y))
    for idx in 1:length(y)
        y[idx] ~ Poisson(τ > idx ? λ₁ : λ₂)
    end
    return τ
end
sample(text_model(rand(1:100, 100)), MH(), 1000)

```

Edit: Man, Turing is neat!! Looking at that Python code above brought back some unpleasant memories!

---

<div class="post-metadata">

**Author:** ![jonathan\_skai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonathan_skai/32/4314_2.png) [@jonathan\_skai](https://discourse.julialang.org/u/jonathan_skai)\
**Post date:** [July 31, 2021, 2:29pm UTC](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880/6 "2021-07-31T14:29:14Z")

</div>

Your results will match approximately.  
Change to α = mean(y) in your code, since Exponential(α) in Distributions.jl is different from the PyMC3 version.