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 77c4e75d authored by Simon Byrne's avatar Simon Byrne

Inverse Gaussian distribution

parent 306b0e4b
......@@ -51,6 +51,7 @@ export
Gumbel,
HyperGeometric,
InverseWishart,
InverseGaussian,
InvertedGamma,
Kolmogorov,
KSDist,
......@@ -211,6 +212,7 @@ include(joinpath("univariate", "erlang.jl"))
include(joinpath("univariate", "geometric.jl"))
include(joinpath("univariate", "gumbel.jl"))
include(joinpath("univariate", "hypergeometric.jl"))
include(joinpath("univariate", "inversegaussian.jl"))
include(joinpath("univariate", "invertedgamma.jl"))
include(joinpath("univariate", "kolmogorov.jl"))
include(joinpath("univariate", "ksdist.jl"))
......
# Special functions
#import Base.Math.@horner
macro horner(x, p...)
ex = p[end]
for i = length(p)-1:-1:1
ex = :($(p[i]) + $x * $ex)
end
ex
end
φ(z::Real) = exp(-0.5*z*z)/√2π
logφ(z::Real) = -0.5*(z*z + log2π)
Φ(z::Real) = 0.5*erfc(-z/√2)
Φc(z::Real) = 0.5*erfc(z/√2)
logΦ(z::Real) = z < -1.0 ? log(0.5*erfcx(-z/√2)) - 0.5*z*z : log1p(-0.5*erfc(z/√2))
logΦc(z::Real) = z > 1.0 ? log(0.5*erfcx(z/√2)) - 0.5*z*z : log1p(-0.5*erfc(-z/√2))
import Base.Math.@horner
# Rational approximations for the inverse cdf, from:
# Wichura, M.J. (1988) Algorithm AS 241: The Percentage Points of the Normal Distribution
......
immutable InverseGaussian <: ContinuousUnivariateDistribution
mu::Float64 # mean
lambda::Float64 # shape
function InverseGaussian(mu::Real, lambda::Real)
mu > zero(mu) || error("mu must be positive")
lambda > zero(lambda) || error("lambda must be positive")
new(float64(mu), float64(lambda))
end
end
InverseGaussian() = InverseGaussian(1.0, 1.0)
insupport(::InverseGaussian, x::Real) = zero(x) <= x < Inf
insupport(::Type{InverseGaussian}, x::Real) = zero(x) <= x < Inf
mean(d::InverseGaussian) = d.mu
mode(d::InverseGaussian) = (r=d.mu/d.lambda; d.mu*(sqrt(1.0+2.25*r*r)-1.5*r))
var(d::InverseGaussian) = d.mu*d.mu*d.mu/d.lambda
skewness(d::InverseGaussian) = 3.0*sqrt(d.mu/d.lambda)
kurtosis(d::InverseGaussian) = 15.0*d.mu/d.lambda
function pdf(d::InverseGaussian, x::Real)
if x <= 0.0
return 0.0
end
dd = x-d.mu
sqrt(d.lambda/(twoπ*x*x*x))*exp(-d.lambda*dd*dd/(2.0*d.mu*d.mu*x))
end
function logpdf(d::InverseGaussian, x::Real)
dd = x-d.mu
0.5*log(d.lambda/(twoπ*x*x*x))-d.lambda*dd*dd/(2.0*d.mu*d.mu*x)
end
function cdf(d::InverseGaussian, x::Real)
u = sqrt(d.lambda/x)
v = x/d.mu
Φ(u*(v-1.0)) + exp(2.0*d.lambda/d.mu)*Φ(-u*(v+1.0))
end
function quantile(d::InverseGaussian, p::Real)
if p <= 0.0 || p >= 1.0
if p == 1.0
return inf(Float64)
elseif p == 0.0
return 0.0
else
return nan(Float64)
end
end
# Whitmore and Yalovsky (1978) approximation
w = Φinv(p)
phi = d.lambda/d.mu
x = d.mu*exp(w/sqrt(phi) - 0.5/phi)
# one multiplicitive Newton iteration
xn = x*exp((p-cdf(d,x))/(pdf(d,x)*exp(x)))
# additive Newton iterations
while !isapprox(xn,x)
x = xn
xn = x + (p-cdf(d,x))/pdf(d,x)
end
xn
end
# rand method from:
# John R. Michael, William R. Schucany and Roy W. Haas (1976)
# Generating Random Variates Using Transformations with Multiple Roots
# The American Statistician , Vol. 30, No. 2, pp. 88-90
function rand(d::InverseGaussian)
z = randn()
v = z*z
w = d.mu*v
x1 = d.mu + d.mu/(2.0*d.lambda)*(w - sqrt(w*(4.0*d.lambda + w)))
p1 = d.mu / (d.mu + x1)
u = rand()
u >= p1 ? d.mu*d.mu/x1 : x1
end
......@@ -16,22 +16,12 @@ const Gaussian = Normal
zval(d::Normal, x::Real) = (x - d.μ)/d.σ
xval(d::Normal, z::Real) = d.μ + d.σ * z
φ(z::Real) = exp(-0.5*z*z)/√2π
pdf(d::Normal, x::Real) = φ(zval(d,x))/d.σ
logφ(z::Real) = -0.5*(z*z + log2π)
logpdf(d::Normal, x::Real) = logφ(zval(d,x)) - log(d.σ)
Φ(z::Real) = 0.5*erfc(-z/√2)
cdf(d::Normal, x::Real) = Φ(zval(d,x))
Φc(z::Real) = 0.5*erfc(z/√2)
ccdf(d::Normal, x::Real) = Φc(zval(d,x))
logΦ(z::Real) = z < -1.0 ? log(0.5*erfcx(-z/√2)) - 0.5*z*z : log1p(-0.5*erfc(z/√2))
logcdf(d::Normal, x::Real) = logΦ(zval(d,x))
logΦc(z::Real) = z > 1.0 ? log(0.5*erfcx(z/√2)) - 0.5*z*z : log1p(-0.5*erfc(-z/√2))
logccdf(d::Normal, x::Real) = logΦc(zval(d,x))
quantile(d::Normal, p::Real) = xval(d, Φinv(p))
......
......@@ -57,6 +57,8 @@ for d in [Arcsine(),
# HyperGeometric(3.0, 2.0, 2.0),
# HyperGeometric(2.0, 3.0, 2.0),
# HyperGeometric(2.0, 2.0, 3.0),
InverseGaussian(1.0,1.0),
InverseGaussian(2.0,7.0),
InvertedGamma(1.0,1.0),
InvertedGamma(2.0,3.0),
Laplace(0.0, 1.0),
......
......@@ -77,6 +77,8 @@ for d in [Arcsine(),
# HyperGeometric(3.0, 2.0, 2.0),
# HyperGeometric(2.0, 3.0, 2.0),
# HyperGeometric(2.0, 2.0, 3.0),
InverseGaussian(1.0,1.0),
InverseGaussian(2.0,7.0),
InvertedGamma(1.0,1.0),
InvertedGamma(2.0,3.0),
Laplace(0.0, 1.0),
......
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