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 469779f9 authored by Sergey Bartunov's avatar Sergey Bartunov

entropy function for Wishart distribution

parent 1115ff27
......@@ -165,7 +165,8 @@ export
std, # standard deviation of distribution
suffstats, # compute sufficient statistics
var, # variance of distribution
wsample # weighted sampling from a source array
wsample, # weighted sampling from a source array
expected_log_det # expected logarithm of random matrix determinant
import Base.mean, Base.median, Base.quantile, Base.scale
import Base.max, Base.min, Base.maximum, Base.minimum
......
......@@ -36,15 +36,31 @@ mean(w::Wishart) = w.nu * w.Schol[:U]' * w.Schol[:U]
pdf(W::Wishart, X::Matrix{Float64}) = exp(logpdf(W, X))
p(W::Wishart) = size(W.Schol, 1)
function expected_log_det(W::Wishart)
logd = 0.
for i=1:p(W)
logd += digamma(0.5 * (W.nu + 1 - i))
end
logd += p(W) * log(2)
logd += log(det(W.Schol))
return logd
end
function log_norm(W::Wishart)
return (W.nu / 2) * log(det(W.Schol)) + (p(W) * W.nu / 2) * log(2) + lpgamma(p(W), W.nu / 2)
end
function logpdf(W::Wishart, X::Matrix{Float64})
if !insupport(W, X)
return -Inf
else
p = size(X, 1)
logd::Float64 = W.nu * p / 2.0 * log(2.0)
logd += W.nu / 2.0 * log(det(W.Schol))
logd += lpgamma(p, W.nu / 2.0)
logd = -logd
logd = -log_norm(W)
logd += 0.5 * (W.nu - p - 1.0) * logdet(X)
logd -= 0.5 * trace(W.Schol \ X)
return logd
......@@ -68,4 +84,8 @@ function rand(w::Wishart)
return Z' * Z
end
function entropy(W::Wishart)
return log_norm(W) - (W.nu - p(W) - 1) / 2 * expected_log_det(W) + W.nu * p(W) / 2
end
var(w::Wishart) = error("Not yet implemented")
# Tests on Wishart distributions
require("../src/Distributions.jl")
using Distributions
using Base.Test
V = [[2. 1.], [1. 2.]]
W = Wishart(3., V)
# logdet
@test_approx_eq expected_log_det(W) 1.9441809588650447
# entropy
@test_approx_eq entropy(W) 7.178942679971454
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