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

Make NormalCanon to be a full-fledged distribution

parent c354e09b
......@@ -73,6 +73,7 @@ export
NoncentralF,
NoncentralT,
Normal,
NormalCanon,
NormalGamma,
NormalInverseGamma,
NormalInverseWishart,
......@@ -92,6 +93,7 @@ export
# methods
binaryentropy, # entropy of distribution in bits
canonform, # get canonical form of a distribution
ccdf, # complementary cdf, i.e. 1 - cdf
cdf, # cumulative distribution function
cf, # characteristic function
......@@ -250,6 +252,7 @@ include(joinpath("univariate", "noncentralchisq.jl"))
include(joinpath("univariate", "noncentralf.jl"))
include(joinpath("univariate", "noncentralt.jl"))
include(joinpath("univariate", "normal.jl"))
include(joinpath("Univariate", "normalcanon.jl"))
include(joinpath("univariate", "pareto.jl"))
include(joinpath("univariate", "poisson.jl"))
include(joinpath("univariate", "rayleigh.jl"))
......
......@@ -17,6 +17,25 @@ end
const Gaussian = Normal
# Properties
mean(d::Normal) = d.μ
median(d::Normal) = d.μ
mode(d::Normal) = d.μ
modes(d::Normal) = [d.μ]
var(d::Normal) = abs2(d.σ)
std(d::Normal) = d.σ
skewness(d::Normal) = 0.0
kurtosis(d::Normal) = 0.0
entropy(d::Normal) = 0.5 * (log2π + 1.) + log(d.σ)
mgf(d::Normal, t::Real) = exp(t * d.μ + 0.5 * d.σ^2 * t^2)
cf(d::Normal, t::Real) = exp(im * t * d.μ - 0.5 * d.σ^2 * t^2)
# Evaluation
zval(d::Normal, x::Real) = (x - d.μ)/d.σ
xval(d::Normal, z::Real) = d.μ + d.σ * z
......@@ -33,44 +52,10 @@ cquantile(d::Normal, p::Real) = xval(d, -Φinv(p))
invlogcdf(d::Normal, p::Real) = xval(d, logΦinv(p))
invlogccdf(d::Normal, p::Real) = xval(d, -logΦinv(p))
entropy(d::Normal) = 0.5 * (log2π + 1.) + log(d.σ)
kurtosis(d::Normal) = 0.0
mean(d::Normal) = d.μ
median(d::Normal) = d.μ
mgf(d::Normal, t::Real) = exp(t * d.μ + 0.5 * d.σ^2 * t^2)
cf(d::Normal, t::Real) = exp(im * t * d.μ - 0.5 * d.σ^2 * t^2)
mode(d::Normal) = d.μ
modes(d::Normal) = [d.μ]
# Sampling
rand(d::Normal) = d.μ + d.σ * randn()
skewness(d::Normal) = 0.0
std(d::Normal) = d.σ
var(d::Normal) = d.σ^2
## Canonical Form
immutable NormalCanon
h::Float64 # σ^(-2) * μ
prec::Float64 # σ^(-2)
end
Base.convert(::Type{Normal}, cf::NormalCanon) = (σ2 = 1.0 / cf.prec; Normal(σ2 * cf.h, sqrt(σ2)))
mean(cf::NormalCanon) = cf.h / cf.prec
mode(cf::NormalCanon) = cf.h / cf.prec
rand(cf::NormalCanon) = (σ2 = 1.0 / cf.prec; σ2 * cf.h + randn() * sqrt(σ2))
rand!{T<:FloatingPoint}(cf::NormalCanon, r::Array{T}) = rand!(convert(Normal, cf), r)
## Fit model
......
## Canonical Form of Normal distribution
immutable NormalCanon <: ContinuousUnivariateDistribution
h::Float64 # σ^(-2) * μ
prec::Float64 # σ^(-2)
μ::Float64 # μ
function NormalCanon(h::Float64, prec::Float64)
prec > 0. || throw(ArgumentError("prec should be positive."))
new(h, prec, h / prec)
end
NormalCanon(h::Real, prec::Real) = NormalCanon(float64(h), float64(prec))
NormalCanon() = new(0.0, 1.0, 0.0)
end
## conversion between Normal and NormalCanon
Base.convert(::Type{Normal}, cf::NormalCanon) = Normal(cf.μ, 1.0 / sqrt(cf.prec))
Base.convert(::Type{NormalCanon}, d::Normal) = (J = 1.0 / abs2(σ); NormalCanon(J * d.μ, J))
canonform(d::Normal) = convert(NormalCanon, d)
## Basic properties
@continuous_distr_support NormalCanon -Inf Inf
mean(cf::NormalCanon) = cf.μ
median(cf::NormalCanon) = mean(cf)
mode(cf::NormalCanon) = mean(cf)
modes(cf::NormalCanon) = [mean(cf)]
skewness(cf::NormalCanon) = 0.0
kurtosis(cf::NormalCanon) = 0.0
var(cf::NormalCanon) = 1.0 / cf.prec
std(cf::NormalCanon) = sqrt(var(cf))
entropy(cf::Normal) = 0.5 * (log2π + 1. - log(cf.prec))
# Evaluation
pdf(d::NormalCanon, x::Real) = (sqrt(d.prec) / √2π) * exp(-0.5 * d.prec * abs2(x - d.μ))
logpdf(d::NormalCanon, x::Real) = 0.5 * (log(d.prec) - log2π - d.prec * abs2(x - d.μ))
zval(d::NormalCanon, x::Real) = (x - d.μ) * sqrt(d.prec)
xval(d::NormalCanon, z::Real) = d.μ + z / sqrt(d.prec)
cdf(d::NormalCanon, x::Real) = Φ(zval(d,x))
ccdf(d::NormalCanon, x::Real) = Φc(zval(d,x))
logcdf(d::NormalCanon, x::Real) = logΦ(zval(d,x))
logccdf(d::NormalCanon, x::Real) = logΦc(zval(d,x))
quantile(d::NormalCanon, p::Real) = xval(d, Φinv(p))
cquantile(d::NormalCanon, p::Real) = xval(d, -Φinv(p))
invlogcdf(d::NormalCanon, p::Real) = xval(d, logΦinv(p))
invlogccdf(d::NormalCanon, p::Real) = xval(d, -logΦinv(p))
## Sampling
rand(cf::NormalCanon) = cf.μ + randn() / sqrt(cf.prec)
rand!{T<:FloatingPoint}(cf::NormalCanon, r::Array{T}) = rand!(convert(Normal, cf), r)
......@@ -84,9 +84,12 @@ for d in [Arcsine(),
NoncentralF(8,10,5),
NoncentralT(2,2),
NoncentralT(10,2),
Normal(0.0, 1.0),
Normal(),
Normal(-1.0, 10.0),
Normal(1.0, 10.0),
NormalCanon(),
NormalCanon(-1.0, 0.5),
NormalCanon(2.0, 0.8),
Pareto(),
Pareto(5.0,2.0),
Pareto(2.0,5.0),
......
......@@ -30,8 +30,7 @@ macro ignore_methoderror(ex)
end
for d in [
Arcsine(),
for d in [Arcsine(),
Bernoulli(0.1),
Bernoulli(0.5),
Bernoulli(0.9),
......@@ -108,6 +107,9 @@ for d in [
Normal(0.0, 1.0),
Normal(-1.0, 10.0),
Normal(1.0, 10.0),
NormalCanon(),
NormalCanon(-1.0, 0.5),
NormalCanon(2.0, 0.8),
Pareto(),
Pareto(5.0,2.0),
Pareto(2.0,5.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