# using Distributions, Bijectors

# Define a deterministic distribution
struct Determin{T<:Real} <: ContinuousUnivariateDistribution
  val::T
end

Distributions.rand(rng::AbstractRNG, d::Determin) = d.val
Distributions.logpdf(d::Determin, x::T) where T<:Real = zero(x)

Bijectors.bijector(d::Determin) = Identity{0}() # <= `0` indicates 0-dimensional, i.e. univariate 



# # useful functions
function compute_stats(chain::Chains, name::Symbol)
    # x = Array(x)
    # x = DataFrame(x)
    x = Array(group(chain, name))
    mean_x = mean(x, dims=1)
    xlb = mean_x - mapslices(x -> quantile(x, 0.025), x, dims=1)
    xub = mean_x - mapslices(x -> quantile(x, 0.975), x, dims=1)
    return (μ = mean_x, lb = xlb, ub = xub)
end

function compute_stats(x::Array)
    mean_x = mean(x, dims=1)
    xlb = mean_x - mapslices(x -> quantile(x, 0.025), x, dims=1)
    xub = mean_x - mapslices(x -> quantile(x, 0.975), x, dims=1)
    return (μ = mean_x, lb = xlb, ub = xub)
end

function get_omega(chain::Chains)   
    w1 = Array(group(chain, :ω1))
    w2 = Array(group(chain, :ω2))
    w = w1 .+ w2
    return (w = w, w1 = w1, w2 = w2)
end


function summarize_omega(chain::Chains)
    w1_array = Array(group(chain, :ω1))
    w2_array = Array(group(chain, :ω2))
    w_array = Array(group(chain, :ω))

    w = compute_stats(w_array)
    w1 = compute_stats(w1_array)
    w2 = compute_stats(w2_array)
    return (w = w, w1 = w1, w2 = w2)
end

invlogit(x) = 1/(1 + exp(-x))

function beta_params_mode(mode, k)
    #= http://doingbayesiandataanalysis.blogspot.com/2012/06/beta-distribution-parameterized-by-mode.html
    mode = (α - 1)/(α + β - 2)
    k = a + b # sample size
    =#
    @assert 0 < mode < 1
    @assert k > 2

    α = mode * (k-2) + 1
    β = (1-mode) * (k-2) + 1
    return α, β
end

function beta_params(μ, σ)
    @assert 0 < μ < 1
    @assert σ > 0

    α = ((1-μ)/σ^2 - (1/μ)) * μ^2
    β = α * (1 / μ - 1)
    return α, β
end

