Ce serveur Gitlab sera éteint le 30 juin 2020, pensez à migrer vos projets vers les serveurs gitlab-research.centralesupelec.fr et gitlab-student.centralesupelec.fr !

Commit cf317148 authored by Dahua Lin's avatar Dahua Lin

re-organized source files for conjugates

parent 6f2b8a2a
......@@ -279,15 +279,20 @@ include("kde.jl")
include("functionals.jl")
# Posteriors and conjugate priors
include("conjugates.jl")
include(joinpath("conjugate-normal", "normalgamma.jl"))
include(joinpath("conjugate-normal", "normalinversegamma.jl"))
include(joinpath("conjugate-normal", "normalwishart.jl"))
include(joinpath("conjugate-normal", "normalinversewishart.jl"))
include(joinpath("conjugate-normal", "normalknowncov.jl"))
include(joinpath("conjugates", "fallbacks.jl"))
include(joinpath("conjugates", "beta_binom.jl"))
include(joinpath("conjugates", "dirichlet_multi.jl"))
include(joinpath("conjugates", "gamma_exp.jl"))
include(joinpath("conjugates", "normalgamma.jl"))
include(joinpath("conjugates", "normalinversegamma.jl"))
include(joinpath("conjugates", "normalwishart.jl"))
include(joinpath("conjugates", "normalinversewishart.jl"))
include(joinpath("conjugates", "normal.jl"))
include(joinpath("conjugates", "mvnormal.jl"))
# other stuff
include("qq.jl")
include("estimators.jl")
end # module
# TODO: The non-conjugate normal prior on mean and inverse-gamma/Wishart prior
# on (co)variance.
#
# NormalInverseGammaSemiConjugate: p(\mu, \sigma) = N(\mu | m_0, V_0)IG(\sigma| a0, b0)
# NormalInverseWishartSemiConjugate: p(\mu,\Sigma) = N(\mu | m_0, V_0)IW(\Sigma|S_0,\nu_0)
# Conjugates for
#
# Beta - Bernoulli
# Beta - Binomial
#
posterior(prior::Beta, ss::BernoulliStats) = Beta(prior.alpha + ss.cnt1, prior.beta + ss.cnt0)
posterior(prior::Beta, ss::BinomialStats) = Beta(prior.alpha + ss.ns, prior.beta + (ss.ne * ss.n - ss.ns))
function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{T})
posterior(prior, suffstats(Binomial, n, x))
end
function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{T}, w::Array{Float64})
posterior(prior, suffstats(Binomial, n, x, w))
end
# Conjugates for
#
# Dirichlet - Categorical
# Dirichlet - Multinomial
#
function dirichlet_posupdate!{T<:Integer}(α::Vector{Float64}, ::Type{Categorical}, x::Array{T})
for i in 1:length(x)
@inbounds xi = x[i]
α[xi] += 1.0 # cannot put @inbounds here, as no guarantee xi is inbound
end
return α
end
function dirichlet_posupdate!{T<:Integer}(α::Vector{Float64}, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
if length(x) != length(w)
throw(ArgumentError("Inconsistent argument dimensions."))
end
for i in 1:length(x)
@inbounds xi = x[i]
@inbounds wi = w[i]
α[xi] += wi
end
return α
end
# each column are counts for one experiment
function dirichlet_posupdate!{T<:Real}(α::Vector{Float64}, ::Type{Multinomial}, X::Matrix{T})
d::Int = length(α)
if d != size(X, 1)
throw(ArgumentError("Inconsistent argument dimensions."))
end
n = size(X, 2)
o = 0
for j = 1:n
for i = 1:d
@inbounds α[i] += X[o+i]
end
o += d
end
return α
end
function dirichlet_posupdate!{T<:Real}(α::Vector{Float64}, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
d::Int = length(α)
if d != size(X, 1)
throw(ArgumentError("Inconsistent argument dimensions."))
end
n = size(X, 2)
if n != length(w)
throw(ArgumentError("Inconsistent argument dimensions."))
end
o = 0
for j = 1:n
@inbounds wj = w[j]
for i = 1:d
@inbounds α[i] += X[o+i] * wj
end
o += d
end
return α
end
function posterior{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Categorical, x))
end
function posterior{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Categorical, x, w))
end
function posterior_mode{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T})
α = dirichlet_posupdate!(copy(prior.alpha), Categorical, x)
dirichlet_mode!(α, α, sum(α))
end
function posterior_mode{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
α = dirichlet_posupdate!(copy(prior.alpha), Categorical, x, w)
dirichlet_mode!(α, α, sum(α))
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, x::Vector{T})
return Dirichlet(prior.alpha + x)
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Multinomial, X))
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Multinomial, X, w))
end
function posterior_mode{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, x::Vector{T})
α = prior.alpha + x
dirichlet_mode!(α, α, sum(α))
end
function posterior_mode{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T})
α = dirichlet_posupdate!(copy(prior.alpha), Multinomial, X)
dirichlet_mode!(α, α, sum(α))
end
function posterior_mode{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
α = dirichlet_posupdate!(copy(prior.alpha), Multinomial, X, w)
dirichlet_mode!(α, α, sum(α))
end
# Fallback functions for conjugates
posterior{D<:Distribution}(pri::Distribution, G::Type{D}, x) = posterior(pri, suffstats(G, x))
posterior(pri::Distribution, G::GenerativeFormulation, x) = posterior(pri, suffstats(G, x))
posterior{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = posterior(pri, suffstats(G, x, w))
posterior(pri::Distribution, G::GenerativeFormulation, x, w) = posterior(pri, suffstats(G, x, w))
posterior_rand(pri::Distribution, s::SufficientStats) = rand(posterior(pri, s))
posterior_rand{D<:Distribution}(pri::Distribution, G::Type{D}, x) = rand(posterior(pri, G, x))
posterior_rand{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = rand(posterior(pri, G, x, w))
posterior_rand(pri::Distribution, G::GenerativeFormulation, x) = rand(posterior(pri, G, x))
posterior_rand(pri::Distribution, G::GenerativeFormulation, x, w) = rand(posterior(pri, G, x, w))
posterior_rand!(r::Array, pri::Distribution, s::SufficientStats) = rand!(posterior(pri, s), r)
posterior_rand!{D<:Distribution}(r::Array, pri::Distribution, G::Type{D}, x) = rand!(posterior(pri, G, x), r)
posterior_rand!{D<:Distribution}(r::Array, pri::Distribution, G::Type{D}, x, w) = rand!(posterior(pri, G, x, w), r)
posterior_rand!(r::Array, pri::Distribution, G::GenerativeFormulation, x) = rand!(posterior(pri, G, x), r)
posterior_rand!(r::Array, pri::Distribution, G::GenerativeFormulation, x, w) = rand!(posterior(pri, G, x, w), r)
posterior_mode(pri::Distribution, s::SufficientStats) = mode(posterior(pri, s))
posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x) = mode(posterior(pri, G, x))
posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = mode(posterior(pri, G, x, w))
posterior_mode(pri::Distribution, G::GenerativeFormulation, x) = mode(posterior(pri, G, x))
posterior_mode(pri::Distribution, G::GenerativeFormulation, x, w) = mode(posterior(pri, G, x, w))
posterior_make{D<:Distribution}(::Type{D}, θ) = D(θ)
fit_map{D<:Distribution}(pri::Distribution, ::Type{D}, x) = posterior_make(D, posterior_mode(pri, D, x))
fit_map{D<:Distribution}(pri::Distribution, ::Type{D}, x, w) = posterior_make(D, posterior_mode(pri, D, x, w))
posterior_sample{D<:Distribution}(pri::Distribution, G::Type{D}, x) = D(rand(posterior(pri, G, x))...)
posterior_sample{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = D(rand(posterior(pri, G, x, w))...)
posterior_sample{D<:Distribution,G<:Distribution}(post::D, ::Type{G}) = G(rand(post)...)
# Conjuagtes for
#
# Gamma - Exponential
#
function posterior(prior::Gamma, ss::ExponentialStats)
return Gamma(prior.shape + ss.sw, 1.0 / (rate(prior) + ss.sx))
end
posterior_make(::Type{Exponential}, θ::Float64) = Exponential(1.0 / θ)
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference. Note that there were some typos in that document so the code
# here may not correspond exactly.
# Conjugates for Multivariate Normal
#
# MvNormal -- MvNormal (with known covariance)
#
#
# TODO: The non-conjugate normal prior on mean and inverse-gamma/Wishart prior
# on (co)variance.
#
import NumericExtensions.sum
......
### Beta -- Bernoulli or Binomial
posterior(prior::Beta, ss::BernoulliStats) = Beta(prior.alpha + ss.cnt1, prior.beta + ss.cnt0)
posterior(prior::Beta, ss::BinomialStats) = Beta(prior.alpha + ss.ns, prior.beta + (ss.ne * ss.n - ss.ns))
function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{T})
posterior(prior, suffstats(Binomial, n, x))
end
function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{T}, w::Array{Float64})
posterior(prior, suffstats(Binomial, n, x, w))
end
### Dirichlet -- Categorical or Multinomial
function dirichlet_posupdate!{T<:Integer}(α::Vector{Float64}, ::Type{Categorical}, x::Array{T})
for i in 1:length(x)
@inbounds xi = x[i]
α[xi] += 1.0 # cannot put @inbounds here, as no guarantee xi is inbound
end
return α
end
function dirichlet_posupdate!{T<:Integer}(α::Vector{Float64}, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
if length(x) != length(w)
throw(ArgumentError("Inconsistent argument dimensions."))
end
for i in 1:length(x)
@inbounds xi = x[i]
@inbounds wi = w[i]
α[xi] += wi
end
return α
end
# each column are counts for one experiment
function dirichlet_posupdate!{T<:Real}(α::Vector{Float64}, ::Type{Multinomial}, X::Matrix{T})
d::Int = length(α)
if d != size(X, 1)
throw(ArgumentError("Inconsistent argument dimensions."))
end
n = size(X, 2)
o = 0
for j = 1:n
for i = 1:d
@inbounds α[i] += X[o+i]
end
o += d
end
return α
end
function dirichlet_posupdate!{T<:Real}(α::Vector{Float64}, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
d::Int = length(α)
if d != size(X, 1)
throw(ArgumentError("Inconsistent argument dimensions."))
end
n = size(X, 2)
if n != length(w)
throw(ArgumentError("Inconsistent argument dimensions."))
end
o = 0
for j = 1:n
@inbounds wj = w[j]
for i = 1:d
@inbounds α[i] += X[o+i] * wj
end
o += d
end
return α
end
function posterior{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Categorical, x))
end
function posterior{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Categorical, x, w))
end
function posterior_mode{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T})
α = dirichlet_posupdate!(copy(prior.alpha), Categorical, x)
dirichlet_mode!(α, α, sum(α))
end
function posterior_mode{T<:Integer}(prior::Dirichlet, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
α = dirichlet_posupdate!(copy(prior.alpha), Categorical, x, w)
dirichlet_mode!(α, α, sum(α))
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, x::Vector{T})
return Dirichlet(prior.alpha + x)
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Multinomial, X))
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
Dirichlet(dirichlet_posupdate!(copy(prior.alpha), Multinomial, X, w))
end
function posterior_mode{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, x::Vector{T})
α = prior.alpha + x
dirichlet_mode!(α, α, sum(α))
end
function posterior_mode{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T})
α = dirichlet_posupdate!(copy(prior.alpha), Multinomial, X)
dirichlet_mode!(α, α, sum(α))
end
function posterior_mode{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
α = dirichlet_posupdate!(copy(prior.alpha), Multinomial, X, w)
dirichlet_mode!(α, α, sum(α))
end
### Gamma -- Exponential (rate)
function posterior(prior::Gamma, ss::ExponentialStats)
return Gamma(prior.shape + ss.sw, 1.0 / (rate(prior) + ss.sx))
end
posterior_make(::Type{Exponential}, θ::Float64) = Exponential(1.0 / θ)
### For Normal distributions
# Conjugates related to normal distribution
#
# Normal - Normal (with known sigma)
# InverseGamma - Normal (with known mu)
# NormalInverseGamma - Normal
#
# known sigma (prior on mu)
......@@ -232,6 +106,7 @@ end
# The NormalInverseGamma version is in normalinversegamma.jl
function posterior_sample{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T})
mu = rand(posterior(prior, suffstats(Normal, x)))
return Normal(mu, prior[2])
......@@ -261,3 +136,4 @@ function posterior_sample{T<:Real}(prior::(Float64, Gamma), ::Type{Normal}, x::A
tau2 = rand(posterior(prior, suffstats(Normal, x, w)))
return Normal(prior[1], 1./sqrt(tau2))
end
......@@ -261,36 +261,3 @@ fit_mle{D<:MultivariateDistribution}(dt::Type{D}, x::Matrix, w::Array) = fit_mle
fit{D<:Distribution}(dt::Type{D}, x) = fit_mle(D, x)
fit{D<:Distribution}(dt::Type{D}, args...) = fit_mle(D, args...)
# Conjugates
posterior{D<:Distribution}(pri::Distribution, G::Type{D}, x) = posterior(pri, suffstats(G, x))
posterior(pri::Distribution, G::GenerativeFormulation, x) = posterior(pri, suffstats(G, x))
posterior{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = posterior(pri, suffstats(G, x, w))
posterior(pri::Distribution, G::GenerativeFormulation, x, w) = posterior(pri, suffstats(G, x, w))
posterior_rand(pri::Distribution, s::SufficientStats) = rand(posterior(pri, s))
posterior_rand{D<:Distribution}(pri::Distribution, G::Type{D}, x) = rand(posterior(pri, G, x))
posterior_rand{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = rand(posterior(pri, G, x, w))
posterior_rand(pri::Distribution, G::GenerativeFormulation, x) = rand(posterior(pri, G, x))
posterior_rand(pri::Distribution, G::GenerativeFormulation, x, w) = rand(posterior(pri, G, x, w))
posterior_rand!(r::Array, pri::Distribution, s::SufficientStats) = rand!(posterior(pri, s), r)
posterior_rand!{D<:Distribution}(r::Array, pri::Distribution, G::Type{D}, x) = rand!(posterior(pri, G, x), r)
posterior_rand!{D<:Distribution}(r::Array, pri::Distribution, G::Type{D}, x, w) = rand!(posterior(pri, G, x, w), r)
posterior_rand!(r::Array, pri::Distribution, G::GenerativeFormulation, x) = rand!(posterior(pri, G, x), r)
posterior_rand!(r::Array, pri::Distribution, G::GenerativeFormulation, x, w) = rand!(posterior(pri, G, x, w), r)
posterior_mode(pri::Distribution, s::SufficientStats) = mode(posterior(pri, s))
posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x) = mode(posterior(pri, G, x))
posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = mode(posterior(pri, G, x, w))
posterior_mode(pri::Distribution, G::GenerativeFormulation, x) = mode(posterior(pri, G, x))
posterior_mode(pri::Distribution, G::GenerativeFormulation, x, w) = mode(posterior(pri, G, x, w))
posterior_make{D<:Distribution}(::Type{D}, θ) = D(θ)
fit_map{D<:Distribution}(pri::Distribution, ::Type{D}, x) = posterior_make(D, posterior_mode(pri, D, x))
fit_map{D<:Distribution}(pri::Distribution, ::Type{D}, x, w) = posterior_make(D, posterior_mode(pri, D, x, w))
posterior_sample{D<:Distribution}(pri::Distribution, G::Type{D}, x) = D(rand(posterior(pri, G, x))...)
posterior_sample{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = D(rand(posterior(pri, G, x, w))...)
posterior_sample{D<:Distribution,G<:Distribution}(post::D, ::Type{G}) = G(rand(post)...)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment