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 00dd06fb authored by John Myles White's avatar John Myles White

Fill in missing methods

Add entropy(), mgf(), cf(), etc.
parent 05539ef1
......@@ -158,10 +158,10 @@ include(joinpath("univariate", "cauchy.jl"))
include(joinpath("univariate", "chi.jl"))
include(joinpath("univariate", "chisq.jl"))
include(joinpath("univariate", "discreteuniform.jl"))
include(joinpath("univariate", "erlang.jl"))
include(joinpath("univariate", "exponential.jl"))
include(joinpath("univariate", "fdist.jl"))
include(joinpath("univariate", "gamma.jl"))
include(joinpath("univariate", "erlang.jl"))
include(joinpath("univariate", "geometric.jl"))
include(joinpath("univariate", "gumbel.jl"))
include(joinpath("univariate", "hypergeometric.jl"))
......
using Distributions
abstract TruncatedContinuousUnivariateDistribution <: ContinuousUnivariateDistribution
abstract TruncatedDiscreteUnivariateDistribution <: DiscreteUnivariateDistribution
typealias TruncatedUnivariateDistribution Union(TruncatedContinuousUnivariateDistribution, TruncatedDiscreteUnivariateDistribution)
......
......@@ -19,6 +19,8 @@ function cdf(d::Arcsine, x::Number)
end
end
entropy(d::Arcsine) = -log(2.0) / pi
function insupport(d::Arcsine, x::Number)
if -1.0 <= x <= 1.0
return true
......
......@@ -17,9 +17,9 @@ Beta() = Beta(1.0) # uniform
function entropy(d::Beta)
o = lbeta(d.alpha, d.beta)
o = o - (d.alpha - 1.0) * digamma(d.alpha)
o = o - (d.beta - 1.0) * digamma(d.beta)
o = o + (d.alpha + d.beta - 2.0) * digamma(d.alpha + d.beta)
o -= (d.alpha - 1.0) * digamma(d.alpha)
o -= (d.beta - 1.0) * digamma(d.beta)
o += (d.alpha + d.beta - 2.0) * digamma(d.alpha + d.beta)
return o
end
......@@ -47,13 +47,16 @@ function rand(d::Beta)
return u / (u + rand(Gamma(d.beta)))
end
# TODO: Don't create temporaries here
function rand(d::Beta, dims::Dims)
u = rand(Gamma(d.alpha), dims)
return u ./ (u + rand(Gamma(d.beta), dims))
end
function rand!(d::Beta, A::Array{Float64})
A[:] = randbeta(d.alpha, d.beta, size(A))
for i in 1:length(A)
A[i] = rand(d)
end
return A
end
......
......@@ -16,7 +16,7 @@ immutable BetaPrime <: ContinuousUnivariateDistribution
end
end
BetaPrime() = BetaPrime(2.0, 1.0)
BetaPrime() = BetaPrime(1.0, 1.0)
cdf(d::BetaPrime, q::Real) = inc_beta(q / 1.0 + q, d.alpha, d.beta)
......@@ -24,16 +24,42 @@ insupport(d::BetaPrime, x::Number) = isreal(x) && x > 0.0 ? true : false
function mean(d::BetaPrime)
if d.beta > 1.0
d.alpha / (d.beta + 1.0)
return d.alpha / (d.beta + 1.0)
else
error("mean not defined when beta <= 1")
end
end
function modes(d::BetaPrime)
if d.alpha > 1.0
return [(d.alpha - 1.0) / (d.beta + 1.0)]
else
return [0.0]
end
end
function pdf(d::BetaPrime, x::Real)
a, b = d.alpha, d.beta
# TODO: Check the 10.0 here
return (x^(a - 1.0) * (10.0 + x)^(-(a + b))) / beta(a, b)
return (x^(a - 1.0) * (1.0 + x)^(-(a + b))) / beta(a, b)
end
rand(d::BetaPrime) = 1.0 / rand(Beta(d.alpha, d.beta))
function skewness(d::BetaPrime)
a, b = d.alpha, d.beta
if b > 3.0
return (2.0 * (2.0 * a + b - 1)) / (b - 3.0) *
sqrt((b - 2.0) / (a * (a + b - 1.0)))
else
error("skewness not defined when beta <= 3")
end
end
rand(d::BetaPrime) = 1.0 / randbeta(d.alpha, d.beta)
function var(d::BetaPrime)
a, b = d.alpha, d.beta
if b > 2.0
return (a * (a + b - 1.0)) / ((b - 2.0) * (b - 1.0)^2)
else
error("var not defined when beta <= 2")
end
end
......@@ -19,26 +19,33 @@ Binomial() = Binomial(1, 0.5)
@_jl_dist_2p Binomial binom
entropy(d::Binomial) = d.size * (-xlogx(1.0 - d.prob) - xlogx(d.prob))
insupport(d::Binomial, x::Number) = isinteger(x) && 0 <= x <= d.size
kurtosis(d::Binomial) = (1.0 - 2.0 * d.prob * (1.0 - d.prob)) / var(d)
mean(d::Binomial) = d.size * d.prob
median(d::Binomial) = iround(d.size * d.prob)
# TODO: May need to subtract 1 sometimes
modes(d::Binomial) = iround((d.size + 1.0) * d.prob)
function mgf(d::Binomial, t::Real)
n = d.size
p = d.prob
n, p = d.size, d.prob
return (1.0 - p + p * exp(t))^n
end
function cf(d::Binomial, t::Real)
n = d.size
p = d.prob
n, p = d.size, d.prob
return (1.0 - p + p * exp(im * t))^n
end
modes(d::Binomial) = iround([d.size * d.prob])
# TODO: rand() is totally screwed up
skewness(d::Binomial) = (1.0 - 2.0 * d.prob) / std(d)
var(d::Binomial) = d.size * d.prob * (1.0 - d.prob)
......@@ -40,11 +40,10 @@ function cdf(d::Categorical, x::Integer)
end
end
entropy(d::Categorical) = entropy(d.prob)
entropy(d::Categorical) = pventropy(d.prob)
function insupport(d::Categorical, x::Real)
return isinteger(x) && 1 <= x <= length(d.prob) &&
d.prob[x] != 0.0
return isinteger(x) && 1 <= x <= length(d.prob) && d.prob[x] != 0.0
end
function mean(d::Categorical)
......@@ -55,6 +54,34 @@ function mean(d::Categorical)
return s
end
function median(d::Categorical)
p, n = 0.0, length(d.prob)
i = 0
while p < 0.5 && i <= n
i += 1
p += d.prob[i]
end
return i
end
function mgf(d::Categorical, t::AbstractVector)
s = 0.0
for i in 1:length(d.prob)
s += d.prob[i] * exp(t[i])
end
return s
end
function cf(d::Categorical, t::AbstractVector)
s = 0.0 + 0.0im
for i in 1:length(d.prob)
s += d.prob[i] * exp(im * t[i])
end
return s
end
modes(d::Categorical) = [indmax(d.prob)]
pdf(d::Categorical, x::Real) = !insupport(d, x) ? 0.0 : d.prob[x]
rand(d::Categorical) = draw(d.drawtable)
......
......@@ -13,6 +13,16 @@ end
DiscreteUniform(b::Integer) = DiscreteUniform(0, b)
DiscreteUniform() = DiscreteUniform(0, 1)
function cdf(d::DiscreteUniform, k::Real)
if k < d.a
return 0.0
elseif <= d.b
return (ifloor(k) - d.a + 1.0) / (d.b - d.a + 1.0)
else
return 1.0
end
end
entropy(d::DiscreteUniform) = log(d.b - d.a + 1.0)
insupport(d::DiscreteUniform, x::Number) = isinteger(x) && d.a <= x <= d.b
......@@ -42,20 +52,15 @@ modes(d::DiscreteUniform) = [d.a:d.b]
function pdf(d::DiscreteUniform, x::Real)
if insupport(d, x)
return (1.0 / (d.b - d.a + 1.))
return (1.0 / (d.b - d.a + 1))
else
return 0.0
end
end
function quantile(d::DiscreteUniform, k::Real)
if k < d.a
return 0.0
elseif <= d.b
return (floor(k) - d.a + 1.) / (d.b - d.a + 1.)
else
return 1.0
end
function quantile(d::DiscreteUniform, p::Real)
n = d.b - d.a + 1
return d.a + ifloor(p * n)
end
rand(d::DiscreteUniform) = d.a + rand(0:(d.b - d.a))
......
......@@ -5,22 +5,38 @@
##############################################################################
immutable Erlang <: ContinuousUnivariateDistribution
shape::Float64
rate::Float64
shape::Int
scale::Float64
nested_gamma::Gamma
end
Erlang(scale::Real) = Erlang(1.0, 1.0)
Erlang() = Erlang(1.0, 1.0)
Erlang(scale::Real) = Erlang(1, 1.0)
Erlang() = Erlang(1, 1.0)
cdf(d::Erlang, x::Real) = cdf(d.nested_gamma, x)
entropy(d::Erlang) = entropy(d.nested_gamma)
insupport(d::Erlang, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x
mean(d::Erlang) = d.scale * d.shape
kurtosis(d::Erlang) = kurtosis(d.nested_gamma)
mean(d::Erlang) = d.shape * d.scale
median(d::Erlang) = median(d.nested_gamma)
mgf(d::Erlang, t::Real) = mgf(d.nested_gamma, t)
cf(d::Erlang, t::Real) = cf(d.nested_gamma, t)
modes(d::Erlang) = modes(d.nested_gamma)
function pdf(d::Erlang, x::Real)
b, c = d.scale, d.shape
return ((x / b)^(c - 1.0) * exp(-x / b)) / (b * gamma(c))
end
quantile(d::Erlang, p::Real) = quantile(d.nested_gamma, p)
function rand(d::Erlang)
b, c = d.scale, d.shape
z = 0.0
......@@ -30,4 +46,6 @@ function rand(d::Erlang)
return -b * log(z)
end
skewness(d::Erlang) = skewness(d.nested_gamma)
var(d::Erlang) = d.scale^2 * d.shape
......@@ -46,7 +46,7 @@ function invlogccdf(d::Exponential, lp::Real)
lp <= 0.0 ? -d.scale * lp : NaN
end
entropy(d::Exponential) = 1 - log(d.scale)
entropy(d::Exponential) = 1.0 - log(d.scale)
insupport(d::Exponential, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x
......@@ -66,6 +66,8 @@ function cf(d::Exponential, t::Real)
return (1.0 - t * im * s)^(-1)
end
modes(d::Exponential) = [0.0]
function pdf(d::Exponential, x::Real)
x <= 0.0 ? 0.0 : exp(-x / d.scale) / d.scale
end
......
......@@ -16,13 +16,15 @@ Gamma() = Gamma(1.0, 1.0) # Standard exponential distribution
@_jl_dist_2p Gamma gamma
function entropy(d::Gamma)
x = d.shape + log(d.scale) + lgamma(d.shape)
x = x + (1.0 - d.shape) * digamma(d.shape)
x = (1.0 - d.shape) * digamma(d.shape)
x += lgamma(d.shape) + log(d.scale) + d.shape
return x
end
insupport(d::Gamma, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x
kurtosis(d::Gamma) = 6.0 / d.shape
mean(d::Gamma) = d.shape * d.scale
function mgf(d::Gamma, t::Real)
......@@ -35,6 +37,12 @@ function cf(d::Gamma, t::Real)
return (1.0 - im * t * theta)^(-k)
end
function modes(d::Gamma)
if d.shape >= 1.0
d.scale * (d.shape - 1.0)
end
end
# rand()
#
# A simple method for generating gamma variables - Marsaglia and Tsang (2000)
......
......@@ -29,6 +29,12 @@ kurtosis(d::Geometric) = 6.0 + d.prob^2 / (1.0 - d.prob)
mean(d::Geometric) = (1.0 - d.prob) / d.prob
function median(d::Geometric)
iceil(-1.0 / log(2.0, 1.0 - d.prob)) - 1
end
modes(d::Geometric) = [0]
function mgf(d::Geometric, t::Real)
p = d.prob
if t >= -log(1.0 - p)
......
......@@ -21,8 +21,23 @@ end
cdf(d::InvertedGamma, x::Real) = 1.0 - cdf(Gamma(d.shape, d.scale), 1.0 / x)
function entropy(d::InvertedGamma)
a = d.shape
b = d.scale
return a + log(b) + lgamma(a) - (1.0 + a) * digamma(a)
end
insupport(d::InvertedGamma, x::Number) = isreal(x) && isfinite(x) && 0 <= x
function kurtosis(d::InvertedGamma)
a = d.shape
if a > 4.0
return (30.0 * a - 66.0) / ((a - 3.0) * (a - 4.0))
else
error("kurtosis not defined unless shape > 4")
end
end
function mean(d::InvertedGamma)
if d.shape > 1.0
return (1.0 / d.scale) / (d.shape - 1.0)
......@@ -31,6 +46,18 @@ function mean(d::InvertedGamma)
end
end
function mgf(d::InvertedGamma, t::Real)
a, b = d.shape, d.scale
return (2 * (-b * t)^(a / 2)) / gamma(a) *
besselk(a, sqrt(-4.0 * b * t))
end
function cf(d::InvertedGamma, t::Real)
a, b = d.shape, d.scale
return (2 * (-im * b * t)^(a / 2)) / gamma(a) *
besselk(a, sqrt(-4.0 * im * b * t))
end
modes(d::InvertedGamma) = [(1.0 / d.scale) / (d.shape + 1.0)]
pdf(d::InvertedGamma, x::Real) = exp(logpdf(d, x))
......@@ -54,6 +81,15 @@ function rand!(d::InvertedGamma, A::Array{Float64})
return A
end
function skewness(d::InvertedGamma)
a = d.shape
if a > 3.0
return (4.0 * sqrt(a - 2.0)) / (a - 3.0)
else
error("skewness not defined unless shape > 3")
end
end
function var(d::InvertedGamma)
if d.shape > 2.0
(1.0 / d.scale)^2 / ((d.shape - 1.0)^2 * (d.shape - 2.0))
......
@truncate Normal
function entropy(d::TruncatedNormal)
s = std(d.untruncated)
a = d.lower
b = d.upper
phi_a = pdf(d.untruncated, a) * s
phi_b = pdf(d.untruncated, b) * s
a_phi_a = a == -Inf ? 0.0 : a * phi_a
b_phi_b = b == Inf ? 0.0 : b * phi_b
z = d.nc
return entropy(d.untruncated) + log(z) +
0.5 * (a_phi_a - b_phi_b) / z - 0.5 * ((phi_a - phi_b) / z)^2
end
function mean(d::TruncatedNormal)
delta = pdf(d.untruncated, d.lower) - pdf(d.untruncated, d.upper)
return mean(d.untruncated) + delta * var(d.untruncated) / d.nc
......@@ -16,19 +29,14 @@ function modes(d::TruncatedNormal)
end
end
function var(d::TruncatedNormal)
s = std(d.untruncated)
a = d.lower
b = d.upper
phi_a = pdf(d.untruncated, a) * s
phi_b = pdf(d.untruncated, b) * s
a_phi_a = a == -Inf ? 0.0 : a * phi_a
b_phi_b = b == Inf ? 0.0 : b * phi_b
z = d.nc
return s^2 * (1 + (a_phi_a - b_phi_b) / z - ((phi_a - phi_b) / z)^2)
function rand(d::TruncatedNormal)
mu = mean(d.untruncated)
sigma = std(d.untruncated)
z = randnt((d.lower - mu) / sigma, (d.upper - mu) / sigma)
return mu + sigma * z
end
function entropy(d::TruncatedNormal)
function var(d::TruncatedNormal)
s = std(d.untruncated)
a = d.lower
b = d.upper
......@@ -37,8 +45,7 @@ function entropy(d::TruncatedNormal)
a_phi_a = a == -Inf ? 0.0 : a * phi_a
b_phi_b = b == Inf ? 0.0 : b * phi_b
z = d.nc
return entropy(d.untruncated) + log(z) +
0.5 * (a_phi_a - b_phi_b) / z - 0.5 * ((phi_a - phi_b) / z)^2
return s^2 * (1 + (a_phi_a - b_phi_b) / z - ((phi_a - phi_b) / z)^2)
end
# Rejection sampler based on algorithm from Robert (1992)
......@@ -88,10 +95,3 @@ function randnt(lower::Real, upper::Real)
end
end
end
function rand(d::TruncatedNormal)
mu = mean(d.untruncated)
sigma = std(d.untruncated)
z = randnt((d.lower - mu) / sigma, (d.upper - mu) / sigma)
return mu + sigma * z
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