How to extract the Standard Deviation out of Turing.jl

Specific Domains Statistics
1

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][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.

# 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)")
2

Looks like the info on Chains is found somewhere else

https://github.com/TuringLang/MCMCChains.jl

# 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[1][][2][1][:mean][]
0.5556262828496785

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

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 │
2 Likes
4

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