diff --git a/Project.toml b/Project.toml index 8c41eb8..1d7162d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MarkovChainHammer" uuid = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9" authors = ["Andre Souza "] -version = "0.0.8" +version = "0.0.9" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/BayesianMatrix/BayesianMatrix.jl b/src/BayesianMatrix/BayesianMatrix.jl index a1c419c..8216555 100644 --- a/src/BayesianMatrix/BayesianMatrix.jl +++ b/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") diff --git a/src/BayesianMatrix/constructors.jl b/src/BayesianMatrix/constructors.jl index b73ab82..9f25039 100644 --- a/src/BayesianMatrix/constructors.jl +++ b/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) \ No newline at end of file +""" + 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 diff --git a/src/BayesianMatrix/extensions.jl b/src/BayesianMatrix/extensions.jl index 2ab5c0b..1fd4e18 100644 --- a/src/BayesianMatrix/extensions.jl +++ b/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 diff --git a/test/test_BayesianMatrix.jl b/test/test_BayesianMatrix.jl index 76167f1..807aaba 100644 --- a/test/test_BayesianMatrix.jl +++ b/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,7 +83,6 @@ 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) @@ -92,3 +90,18 @@ end 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 \ No newline at end of file