# R equivalent Turing model from Statistical Rethkinking (model 8.3)

**URL:** <https://discourse.julialang.org/t/r-equivalent-turing-model-from-statistical-rethkinking-model-8-3/69988>\
**Category:** General Usage\
**Created:** [October 18, 2021, 5:09pm UTC](https://discourse.julialang.org/t/r-equivalent-turing-model-from-statistical-rethkinking-model-8-3/69988 "2021-10-18T17:09:09Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![amit1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amit1/32/25318_2.png) [@amit1](https://discourse.julialang.org/u/amit1)\
**Post date:** [October 18, 2021, 5:09pm UTC](https://discourse.julialang.org/t/r-equivalent-turing-model-from-statistical-rethkinking-model-8-3/69988/1 "2021-10-18T17:09:09Z")

</div>

I am trying to convert the following example from the book Statistical Rethinking book to a Turing model:

```r
m8.3 <- quap(
    alist(
        log_gdp_std ~ dnorm( mu , sigma ) ,
        mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) ,
        a[cid] ~ dnorm( 1 , 0.1 ) ,
        b[cid] ~ dnorm( 0 , 0.3 ) ,
        sigma ~ dexp( 1 )
    ) , data=dd )

```

My attempt is as follows:

```julia
@model function m8_3(log_gdp, cid, rugged_st)

	for i = 1 : length(log_gdp)
		a1 ~ Normal(1.0, 0.1)
		b1 ~ Normal(0.0, 0.3)
		a2 ~ Normal(1.0, 0.1)
		b2 ~ Normal(0.0, 0.3)
		σ ~ Exponential(1.0)
		if cid[i] == 1
			μ = a1 + b1 * (rugged_st[i] - 0.215)
		else
			μ = a2 + b2 * (rugged_st[i] - 0.215)
		end
		log_gdp[i] ~ Normal(μ, σ)
	end
end

chain = sample( m8_3(
                      dd.log_gdp_std, 
                      dd.cont_africa, 
                      dd.rugged_std), 
                  NUTS(),10000)

```

The expected output in the book is

```julia
       mean sd 5.5% 94.5%
a[1] 0.89 0.02 0.86 0.91
a[2] 1.05 0.01 1.03 1.07
b[1] 0.13 0.07 0.01 0.25
b[2] -0.14 0.05 -0.23 -0.06
sigma 0.11 0.01 0.10 0.12

```

However my output are somewhat different

```julia
Param a1 a2 b1 b2 sigma
Mean 0.9954 1.0041 -0.0113 0.0064 0.2646

```

Is my conversion correct?

Is there a better way to do it?

* * *

My original version was using an vector for a and b, but it ran in the [following error](https://discourse.julialang.org/t/turing-typeerror-in-typeassert-expected-float64-got-a-value-of-type-forwarddiff-dual-nothing-float64-5/63986). And I do not know how to solve it.

---

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [October 18, 2021, 5:22pm UTC](https://discourse.julialang.org/t/r-equivalent-turing-model-from-statistical-rethkinking-model-8-3/69988/2 "2021-10-18T17:22:05Z")

</div>

I’d do a couple of things:

(1) Lift your alphas, betas, and sigmas out of the loop. I assume it’d be slower, but depending on how Turing translates the model to an actual logpdf, it might also be incorrect (e.g., weighting the prior N times instead of once.) I’m not sure, but I’d be worried.

(2) Since you assign values to mu within the loop, use \mu[i] instead of just \mu. What you are doing now is effectively just creating a single scalar \mu from whatever your last observation i is. Note that this will also require that (outside of the loop) you defined an empty vector \mu, that you then assign to in the loop.

I’m sure a better Turing user than me might have better advice, but that’s what jumps out to me.

---

<div class="post-metadata">

**Author:** ![amit1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amit1/32/25318_2.png) [@amit1](https://discourse.julialang.org/u/amit1)\
**Post date:** [October 18, 2021, 6:05pm UTC](https://discourse.julialang.org/t/r-equivalent-turing-model-from-statistical-rethkinking-model-8-3/69988/3 "2021-10-18T18:05:06Z")

</div>

Your suggestion (1) did indeed make results same as the book

```julia
Params a1 b1 a2 b2 σ   
Mean 0.865713 0.0273194 1.07482 -0.0310709 0.111679

```

I used scalar \mu as thats what was there in the R code, instead of a vector. Besides if i use vector inside Turing model I get the same error

> TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 5}

as I linked in my original post.

My attempted model :

```julia
@model function m8_3(log_gdp, cid, rugged_st)
	a1 ~ Normal(1.0, 0.1)
	b1 ~ Normal(0.0, 0.3)
	a2 ~ Normal(1.0, 0.1)
	b2 ~ Normal(0.0, 0.3)
	μ = zero(log_gdp)
	σ ~ Exponential(1.0)
	for i = 1 : length(log_gdp)
		if cid[i] == 1
			μ[i] = a1 + b1 * (rugged_st[i] - 0.215)
		else
			μ[i] = a2 + b2 * (rugged_st[i] - 0.215)
		end
		log_gdp[i] ~ Normal(μ[i], σ)
	end
end

```

---

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [October 18, 2021, 6:23pm UTC](https://discourse.julialang.org/t/r-equivalent-turing-model-from-statistical-rethkinking-model-8-3/69988/4 "2021-10-18T18:23:22Z")

</div>

Ah, yeah. Since it’s within the for loop, \mu doesn’t need to be a vector, sorry.

But just for future reference, [this part](https://turing.ml/dev/docs/using-turing/performancetips#ensure-that-types-in-your-model-can-be-inferred) of the docs shows how to best initialize an empty array of parameters.