Merge pull request #56 from sandreza/save_param

params() functionality

Project.toml

@@ -1,7 +1,7 @@
 name = "MarkovChainHammer"
 uuid = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9"
 authors = ["Andre Souza <[email protected]>"]
-version = "0.0.8"
+version = "0.0.9"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/BayesianMatrix/BayesianMatrix.jl

@@ -8,12 +8,12 @@ import MarkovChainHammer.TransitionMatrix: holding_times, count_operator
 # export functions
 export BayesianGenerator, GeneratorParameterDistributions
 export eigen_distribution, eigvals_distribution
-export uninformative_prior
+export uninformative_prior 
+export params, unpack
 
 # general abstractions
-struct BayesianGenerator{PB,D,PA,PP}
+struct BayesianGenerator{PB,PA,PP}
     prior::PB
-    data::D
     posterior::PA
     predictive::PP
 end
@@ -29,7 +29,6 @@ struct GeneratorPredictiveDistributions{H,P}
 end
 
 include("constructors.jl")
-export uninformative_prior
 
 include("extensions.jl")
 

src/BayesianMatrix/constructors.jl

@@ -48,19 +48,52 @@ The covariance is CoVar[X⃗ ⊗ X⃗] = Diagonal(α̃) - α̃ ⊗ α̃ / (α₀
 The variance Var(Xᵢ) = α̃ᵢ(1-α̃ᵢ) / (α₀ + 1) where α₀ = sum(αs).
 
 """
-function GeneratorParameterDistributions(number_of_states::Int; α=1, β=1, αs=ones(number_of_states - 1))
-    @assert length(αs) == number_of_states - 1
+function GeneratorParameterDistributions(number_of_states::Int;  α=ones(number_of_states), β=ones(number_of_states), αs=ones(number_of_states - 1, number_of_states))
+    if length(α) == 1 && length(β) == 1 
+        α = fill(α, number_of_states)
+        β = fill(β, number_of_states)
+        αs = repeat(αs, 1, number_of_states)
+    end
+    @assert size(αs)[1] == number_of_states - 1
+    @assert size(αs)[2] == size(α)[1] == size(β)[1] == number_of_states
     @assert all(αs .> 0)
-    @assert α > 0
-    @assert β > 0
+    @assert all(α .> 0)
+    @assert all(β .> 0)
     α = α
     β = β
-    θ = 1 / β
-    rates = [Gamma(α, θ) for i in 1:number_of_states]
-    exit_probabilities = [Dirichlet(αs) for i in 1:number_of_states]
+    θ = 1 ./ β
+    rates = [Gamma(α[i], θ[i]) for i in 1:number_of_states]
+    exit_probabilities = [Dirichlet(αs[:,i]) for i in 1:number_of_states]
     return GeneratorParameterDistributions(rates, exit_probabilities)
 end
 
+function GeneratorParameterDistributions(parameters::NamedTuple)
+    @assert haskey(parameters, :α)
+    @assert haskey(parameters, :β)
+    @assert haskey(parameters, :αs)
+    @assert length(parameters.α) == length(parameters.β) == length(parameters.αs[1,:]) == (length(parameters.αs[:,1])+1)
+    return GeneratorParameterDistributions(length(parameters.α), α=parameters.α, β=parameters.β, αs=parameters.αs)
+end
+
+function GeneratorPredictiveDistributions(number_of_states::Int; μ=ones(number_of_states), σ=ones(number_of_states), ξ=ones(number_of_states), n=ones(number_of_states), αs=ones(number_of_states - 1, number_of_states))
+    if length(μ) == 1 && length(σ) == 1
+        μ = fill(α, number_of_states)
+        σ = fill(β, number_of_states)
+        ξ = fill(ξ, number_of_states)
+        n = fill(n, number_of_states)
+        αs = repeat(αs, 1, number_of_states)
+    end
+    @assert size(αs)[1] == number_of_states - 1
+    @assert size(αs)[2] == size(μ)[1] == size(ξ)[1] == size(n)[1] == number_of_states
+    @assert all(αs .> 0)
+    @assert all(n .> 0)
+    @assert all(σ .> 0)
+
+    holding_times = [GeneralizedPareto(μ[i], σ[i], ξ[i]) for i in 1:number_of_states]
+    exit_counts = [DirichletMultinomial(n[i], αs[:,i]) for i in 1:number_of_states]
+    return GeneratorPredictiveDistributions(holding_times, exit_counts)
+end
+
 """
     BayesianGenerator(data, prior::GeneratorParameterDistributions; dt=1)
 
@@ -124,10 +157,10 @@ function BayesianGenerator(data, prior::GeneratorParameterDistributions; dt=1)
     end
     posterior = GeneratorParameterDistributions(posterior_rates, posterior_exit_probabilities)
     predictive = GeneratorPredictiveDistributions(predictive_holding_times, predictive_exit_counts)
-    return BayesianGenerator(prior, data, posterior, predictive)
+    return BayesianGenerator(prior, posterior, predictive)
 end
 
-BayesianGenerator(prior::GeneratorParameterDistributions) = BayesianGenerator(prior, nothing, prior, nothing)
+BayesianGenerator(prior::GeneratorParameterDistributions) = BayesianGenerator(prior, prior, nothing)
 BayesianGenerator(data; dt = 1) = BayesianGenerator(data,  uninformative_prior(maximum(data)) ; dt=dt)
 
 """
@@ -157,5 +190,63 @@ function BayesianGenerator(data::Vector{Vector{Int64}}, prior::GeneratorParamete
 end
 BayesianGenerator(data::Vector{Vector{Int64}}; dt=1) = BayesianGenerator(data, uninformative_prior(maximum(reduce(vcat, data))); dt=dt)
 
+"""
+    BayesianGenerator(parameters::NamedTuple)
+
+Construct a BayesianGenerator object from a NamedTuple of parameters.
+
+# Arguments
+- `parameters::NamedTuple`: A NamedTuple containing the parameters for the BayesianGenerator object. Must contain the fields `prior`, `posterior` and `predictive`, each of which must be a NamedTuple containing the parameters for the respective distributions.
+
+# Returns
+- `BayesianGenerator`: A BayesianGenerator object constructed from the parameters. Contains the posterior distributions for the rates and exit probabilities, as well as the predictive distributions for the holding times and exit counts.
+
+"""
+function BayesianGenerator(parameters::NamedTuple) 
+    @assert haskey(parameters, :prior)
+    @assert haskey(parameters, :posterior)
+    @assert haskey(parameters, :predictive)
+
+    number_of_states = length(parameters.prior.α)
+    prior = GeneratorParameterDistributions(number_of_states; α=parameters.prior.α, β=parameters.prior.β, αs=parameters.prior.αs)
+    posterior = GeneratorParameterDistributions(number_of_states; α=parameters.posterior.α, β=parameters.posterior.β, αs=parameters.posterior.αs)
+    predictive = GeneratorPredictiveDistributions(number_of_states; μ=parameters.predictive.μ, σ=parameters.predictive.σ, ξ=parameters.predictive.ξ, n=parameters.predictive.n, αs=parameters.predictive.αs)
+    return BayesianGenerator(prior, posterior, predictive)
+end
+
+"""
+    uninformative_prior(number_of_states; scale=eps(1e2))
+
+Construct an uninformative prior distribution for a BayesianGenerator object.
+
+# Arguments
+- `number_of_states::Int`: The number of states in the BayesianGenerator object.
+
+# Keyword Arguments
+
+- `scale::Number=eps(1e2)`: The scale parameter for the Gamma and Dirichlet distributions.
+
+# Returns
+- `GeneratorParameterDistributions`: An uninformative prior distribution for a BayesianGenerator object.
+
+"""
+uninformative_prior(number_of_states; scale=eps(1e2)) = GeneratorParameterDistributions(number_of_states::Int; α=scale, β=scale, αs=ones(number_of_states - 1) * scale)
 
-uninformative_prior(number_of_states; scale=eps(1e2)) = GeneratorParameterDistributions(number_of_states::Int; α=scale, β=scale, αs=ones(number_of_states - 1) * scale)
+"""
+    unpack(Q::BayesianGenerator)
+
+Unpack a BayesianGenerator object into its parameters.
+
+# Arguments
+- `Q::BayesianGenerator`: The BayesianGenerator object to unpack.
+
+# Returns
+- `Tuple{NamedTuple{(:prior, :posterior, :predictive),Tuple{NamedTuple{(:α, :β, :αs),Tuple{Vector{Float64},Vector{Float64},Vector{Float64}}},NamedTuple{(:α, :β, :αs),Tuple{Vector{Float64},Vector{Float64},Vector{Float64}}},NamedTuple{(:μ, :σ, :ξ, :n, :αs),Tuple{Vector{Float64},Vector{Float64},Vector{Float64},Vector{Int64},Vector{Vector{Float64}}}}}}}`: A tuple containing the parameters for the prior, posterior and predictive distributions.
+
+"""
+function unpack(Q::BayesianGenerator)
+    prior = params(Q.prior)
+    posterior = params(Q.posterior)
+    predictive = params(Q.predictive)
+    return (; prior, posterior, predictive)
+end

src/BayesianMatrix/extensions.jl

@@ -12,6 +12,36 @@ rand(Q::BayesianGenerator) = rand(Q.posterior)
 rand(Q::GeneratorParameterDistributions, n::Int) = [rand(Q) for i in 1:n]
 rand(Q::BayesianGenerator, n::Int) = [rand(Q) for i in 1:n]
 
+# extend Distributions 
+import Distributions.params
+
+function params(dist::GeneratorParameterDistributions)
+    number_of_states = size(dist.rates)[1]
+    α = zeros(number_of_states)
+    β = zeros(number_of_states)
+    αs = zeros(number_of_states-1, number_of_states)
+    for i in 1:number_of_states
+        α[i] = params(dist.rates[i])[1]
+        β[i] = 1/params(dist.rates[i])[2]
+        αs[:,i] = params(dist.exit_probabilities[i])[1]
+    end
+    return (; α, β, αs)
+end
+
+function params(dist::GeneratorPredictiveDistributions)
+    number_of_states = size(dist.holding_times)[1]
+    μ=zeros(number_of_states)
+    σ=zeros(number_of_states)
+    ξ=zeros(number_of_states)
+    n=zeros(Int64, number_of_states)
+    αs=zeros(number_of_states - 1, number_of_states)
+    for i in 1:number_of_states
+        μ[i], σ[i], ξ[i] = params(dist.holding_times[i])
+        n[i], αs[:,i] = params(dist.exit_counts[i])
+    end
+    return (; μ, σ, ξ, n, αs)
+end
+
 # Base.Statistics 
 import Statistics.mean, Statistics.var, Statistics.std
 mean(Q::BayesianGenerator) = mean(Q.posterior)
@@ -59,7 +89,7 @@ eigvals_distribution(Q::BayesianGenerator; samples=100) = eigvals_distribution(Q
 import Base.copy
 copy(Q::GeneratorParameterDistributions) = GeneratorParameterDistributions(copy(Q.rates), copy(Q.exit_probabilities))
 copy(Q::GeneratorPredictiveDistributions) = GeneratorPredictiveDistributions(copy(Q.holding_times), copy(Q.exit_counts))
-copy(Q::BayesianGenerator) = BayesianGenerator(copy(Q.prior), copy(Q.data), copy(Q.posterior), copy(Q.predictive))
+copy(Q::BayesianGenerator) = BayesianGenerator(copy(Q.prior), copy(Q.posterior), copy(Q.predictive))
 
 # Base.show
 import Base.show

test/test_BayesianMatrix.jl

@@ -60,7 +60,6 @@ end
     end
 end
 
-
 @testset "Bayesian Matrix: Prior to Posterior" begin
     # test confidence with respect to same data
     markov_chain = [1 1 1 1 2 2 1 1 1 3 1 1]
@@ -84,11 +83,25 @@ end
     @test sum([sum(abs.((params.(Q21.posterior.rates)[i].-params.(Q12.posterior.rates)[i])[1])) for i in 1:3]) < eps(10.0)
 end
 
-
 @testset "Bayesian Matrix: Ensemble Functionality" begin
     data = [rand([1, 2, 3], 10) for i in 1:2]
     Q = BayesianGenerator(data)
     Q1 = BayesianGenerator(data[1])
     Q2 = BayesianGenerator(data[2], Q1.posterior)
     @test all( abs.(mean(Q) - mean(Q2)) .< eps(10.0))
 end
+
+@testset "Unpacking and Reconstructing" begin
+    # test confidence with respect to same data
+    markov_chain = [1 1 1 1 2 2 1 1 1 3 1 1]
+    Q1 = BayesianGenerator(markov_chain)
+    parameters = unpack(Q1)
+    Q2 = BayesianGenerator(parameters)
+    @test all(mean(Q1) .== mean(Q2))
+    @test all(var(Q1) .== var(Q2))
+    parameters_posterior = params(Q1.posterior)
+    parameters_posterior_2 = params(GeneratorParameterDistributions(parameters_posterior))
+    @test all(parameters_posterior.α .== parameters_posterior_2.α)
+    @test all(parameters_posterior.β .== parameters_posterior_2.β)
+    @test all(parameters_posterior.αs .== parameters_posterior_2.αs)
+end