# Is there a turing alternative to pm.Deterministic from pymc3?

**URL:** <https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667>\
**Category:** New to Julia\
**Tags:** turing, distributions\
**Created:** [May 3, 2020, 11:58am UTC](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667 "2020-05-03T11:58:42Z")\
**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 3, 2020, 11:58am UTC](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/1 "2020-05-03T11:58:42Z")

</div>

Hey all,

I am looking an alternative to `pm.Deterministic` so that I can store and return a computed value.

Thanks

---

<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 5, 2020, 8:55am UTC](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/2 "2020-05-05T08:55:35Z")

</div>

Hi,  
You can do something like this:

```julia
using Turing, Random

# Define a deterministic distribution
struct Determin{T<:Real} <: ContinuousUnivariateDistribution
  val::T
end

Distributions.rand(rng::AbstractRNG, d::Determin) = d.val
Distributions.logpdf(d::Determin, x::T) where T<:Real = zero(x)

@model function testmodel(x)
  y ~ Normal()
  m ~ Determin(y * 2)
  for i in eachindex(x)
    x[i] ~ Normal(m)
  end
end

# instantiate a model
model = testmodel(randn(10))

# inference
chain = sample(model, SMC(), 1000)

```

And if I would recall how to properly define the bijection for a custom distribution, you could also use HMC or NUTS. @torfjelde could you jump in?

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [May 7, 2020, 5:29pm UTC](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/3 "2020-05-07T17:29:00Z")

</div>

As @trappmartin says, you need to define a transformation from unconstrained to constrained, but in this case that’s just going to be the identity operation, so you only the need additional bit of code:

```julia
Bijectors.bijector(d::Determin) = Identity{0}() # <= `0` indicates 0-dimensional, i.e. univariate 

```

and you’ll be able to use HMC and NUTS:)

---

<div class="post-metadata">

**Author:** ![mateuszbaran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mateuszbaran/32/221842_2.png) [@mateuszbaran](https://discourse.julialang.org/u/mateuszbaran)\
**Post date:** [June 23, 2020, 6:42am UTC](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/4 "2020-06-23T06:42:17Z")

</div>

I’m trying to use this solution in my model but so far it ends very poorly. Here is a simplified version:

```julia
using LinearAlgebra
using Turing
using Distributions
using Bijectors
using Random

a1 = [0.05, 0.03, 0.02, 0.03, 0.04, 0.09, 0.13, 0.16]

# Define a deterministic distribution
struct Determin{T<:Real} <: ContinuousUnivariateDistribution
  val::T
end

Distributions.rand(rng::AbstractRNG, d::Determin) = d.val
Distributions.logpdf(d::Determin, x::T) where T<:Real = zero(x)

Bijectors.bijector(d::Determin) = Identity{0}() # <= `0` indicates 0-dimensional, i.e. univariate

@model model1(a, ::Type{T} = Float64) where {T} = begin
    ma ~ Normal(0.0, 0.5)
    for i in 1:length(a)
        a[i] ~ Normal(ma, 0.05)
    end

    b ~ Determin(-exp(-2*ma))
    return ma, b
end

function test_model1()
    # sampler = HMC(0.001, 10)
    # sampler = HMCDA(1000, 0.75, 0.05)
    sampler = NUTS(0.65)
    c1 = sample(model1(a1), sampler, 2000)

    println(c1)
end

```

Results typically look like this:

```julia
julia> test_model1()
┌ Info: Found initial step size
└ ϵ = 0.025
Object of type Chains, with data of type 1000×14×1 Array{Float64,3}

Iterations = 1:1000
Thinning interval = 1
Chains = 1
Samples per chain = 1000
internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
parameters = b, ma

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters mean std naive_se mcse ess r_hat
  ────────── ───────────────────────── ──────────────────────── ─────────────────────── ──────────────────────── ──────── ──────
           b 13816403112283807744.0000 7133491167269512192.0000 225580797572648256.0000 2172001283591027712.0000 4.8969 1.3982
          ma 0.0690 0.0179 0.0006 0.0007 830.2122 0.9998

Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5%
  ────────── ───────────────────────── ───────────────────────── ───────────────────────── ───────────────────────── ─────────────────────────
           b -5026632495806943232.0000 12090906044107700224.0000 15653120955249981440.0000 18317459011682740224.0000 22252367610115112960.0000
          ma 0.0349 0.0565 0.0694 0.0811 0.1047

```

HMC and HMCDA samplers also give similarly bad results. Did I make something stupid here or is this a bug?

EDIT: I just tried SMC and it actually works.

---

<div class="post-metadata">

**Author:** ![mateuszbaran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mateuszbaran/32/221842_2.png) [@mateuszbaran](https://discourse.julialang.org/u/mateuszbaran)\
**Post date:** [June 25, 2020, 8:41am UTC](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/5 "2020-06-25T08:41:51Z")

</div>

I have reported this as an issue in the Turing repository: [https://github.com/TuringLang/Turing.jl/issues/1335](https://github.com/TuringLang/Turing.jl/issues/1335) .