using Random

# Load a plotting library.
using Plots

# Load Turing and MCMCChains.
using Turing, MCMCChains

# Load the distributions library.
using Distributions

# Load StatsPlots for density plots.
using StatsPlots

coin_flips = rand(Bernoulli(0.7), 100)

100-element Vector{Bool}:
 1
 0
 1
 1
 1
 1
 1
 1
 1
 0
 1
 1
 0
 ⋮
 0
 0
 1
 1
 0
 1
 1
 1
 1
 1
 1
 1
@model function coin(coin_flips)
    p ~ Beta(1, 1)
    for i ∈ 1:length(coin_flips)
        coin_flips[i] ~ Bernoulli(p)
    end
end;
  
#\in<tab> = ∈

chain_coin = sample(coin(coin_flips), MH(), 100) #single chain

Sampling: 100%|█████████████████████████████████████████| Time: 0:00:02

Chains MCMC chain (100×2×1 Array{Float64, 3}):

Iterations        = 1:1:100
Number of chains  = 1
Samples per chain = 100
Wall duration     = 6.07 seconds
Compute duration  = 6.07 seconds
parameters        = p
internals         = lp

Summary Statistics
  parameters      mean       std   naive_se      mcse       ess      rhat   es ⋯
      Symbol   Float64   Float64    Float64   Float64   Float64   Float64      ⋯

           p    0.7468    0.0597     0.0060    0.0091   29.1777    1.0600      ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           p    0.6117    0.7381    0.7623    0.7623    0.8284

Multi Parallel

MCMCThreads(): uses multithread stuff with Threads.jl

MCMCDistributed(): uses multiprocesses stuff with Distributed.jl and uses the MPI – Message Passing Interface

#We can Set enviornment variable "export JULIA_NUM_THREADS=4" in bashrc file to control number of threads
Threads.nthreads()

1
nchains = 3
chains = sample(coin(coin_flips), NUTS(), MCMCThreads(), 100, nchains)

┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC /home/patel_zeel/.julia/packages/AbstractMCMC/6aLyN/src/sample.jl:292
┌ Info: Found initial step size
│   ϵ = 0.8
└ @ Turing.Inference /home/patel_zeel/.julia/packages/Turing/rl6ku/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/patel_zeel/.julia/packages/AdvancedHMC/w90s5/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/patel_zeel/.julia/packages/AdvancedHMC/w90s5/src/hamiltonian.jl:47
┌ Info: Found initial step size
│   ϵ = 0.8
└ @ Turing.Inference /home/patel_zeel/.julia/packages/Turing/rl6ku/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/patel_zeel/.julia/packages/AdvancedHMC/w90s5/src/hamiltonian.jl:47
┌ Info: Found initial step size
│   ϵ = 0.8
└ @ Turing.Inference /home/patel_zeel/.julia/packages/Turing/rl6ku/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/patel_zeel/.julia/packages/AdvancedHMC/w90s5/src/hamiltonian.jl:47

Chains MCMC chain (100×13×3 Array{Float64, 3}):

Iterations        = 51:1:150
Number of chains  = 3
Samples per chain = 100
Wall duration     = 0.01 seconds
Compute duration  = 0.01 seconds
parameters        = p
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse       ess      rhat   es ⋯
      Symbol   Float64   Float64    Float64   Float64   Float64   Float64      ⋯

           p    0.6890    0.0480     0.0028    0.0048   91.5667    1.0340      ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           p    0.5949    0.6567    0.6887    0.7228    0.7723

Playing with MCMCChains object

Shape of chain: nsample $\times$ n_para_properties $\times$ nchains

  1. Access first chain
chains[:,:,1]

Chains MCMC chain (100×13×1 Array{Float64, 3}):

Iterations        = 51:1:150
Number of chains  = 1
Samples per chain = 100
Wall duration     = 0.01 seconds
Compute duration  = 0.01 seconds
parameters        = p
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse       ess      rhat   es ⋯
      Symbol   Float64   Float64    Float64   Float64   Float64   Float64      ⋯

           p    0.6833    0.0526     0.0053    0.0074   36.6413    1.0344      ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           p    0.5892    0.6503    0.6811    0.7185    0.7745

methodswith(MCMCChains.Chains) #to know methods name which we can apply on chain object
78-element Vector{Method}: