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

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/

I get the error

``````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

``````# personally extracting the statistics for the parameter p
p_mean = mean(chain[:p])[:mean] # Extracting the mean for parameter p
p_std  =  std(chain[:p])[:std]  # 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.

``````# 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] # Extracting the mean for parameter p
p_std  =  std(chain[:p])[:std]  # Extracting the  std for parameter p

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

``````

Looks like the info on Chains is found somewhere else

``````# 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> p_mean = chain.info[][:mean][]
0.5556262828496785

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

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:

``````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.

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

`sumstats` looks like this:

``````│ 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 │

``````
1 Like

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.

``````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:

``````│ 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> sumstats[:Param2,:mean]|>first
0.5165204741265434
``````