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 1beadf17 authored by Dahua Lin's avatar Dahua Lin

Implement posterior_mode and fit_map for a series of distributions.

parent 1056c476
......@@ -114,7 +114,11 @@ export
logpdf!, # evaluate log pdf to provided storage
logpmf, # log probability mass
logpmf!, # evaluate log pmf to provided storage
posterior, # Bayesian updating
posterior, # get posterior distribution given prior and observed data
posterior_mode, # get the mode of posterior distribution
posterior_rand, # draw samples from the posterior distribution
posterior_rand!,
posterior_make, # create a distribution/model from params obtained from posterior
scale, # scale parameter of a distribution
rate, # rate parameter of a distribution
sqmahal, # squared Mahalanobis distance to Gaussian center
......
......@@ -13,54 +13,48 @@ function posterior{T<:Real}(prior::Beta, ::Type{Binomial}, n::Integer, x::Array{
posterior(prior, suffstats(Binomial, n, x, w))
end
### Dirichlet -- Categorical or Multinomial
function posterior{T<:Real}(prior::Dirichlet, ::Type{Categorical}, x::Array{T})
α = copy(prior.alpha)
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 Dirichlet(α)
return α
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Categorical}, x::Array{T}, w::Array{Float64})
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
α = copy(prior.alpha)
for i in 1:length(x)
@inbounds xi = x[i]
@inbounds wi = w[i]
α[xi] += wi
end
return Dirichlet(α)
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, x::Vector{T})
return Dirichlet(prior.alpha + x)
return α
end
# each column are counts for one experiment
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T})
d::Int = dim(prior)
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)
α = copy(prior.alpha)
o = 0
for j in 1:n
for i in 1:d
for j = 1:n
for i = 1:d
@inbounds α[i] += X[o+i]
end
o += d
end
return Dirichlet(α)
return α
end
function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T}, w::Array{Float64})
d::Int = dim(prior)
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
......@@ -68,29 +62,77 @@ function posterior{T<:Real}(prior::Dirichlet, ::Type{Multinomial}, X::Matrix{T},
if n != length(w)
throw(ArgumentError("Inconsistent argument dimensions."))
end
α = copy(prior.alpha)
o = 0
for j in 1:n
for j = 1:n
@inbounds wj = w[j]
for i in 1:d
for i = 1:d
@inbounds α[i] += X[o+i] * wj
end
o += d
end
return Dirichlet(α)
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
### 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
# known sigma (prior on mu)
function posterior(prior::Normal, ss::NormalKnownSigmaStats)
μ0 = prior.μ
c0 = 1.0 / abs2(prior.σ)
......@@ -102,24 +144,54 @@ function posterior(prior::Normal, ss::NormalKnownSigmaStats)
return Normal(μ1, σ1)
end
function posterior(prior::InverseGamma, ss::NormalKnownMuStats)
α1 = prior.shape + ss.tw / 2
β1 = rate(prior) + ss.s2 / 2
return InverseGamma(α1, 1.0 / β1)
end
function posterior{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T})
pri_μ::Normal = prior[1]
σ::Float64 = prior[2]
pri_μ::Normal, σ::Float64 = prior
posterior(pri_μ, suffstats(NormalKnownSigma(σ), x))
end
function posterior{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T}, w::Array{Float64})
pri_μ::Normal = prior[1]
σ::Float64 = prior[2]
pri_μ::Normal, σ::Float64 = prior
posterior(pri_μ, suffstats(NormalKnownSigma(σ), x, w))
end
function posterior_mode(prior::Normal, ss::NormalKnownSigmaStats)
μ0 = prior.μ
c0 = 1.0 / abs2(prior.σ)
c1 = 1.0 / abs2(ss.σ)
return (μ0 * c0 + ss.s * c1) / (c0 + ss.tw * c1)
end
function posterior_mode{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T})
pri_μ::Normal, σ::Float64 = prior
posterior_mode(pri_μ, suffstats(NormalKnownSigma(σ), x))
end
function posterior_mode{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T}, w::Array{Float64})
pri_μ::Normal, σ::Float64 = prior
posterior_mode(pri_μ, suffstats(NormalKnownSigma(σ), x, w))
end
function fit_map{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T})
pri_μ::Normal, σ::Float64 = prior
μ = posterior_mode(pri_μ, suffstats(NormalKnownSigma(σ), x))
Normal(μ, σ)
end
function fit_map{T<:Real}(prior::(Normal, Float64), ::Type{Normal}, x::Array{T}, w::Array{Float64})
pri_μ::Normal, σ::Float64 = prior
μ = posterior_mode(pri_μ, suffstats(NormalKnownSigma(σ), x, w))
Normal(μ, σ)
end
# known mu (prior on sigma)
function posterior(prior::InverseGamma, ss::NormalKnownMuStats)
α1 = prior.shape + ss.tw / 2
β1 = rate(prior) + ss.s2 / 2
return InverseGamma(α1, 1.0 / β1)
end
function posterior{T<:Real}(prior::(Float64, InverseGamma), ::Type{Normal}, x::Array{T})
μ::Float64 = prior[1]
pri_σ::InverseGamma = prior[2]
......
# uniform interface for model estimation
export Estimator, MLEstimator
export Estimator, MLEstimator, MAPEstimator
export nsamples, estimate, prior_score
abstract Estimator{D<:Distribution}
......@@ -16,5 +16,10 @@ estimate{D<:Distribution}(e::MLEstimator{D}, x, w) = fit_mle(D, x, w)
prior_score{D<:Distribution}(e::MLEstimator{D}, d::D) = 0.
# TODO: add MAPEstimator
immutable MAPEstimator{D<:Distribution,Pri} <: Estimator{D}
pri::Pri
end
MAPEstimator{D<:Distribution,Pri}(::Type{D}, pri::Pri) = MAPEstimator{D,Pri}(pri)
estimate{D<:Distribution,Pri}(e::MAPEstimator{D,Pri}, x) = fit_map(e.pri, D, x)
estimate{D<:Distribution,Pri}(e::MAPEstimator{D,Pri}, x, w) = fit_map(e.pri, D, x, w)
......@@ -287,3 +287,7 @@ posterior_mode{D<:Distribution}(pri::Distribution, G::Type{D}, x, w) = mode(post
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))
......@@ -80,21 +80,21 @@ function entropy(d::Dirichlet)
return en
end
function mode(d::Dirichlet)
k = length(d.alpha)
x = Array(Float64, k)
s = d.alpha0 - k
a = d.alpha
for i in 1:k
@inbounds ai = a[i]
if ai <= 1.0
error("Dirichlet has a mode only when alpha[i] > 1 for all i")
function dirichlet_mode!(r::Vector{Float64}, α::Vector{Float64}, α0::Float64)
k = length(α)
s = α0 - k
for i = 1:k
@inbounds ai = α[i]
if ai <= 1.
error("Dirichlet has a mode only when alpha[i] > 1 for all i" )
end
x[i] = (ai - 1.0) / s
@inbounds r[i] = (ai - 1.0) / s
end
return x
r
end
mode(d::Dirichlet) = dirichlet_mode!(zeros(dim(d)), d.alpha, d.alpha0)
modes(d::Dirichlet) = [mode(d)]
......
......@@ -33,11 +33,19 @@ p = posterior(pri, Bernoulli, x)
@test_approx_eq p.alpha pri.alpha + sum(x)
@test_approx_eq p.beta pri.beta + (n - sum(x))
f = fit_map(pri, Bernoulli, x)
@test isa(f, Bernoulli)
@test_approx_eq f.p1 mode(p)
p = posterior(pri, Bernoulli, x, w)
@test isa(p, Beta)
@test_approx_eq p.alpha pri.alpha + sum(x .* w)
@test_approx_eq p.beta pri.beta + (sum(w) - sum(x .* w))
f = fit_map(pri, Bernoulli, x, w)
@test isa(f, Bernoulli)
@test_approx_eq f.p1 mode(p)
x = rand(Binomial(10, 0.3), n)
p = posterior(pri, Binomial, 10, x)
@test isa(p, Beta)
......@@ -58,24 +66,42 @@ p = posterior(pri, Categorical, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + ccount(3, x)
f = fit_map(pri, Categorical, x)
@test isa(f, Categorical)
@test_approx_eq f.prob mode(p)
p = posterior(pri, Categorical, x, w)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + ccount(3, x, w)
f = fit_map(pri, Categorical, x, w)
@test isa(f, Categorical)
@test_approx_eq f.prob mode(p)
x = rand(Multinomial(100, [0.2, 0.3, 0.5]), 1)
p = posterior(pri, Multinomial, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + x
r = posterior_mode(pri, Multinomial, x)
@test_approx_eq r mode(p)
x = rand(Multinomial(10, [0.2, 0.3, 0.5]), n)
p = posterior(pri, Multinomial, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + vec(sum(x, 2))
r = posterior_mode(pri, Multinomial, x)
@test_approx_eq r mode(p)
p = posterior(pri, Multinomial, x, w)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + vec(x * w)
r = posterior_mode(pri, Multinomial, x, w)
@test_approx_eq r mode(p)
# Gamma - Exponential
......@@ -87,14 +113,24 @@ p = posterior(pri, Exponential, x)
@test_approx_eq p.shape pri.shape + n
@test_approx_eq rate(p) rate(pri) + sum(x)
f = fit_map(pri, Exponential, x)
@test isa(f, Exponential)
@test_approx_eq rate(f) mode(p)
p = posterior(pri, Exponential, x, w)
@test isa(p, Gamma)
@test_approx_eq p.shape pri.shape + sum(w)
@test_approx_eq rate(p) rate(pri) + sum(x .* w)
f = fit_map(pri, Exponential, x, w)
@test isa(f, Exponential)
@test_approx_eq rate(f) mode(p)
# Normal likelihood
# known sigma
pri = Normal(1.0, 5.0)
x = rand(Normal(2.0, 3.0), n)
......@@ -103,11 +139,29 @@ p = posterior((pri, 3.0), Normal, x)
@test_approx_eq mean(p) (mean(pri) / var(pri) + sum(x) / 9.0) / (1.0 / var(pri) + n / 9.0)
@test_approx_eq var(p) inv(1.0 / var(pri) + n / 9.0)
r = posterior_mode((pri, 3.0), Normal, x)
@test_approx_eq r mode(p)
f = fit_map((pri, 3.0), Normal, x)
@test isa(f, Normal)
@test f.μ == r
@test f.σ == 3.0
p = posterior((pri, 3.0), Normal, x, w)
@test isa(p, Normal)
@test_approx_eq mean(p) (mean(pri) / var(pri) + dot(x, w) / 9.0) / (1.0 / var(pri) + sum(w) / 9.0)
@test_approx_eq var(p) inv(1.0 / var(pri) + sum(w) / 9.0)
r = posterior_mode((pri, 3.0), Normal, x, w)
@test_approx_eq r mode(p)
f = fit_map((pri, 3.0), Normal, x, w)
@test isa(f, Normal)
@test f.μ == r
@test f.σ == 3.0
# known mu
pri = InverseGamma(1.5, 0.5) # β = 2.0
x = rand(Normal(2.0, 3.0), n)
......@@ -120,3 +174,5 @@ p = posterior((2.0, pri), Normal, x, w)
@test isa(p, InverseGamma)
@test_approx_eq p.shape pri.shape + sum(w) / 2
@test_approx_eq rate(p) rate(pri) + dot(w, abs2(x - 2.0)) / 2
using Distributions
using Base.Test
# Fit MLE
N = 10^5
d = fit(DiscreteUniform, rand(DiscreteUniform(10, 15), N))
......@@ -72,3 +74,4 @@ d = fit(Laplace, rand(Laplace(5.0, 3.0), N))
d = fit(Poisson, rand(Poisson(8.2), N))
@test isa(d, Poisson)
@test_approx_eq_eps mean(d) 8.2 0.2
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