# How to extract the Standard Deviation out of Turing.jl

**URL:** <https://discourse.julialang.org/t/how-to-extract-the-standard-deviation-out-of-turing-jl/31519>\
**Category:** Statistics\
**Created:** [November 26, 2019, 6:38am UTC](https://discourse.julialang.org/t/how-to-extract-the-standard-deviation-out-of-turing-jl/31519 "2019-11-26T06:38:39Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![StevenSiew](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevensiew/32/218393_2.png) [@StevenSiew](https://discourse.julialang.org/u/StevenSiew)\
**Post date:** [November 26, 2019, 6:38am UTC](https://discourse.julialang.org/t/how-to-extract-the-standard-deviation-out-of-turing-jl/31519/1 "2019-11-26T06:38:39Z")

</div>

I am trying to extract the Standard Deviation out of Turing for the parameter p but it kept giving me an error. I was tired and frustrated and I could not think straight. I guess I need to take a break.

The guide here is of no help [https://turing.ml/dev/tutorials/](https://turing.ml/dev/tutorials/)

I get the error

```julia
ERROR: LoadError: MethodError: no method matching iterate(::Chains{Union{Missing, Real},Missing,NamedTuple{(:parameters,),Tuple{Array{String,1}}},NamedTuple{(:hashedsummary,),Tuple{Base.RefValue{Tuple{UInt64,Array{ChainDataFrame,1}}}}}})
Closest candidates are:
  iterate(::Core.SimpleVector) at essentials.jl:600
  iterate(::Core.SimpleVector, ::Any) at essentials.jl:600
  iterate(::ExponentialBackOff) at error.jl:218

```

for

```julia
# personally extracting the statistics for the parameter p
p_mean = mean(chain[:p])[:mean][1] # Extracting the mean for parameter p
p_std = std(chain[:p])[:std][1] # Extracting the std for parameter p

println("Parameter p has a mean of $(p_mean) and a standard deviation of $(p_std)")

```

Here is the entire source code to reproduce the error.

```julia
# Import libraries.
using Turing
using Distributions, Random, Plots

# Set the true probability of heads in a coin.
p_true = 0.5

# We has 7 coin toss observations.
NumOfCoinToss = 7

# Draw data from a Bernoulli distribution, i.e. draw heads or tails.
Random.seed!(12345)
data = rand(Bernoulli(p_true), NumOfCoinToss)
@show data

println("The result of ",NumOfCoinToss," coin tosses are:")
commaflag=false
for coinresult in data
    global commaflag
    print(commaflag ? ", " : "")
    print(coinresult == 1 ? "head" : "tail")
    commaflag=true
end
println()

# Now we use Turing to find the unknown probability p

# Declare our Turing model.
@model coinflip(data) = begin
    # Our prior belief about the unknown parameter p which is
    # the probability of heads in a coin.
    p ~ Beta(1, 1)

    # The number of observations.
    N = length(data)
    for n in 1:N
        # Heads or tails of a coin are drawn using
        # a process which have a Bernoulli distribution
        data[n] ~ Bernoulli(p)
    end
end;

# Settings of the Hamiltonian Monte Carlo (HMC) sampler.
# If you never used HMC before then do not change the value
# of ϵ and τ from 0.05 and 10
samples = 1_000
ϵ = 0.05 # leapfrog step size
τ = 10 # leapfrog step numbers

# Start sampling.
chain = sample(coinflip(data), HMC(ϵ, τ), samples);

# Summarise results
describe(chain)

# Plot a summary of the sampling process for the parameter p, i.e. the probability of heads in a coin.
histogram(chain[:p])

# personally extracting the statistics for the parameter p
p_mean = mean(chain[:p])[:mean][1] # Extracting the mean for parameter p
p_std = std(chain[:p])[:std][1] # Extracting the std for parameter p

println("Parameter p has a mean of $(p_mean) and a standard deviation of $(p_std)")

```

---

<div class="post-metadata">

**Author:** ![StevenSiew](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevensiew/32/218393_2.png) [@StevenSiew](https://discourse.julialang.org/u/StevenSiew)\
**Post date:** [November 26, 2019, 9:37am UTC](https://discourse.julialang.org/t/how-to-extract-the-standard-deviation-out-of-turing-jl/31519/2 "2019-11-26T09:37:23Z")

</div>

Looks like the info on Chains is found somewhere else

> **[GitHub - TuringLang/MCMCChains.jl: Types and utility functions for summarizing Markov...](https://github.com/TuringLang/MCMCChains.jl)**
>
> Types and utility functions for summarizing Markov chain Monte Carlo simulations

```julia-auto
# construction of a Chains object with no names
Chains(val::AbstractArray{A,3};
        start::Int=1,
        thin::Int=1,
        evidence = 0.0,
        info=NamedTuple())

# construction of a chains object with new names
Chains(val::AbstractArray{A,3},
        parameter_names::Vector{String},
        name_map = copy(DEFAULT_MAP);
        start::Int=1,
        thin::Int=1,
        evidence = 0.0,
        info=NamedTuple())

# Indexing a Chains object
chn = Chains(...)
chn_param1 = chn[:,2,:] # returns a new Chains object for parameter 2
chn[:,2,:] = ... # set values for parameter 2

```

### The `get` Function

MCMCChains provides a `get` function designed to make it easier to access parameters `get(chn, :P)` returns a `NamedTuple` which can be easy to work with.

Is there an easier way to get the information than to write this strange incantation below?

```julia-auto
julia> p_mean = chain.info[1][][2][1][:mean][]
0.5556262828496785

julia> p_std = chain.info[1][][2][1][:std][]
0.16302793812807917

```

---

<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:** [November 26, 2019, 2:24pm UTC](https://discourse.julialang.org/t/how-to-extract-the-standard-deviation-out-of-turing-jl/31519/3 "2019-11-26T14:24:38Z")

</div>

You can pull the actual values for the chain out using `chain.value` – for a single parameter, the quicker way to extract summary stats is to make them yourself:

```julia
using Statistics

val_p = chain[:p].value
p_mean = mean(val_p)
p_std = std(val_p)

```

If you have multiple prameters you want the stats for, a good way to do this is with `summarize`, which maps functions across parameters.

```julia
chain = Chains(rand(100,3,3))
sumstats = summarize(chain, mean, std)
means = sumstats[:mean]
stdevs = sumstats[:std]

```

`sumstats` looks like this:

```julia
│ Row │ parameters │ mean │ std │
│ │ Symbol │ Float64 │ Float64 │
├─────┼────────────┼──────────┼──────────┤
│ 1 │ Param1 │ 0.513241 │ 0.271461 │
│ 2 │ Param2 │ 0.471095 │ 0.293417 │
│ 3 │ Param3 │ 0.501808 │ 0.290675 │

```

---

<div class="post-metadata">

**Author:** ![StevenSiew](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevensiew/32/218393_2.png) [@StevenSiew](https://discourse.julialang.org/u/StevenSiew)\
**Post date:** [November 27, 2019, 12:35am UTC](https://discourse.julialang.org/t/how-to-extract-the-standard-deviation-out-of-turing-jl/31519/4 "2019-11-27T00:35:18Z")

</div>

Thank you very much, I came up with this to extract the property for for multiple parameters

If you have multiple prameters you want the stats for, a good way to do this is with summarize, which maps functions across parameters.

```julia
chain = Chains(rand(100,3,3))
sumstats = summarize(chain, mean, std)
means = sumstats[:mean]
stdevs = sumstats[:std]
varpos(v::Symbol,vec::Vector) = first(findall(x->x==v,vec))
dict_sumstats=Dict(:mean=>2,:std=>3)
getParamProperty(param,property) = sumstats[varpos(param,sumstats[:,1]),dict_sumstats[property]]

getParamProperty(:Param2,:std)
0.2761685499505822

```

sumstats looks like this:

```julia
│ Row │ parameters │ mean │ std │
│ │ Symbol │ Float64 │ Float64 │
├─────┼────────────┼──────────┼──────────┤
│ 1 │ Param1 │ 0.489869 │ 0.288344 │
│ 2 │ Param2 │ 0.51652 │ 0.276169 │
│ 3 │ Param3 │ 0.496816 │ 0.275307 │

```

Apparently you can do this too

```julia
julia> sumstats[:Param2,:mean]|>first
0.5165204741265434

```