# PyMC Models to Turing.jl

**URL:** https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574
**Category:** Probabilistic Programming
**Tags:** question, turing, bayesian-inference
**Created:** [December 10, 2020, 12:30am UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574 "2020-12-10T00:30:38Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![Stefan\_Bringuier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefan_bringuier/32/20132_2.png) [@Stefan\_Bringuier](https://discourse.julialang.org/u/Stefan_Bringuier)
#### Post date: [December 10, 2020, 12:30am UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574/1 "2020-12-10T00:30:38Z")

</div>

I’m learning (beginner) some Bayesian methods from books like [Bayesian Methods for Hackers, C. Davidson-Pilon](https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers) and [Bayesian Analysis with Python](https://www.packtpub.com/product/bayesian-analysis-with-python-second-edition/9781789341652). However, I prefer to do most of my workflows using Julia and I would like to do things using packages such as **Turing.jl**. Since I’m new to probalistic programming its not clear how to translate many of the examples from PyMC to Turing.jl, are there any guides or ways of thinking to help with this?

---

<div class="post-metadata">

### Author: ![Stefan\_Bringuier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefan_bringuier/32/20132_2.png) [@Stefan\_Bringuier](https://discourse.julialang.org/u/Stefan_Bringuier)
#### Post date: [December 10, 2020, 1:10am UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574/2 "2020-12-10T01:10:29Z")

</div>

I’m adding on here a specific example that I’d like to reproduce with Turing.jl from Bayesian Methods for Hackers:  
[Example: Kaggle contest on \*Observing Dark World](https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter5_LossFunctions/Ch5_LossFunctions_PyMC3.ipynb)

---

<div class="post-metadata">

### Author: ![cpfiffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cpfiffer/32/208747_2.png) [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)
#### Post date: [December 10, 2020, 8:45pm UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574/3 "2020-12-10T20:45:58Z")

</div>

Check out the [tutorials](https://turing.ml/stable/tutorials/) section for Turing.jl to see if anything tickles your fancy.

---

<div class="post-metadata">

### Author: ![Stefan\_Bringuier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefan_bringuier/32/20132_2.png) [@Stefan\_Bringuier](https://discourse.julialang.org/u/Stefan_Bringuier)
#### Post date: [December 11, 2020, 7:21pm UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574/4 "2020-12-11T19:21:14Z")

</div>

Thanks, the PyMC3 notebook example is pretty useful. I’m still having some difficulty though conceptualizing how setup something like this with Turing.jl:

```python
import pymc3 as pm

...

with pm.Model() as model:
    #set the size of the halo's mass
    mass_large = pm.Uniform("mass_large", 40, 180)
    
    #set the initial prior position of the halos, it's a 2-d Uniform dist.
    halo_position = pm.Uniform("halo_position", 0, 4200, shape=(1,2))
    
    mean = pm.Deterministic("mean", mass_large /\
            f_distance(T.as_tensor(data[:,:2]), halo_position, 240)*\
            tangential_distance(T.as_tensor(data[:,:2]), halo_position))
    
    ellpty = pm.Normal("ellipcity", mu=mean, tau=1./0.05, observed=data[:,2:])

```

It’s not clear to me how to define a Turing.jl model with variables that are determined from priors (e.g. `mean` in the code snippet above). Tis could just be that I’m new to probabilistic programming and Turing.jl.

Thanks in advance.

---

<div class="post-metadata">

### Author: ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [December 11, 2020, 9:35pm UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574/5 "2020-12-11T21:35:49Z")

</div>

In general, you can just write Julia code that uses the probabilistic variables as though they were ordinary variables. So `mean` in the model above could be written similarly, minus the `pm.Deterministic` and `T.as_tensor` boilerplate.

I started translating the dark world example, but ended up rewriting it in a way that made more sense to me. Have not checked it for correctness or optimized it, but this at least shows conceptually one way to do it:

```julia
using Turing
using LinearAlgebra: norm
using Plots, StatsPlots

# Copied and pasted from the dark sky example linked above. Columns are
# galaxy x and y positions, and the two components of their eccentricity.
# Transposed to get things column-major.
data = [1.62690000e+02 1.60006000e+03 1.14664000e-01 -1.90326000e-01;
    2.27228000e+03 5.40040000e+02 6.23555000e-01 2.14979000e-01;
    3.55364000e+03 2.69771000e+03 2.83527000e-01 -3.01870000e-01]'
galaxy_positions = data[1:2, :]
galaxy_eccs = data[3:4, :]

f_distance(galaxy_pos, halo_pos, c) = 1 / min(norm(galaxy_pos .- halo_pos), c)

function tangential_distance(galaxy_pos, halo_pos)
    Δ = galaxy_pos .- halo_pos
    t = 2 * atan(Δ[2] / Δ[1])
    return [-cos(t), -sin(t)]
end

function predict_ecc(galaxy_pos, halo_positions, halo_masses, c)
    d = [tangential_distance(galaxy_pos, hp) for hp in eachcol(halo_positions)]
    f = [f_distance(galaxy_pos, hp, c) for hp in eachcol(halo_positions)]
    return(reduce(+, d .* halo_masses .* f))
end
    
@model function DarkSky(galaxy_positions, galaxy_eccs, nhalos=3, c=240)
    n = size(galaxy_positions, 2)
    halo_positions ~ filldist(Uniform(0, 4200), 2, nhalos)
    halo_masses ~ filldist(Uniform(40, 180), nhalos)
    σ2 = 0.05

    # for each galaxy, predict its observed eccentricity and compare it with the data,
    # assuming m.v. normal errors with diagonal covariance
    for i in 1:n
        μ_ecc = predict_ecc(galaxy_positions[:, i], halo_positions, halo_masses, 240)
        galaxy_eccs[:, i] ~ MvNormal(μ_ecc, sqrt(σ2))
    end
end

model = DarkSky(galaxy_positions, galaxy_eccs, 3, 240)
chain = sample(model, NUTS(), 1000)

```

---

<div class="post-metadata">

### Author: ![Stefan\_Bringuier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefan_bringuier/32/20132_2.png) [@Stefan\_Bringuier](https://discourse.julialang.org/u/Stefan_Bringuier)
#### Post date: [December 11, 2020, 9:45pm UTC](https://discourse.julialang.org/t/pymc-models-to-turing-jl/51574/6 "2020-12-11T21:45:14Z")

</div>

Thanks a bunch! This workflow translation was exactly the kind of thing I was needing help with. Very much appreciated.
