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 0bc04c86 authored by theodore's avatar theodore

Added multivariate t-distribution

parent ce14deb1
......@@ -72,6 +72,7 @@ export
MvNormal,
MvNormalCanon,
MvNormalKnownSigma,
MvTDist,
NegativeBinomial,
NoncentralBeta,
NoncentralChisq,
......@@ -275,6 +276,7 @@ include(joinpath("multivariate", "dirichlet.jl"))
include(joinpath("multivariate", "multinomial.jl"))
include(joinpath("multivariate", "mvnormal.jl"))
include(joinpath("multivariate", "mvnormalcanon.jl"))
include(joinpath("multivariate", "mvtdist.jl"))
include(joinpath("multivariate", "vonmisesfisher.jl"))
# Matrix distributions
......
# Multivariate t-distribution
## Generic multivariate t-distribution class
abstract AbstractMvTDist <: ContinuousMultivariateDistribution
immutable GenericMvTDist{Cov<:AbstractPDMat} <: AbstractMvTDist
df::Float64 # non-integer degrees of freedom allowed
dim::Int
zeromean::Bool
μ::Vector{Float64}
Σ::Cov
function GenericMvTDist{Cov}(df::Float64, dim::Int, zmean::Bool, μ::Vector{Float64}, Σ::Cov)
df > zero(df) || error("df must be positive")
new(float64(df), dim, zmean, μ, Σ)
end
end
function GenericMvTDist{Cov<:AbstractPDMat}(df::Float64, μ::Vector{Float64}, Σ::Cov, zmean::Bool)
d = length(μ)
dim(Σ) == d || throw(ArgumentError("The dimensions of μ and Σ are inconsistent."))
GenericMvTDist{Cov}(df, d, zmean, μ, Σ)
end
function GenericMvTDist{Cov<:AbstractPDMat}(df::Float64, μ::Vector{Float64}, Σ::Cov)
d = length(μ)
dim(Σ) == d || throw(ArgumentError("The dimensions of μ and Σ are inconsistent."))
GenericMvTDist{Cov}(df, d, allzeros(μ), μ, Σ)
end
function GenericMvTDist{Cov<:AbstractPDMat}(df::Float64, Σ::Cov)
d = dim(Σ)
GenericMvTDist{Cov}(df, d, true, zeros(d), Σ)
end
## Construction of multivariate normal with specific covariance type
typealias IsoTDist GenericMvTDist{ScalMat}
typealias DiagTDist GenericMvTDist{PDiagMat}
typealias MvTDist GenericMvTDist{PDMat}
MvTDist(df::Float64, μ::Vector{Float64}, C::PDMat) = GenericMvTDist(df, μ, C)
MvTDist(df::Float64, C::PDMat) = GenericMvTDist(df, C)
MvTDist(df::Float64, μ::Vector{Float64}, Σ::Matrix{Float64}) = GenericMvTDist(df, μ, PDMat(Σ))
MvTDist(df::Float64, Σ::Matrix{Float64}) = GenericMvTDist(df, PDMat(Σ))
DiagTDist(df::Float64, μ::Vector{Float64}, C::PDiagMat) = GenericMvTDist(df, μ, C)
DiagTDist(df::Float64, C::PDiagMat) = GenericMvTDist(df, C)
DiagTDist(df::Float64, μ::Vector{Float64}, σ::Vector{Float64}) = GenericMvTDist(df, μ, PDiagMat(abs2(σ)))
IsoTDist(df::Float64, μ::Vector{Float64}, C::ScalMat) = GenericMvTDist(df, μ, C)
IsoTDist(df::Float64, C::ScalMat) = GenericMvTDist(df, C)
IsoTDist(df::Float64, μ::Vector{Float64}, σ::Real) = GenericMvTDist(df, μ, ScalMat(length(μ), abs2(float64(σ))))
IsoTDist(df::Float64, d::Int, σ::Real) = GenericMvTDist(df, ScalMat(d, abs2(float64(σ))))
## convenient function to construct distributions of proper type based on arguments
gmvtdist(df::Float64, μ::Vector{Float64}, C::AbstractPDMat) = GenericMvTDist(df, μ, C)
gmvtdist(df::Float64, C::AbstractPDMat) = GenericMvTDist(df, C)
gmvtdist(df::Float64, μ::Vector{Float64}, σ::Real) = IsoTDist(df, μ, float64(σ))
gmvtdist(df::Float64, d::Int, σ::Float64) = IsoTDist(d, σ)
gmvtdist(df::Float64, μ::Vector{Float64}, σ::Vector{Float64}) = DiagTDist(df, μ, σ)
gmvtdist(df::Float64, μ::Vector{Float64}, Σ::Matrix{Float64}) = MvTDist(df, μ, Σ)
gmvtdist(df::Float64, Σ::Matrix{Float64}) = MvTDist(df, Σ)
# Basic statistics
dim(d::GenericMvTDist) = d.dim
mean(d::GenericMvTDist) = d.df>1 ? d.μ : NaN
mode(d::GenericMvTDist) = d.μ
modes(d::GenericMvTDist) = [mode(d)]
var(d::GenericMvTDist) = d.df>2 ? (d.df/(d.df-2))*diag(d.Σ) : Float64[NaN for i = 1:d.dim]
cov(d::GenericMvTDist) = d.df>2 ? (d.df/(d.df-2))*full(d.Σ) : NaN*ones(d.dim, d.dim)
invcov(d::GenericMvTDist) = d.df>2 ? ((d.df-2)/d.df)*full(inv(d.Σ)) : NaN*ones(d.dim, d.dim)
logdet_cov(d::GenericMvTDist) = d.df>2 ? logdet((d.df/(d.df-2))*d.Σ) : NaN
# For entropy calculations see "Multivariate t Distributions and their Applications", S. Kotz & S. Nadarajah
function entropy(d::GenericMvTDist)
hdf, hdim = 0.5*d.df, 0.5*d.dim
shdfhdim = hdf+hdim
0.5*logdet(d.Σ)+hdim*log(d.df*pi)+lbeta(hdim, hdf)-lgamma(hdim)+shdfhdim*(digamma(shdfhdim)-digamma(hdf))
end
# evaluation (for GenericMvTDist)
function sqmahal(d::GenericMvTDist, x::Vector{Float64})
z::Vector{Float64} = d.zeromean ? x : x - d.μ
invquad(d.Σ, z)
end
function sqmahal!(r::Array{Float64}, d::GenericMvTDist, x::Matrix{Float64})
if !(size(x, 1) == dim(d) && size(x, 2) == length(r))
throw(ArgumentError("Inconsistent argument dimensions."))
end
z::Matrix{Float64} = d.zeromean ? x : bsubtract(x, d.μ, 1)
invquad!(r, d.Σ, z)
end
# generic PDF evaluation (appliable to AbstractMvTDist)
insupport{T<:Real}(d::AbstractMvTDist, x::Vector{T}) = dim(d) == length(x) && allfinite(x)
insupport{T<:Real}(d::AbstractMvTDist, x::Matrix{T}) = dim(d) == size(x, 1) && allfinite(x)
insupport{G<:AbstractMvTDist,T<:Real}(::Type{G}, x::Vector{T}) = allfinite(x)
insupport{G<:AbstractMvTDist,T<:Real}(::Type{G}, x::Matrix{T}) = allfinite(x)
sqmahal(d::AbstractMvTDist, x::Matrix{Float64}) = sqmahal!(Array(Float64, size(x, 2)), d, x)
function logpdf(d::AbstractMvTDist, x::Vector{Float64})
hdf, hdim = 0.5*d.df, 0.5*d.dim
shdfhdim = hdf+hdim
lgamma(shdfhdim)-lgamma(hdf)-hdim*log(d.df)-hdim*log(pi)-0.5*logdet(d.Σ)-shdfhdim*log(1+sqmahal(d, x)/d.df)
end
function logpdf!(r::Array{Float64}, d::AbstractMvTDist, x::Matrix{Float64})
sqmahal!(r, d, x)
hdf, hdim = 0.5*d.df, 0.5*d.dim
shdfhdim = hdf+hdim
for i = 1:size(x, 2)
r[i] = lgamma(shdfhdim)-lgamma(hdf)-hdim*log(d.df)-hdim*log(pi)-0.5*logdet(d.Σ)-shdfhdim*log(1+r[i]/d.df)
end
r
end
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