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

Merge commit 'bc94efaa'

parents 7518365e bc94efaa
......@@ -218,6 +218,7 @@ List of Distributions
- :ref:`pareto`
- :ref:`rayleigh`
- :ref:`tdist`
- :ref:`triangular`
- :ref:`uniform`
- :ref:`vonmises`
- :ref:`weibull`
......@@ -732,6 +733,25 @@ The probability density function of a `Students T distribution <http://en.wikipe
TDist(d) # t-distribution with d degrees of freedom
.. _triangular:
Triangular Distribution
~~~~~~~~~~~~~~~~~~~~~~~~~
The probability density function of a `Triangular distribution <http://en.wikipedia.org/wiki/Triangular_distribution>`_ with lower limit a, upper limit b and mode c is
.. math::
f(x; a, b, c)= \begin{cases}
0 & \mathrm{for\ } x < a, \\
\frac{2(x-a)}{(b-a)(c-a)} & \mathrm{for\ } a \le x \leq c, \\[4pt]
\frac{2(b-x)}{(b-a)(b-c)} & \mathrm{for\ } c < x \le b, \\[4pt]
0 & \mathrm{for\ } b < x,
\end{cases}
.. code-block:: julia
TriangularDist(a, b, c) # Triangular distribution with lower limit a, upper limit b and mode c
.. _uniform:
......
......@@ -111,6 +111,7 @@ export
Skellam,
SymTriangularDist,
TDist,
TriangularDist,
Truncated,
TruncatedNormal,
Uniform,
......
# Triangular distribution
immutable TriangularDist <: ContinuousUnivariateDistribution
a::Float64
b::Float64
c::Float64
function TriangularDist(a::Real, b::Real, c::Real)
a < b || error("a<b must be true")
a <= c <= b || error("a<=c<=b must be true")
new(float64(a), float64(b), float64(c))
end
end
## Support
isupperbounded(::Union(TriangularDist, Type{TriangularDist})) = true
islowerbounded(::Union(TriangularDist, Type{TriangularDist})) = true
isbounded(::Union(TriangularDist, Type{TriangularDist})) = true
minimum(d::TriangularDist) = d.a
maximum(d::TriangularDist) = d.b
insupport(d::TriangularDist, x::Real) = minimum(d) <= x <= maximum(d)
## Properties
mean(d::TriangularDist) = (d.a + d.b + d.c) / 3.0
median(d::TriangularDist) = d.c >= (d.a+d.b)/2.0 ?
d.a + sqrt((d.b-d.a)*(d.c-d.a))/sqrt(2.0) :
d.b - sqrt((d.b-d.a)*(d.b-d.c))/sqrt(2.0)
mode(d::TriangularDist) = d.c
var(d::TriangularDist) = (d.a^2.0 + d.b^2.0 + d.c^2.0-d.a*d.b-d.a*d.c-d.b*d.c)/18.0
skewness(d::TriangularDist) = sqrt(2.0)*(d.a+d.b-2.0d.c)*(2d.a-d.b-d.c)*(d.a-2d.b+d.c)/5.0/(d.a^2.0 + d.b^2.0 + d.c^2.0-d.a*d.b-d.a*d.c-d.b*d.c)^(3.0/2.0)
kurtosis(d::TriangularDist) = -0.6
entropy(d::TriangularDist) = 0.5 + log((d.b-d.a)/2.0)
## Functions
function pdf(d::TriangularDist, x::Real)
if d.a<=x<=d.c
return 2.0*(x-d.a)/(d.b-d.a)/(d.c-d.a)
elseif d.c<x<=d.b
return 2.0*(d.b-x)/(d.b-d.a)/(d.b-d.c)
else
return zero(x)
end
end
function cdf(d::TriangularDist, x::Real)
if x<d.a
return zero(x)
elseif d.a<=x<=d.c
return (x-d.a)^2/(d.b-d.a)/(d.c-d.a)
elseif d.c<x<=d.b
return 1.0-(d.b-x)^2/(d.b-d.a)/(d.b-d.c)
else
return one(x)
end
end
function quantile(d::TriangularDist, p::Real)
@checkquantile p begin
if p <= (d.c-d.a)/(d.b-d.a)
return d.a + sqrt((d.b-d.a)*(d.c-d.a)*p)
else
return d.b - sqrt((d.b-d.a)*(d.b-d.c)*(1-p))
end
end
end
function mgf(d::TriangularDist, t::Real)
if t==zero(t)
return one(t)
else
nominator = (d.b-d.c)*exp(d.a*t)-(d.b-d.a)*exp(d.c*t)+(d.c-d.a)*exp(d.b*t)
denominator = (d.b-d.a)*(d.c-d.a)*(d.b-d.c)*t^2
return 2*nominator/denominator
end
end
function cf(d::TriangularDist, t::Real)
# Is this correct?
if t==zero(t)
return one(t)
else
nominator = (d.b-d.c)*exp(im*d.a*t)-(d.b-d.a)*exp(im*d.c*t)+(d.c-d.a)*exp(im*d.b*t)
denominator = (d.b-d.a)*(d.c-d.a)*(d.b-d.c)*t^2
return -2*nominator/denominator
end
end
## Sampling
function rand(d::TriangularDist)
u = rand()
if u < (d.c-d.a)/(d.b-d.a)
return d.a + sqrt(u*(d.b-d.a)*(d.c-d.a))
else
return d.b - sqrt((1.0-u)*(d.b-d.a)*(d.b-d.c))
end
end
......@@ -152,6 +152,7 @@ for finame in ["arcsine.jl",
"rayleigh.jl",
"skellam.jl",
"tdist.jl",
"triangular.jl",
"symtriangular.jl",
"uniform.jl",
"vonmises.jl",
......
......@@ -50,6 +50,9 @@ Rayleigh(8.0)
TDist(1.2)
TDist(5.0)
TDist(28.0)
TriangularDist(-4, 14, 3),
TriangularDist(2, 2000, 500),
TriangularDist(1, 3, 2),
TruncatedNormal(0, 1, -2, 2)
TruncatedNormal(3, 10, 7, 8)
TruncatedNormal(27, 3, 0, Inf)
......
......@@ -19,7 +19,8 @@ tests = [
"conjugates_mvnormal",
"wishart",
"mixture",
"gradlogpdf"]
"gradlogpdf",
"triangular"]
print_with_color(:blue, "Running tests:\n")
......
using Distributions
using Base.Test
# Test constructors
@test_throws MethodError TriangularDist()
@test_throws MethodError TriangularDist(3.0)
@test_throws MethodError TriangularDist(4.5, 2.4)
@test_throws ErrorException TriangularDist(4.0, 2.0, 3.0)
@test_throws ErrorException TriangularDist(3.0, 3.0, 3.0)
@test_throws ErrorException TriangularDist(3.0, 5.0, 2.0)
@test_throws ErrorException TriangularDist(3.0, 5.0, 6.0)
# Create two distributions for further tests
# There are two so that all code-paths of some functions
# can be tested.
a1 = 3.0
b1 = 5.0
c1 = 7.0/2.0
d1 = TriangularDist(a1,b1,c1)
a2 = 3.0
b2 = 5.0
c2 = 9.0/2.0
d2 = TriangularDist(a2,b2,c2)
@test isupperbounded(d1)==true
@test isupperbounded(TriangularDist)==true
@test islowerbounded(d1)==true
@test islowerbounded(TriangularDist)==true
@test isbounded(d1)==true
@test isbounded(TriangularDist)==true
@test minimum(d1)==a1
@test maximum(d1)==b1
@test insupport(d1, 3.0)==true
@test insupport(d1, 5.0)==true
@test insupport(d1, 3.5)==true
@test insupport(d1, 3.1)==true
@test insupport(d1, 4.999)==true
@test insupport(d1, 1.0)==false
@test insupport(d1, 5.5)==false
@test_approx_eq mean(d1) 23./6.
@test_approx_eq median(d1) 5.-sqrt(3./2.0)
@test_approx_eq median(d2) 3.+sqrt(3./2.0)
@test_approx_eq mode(d1) 7./2.
@test_approx_eq var(d1) 13/72
@test_approx_eq skewness(d1) 14*sqrt(2/13)/13
@test_approx_eq kurtosis(d1) 12/5-3
# Double check whether this test is correct
@test_approx_eq entropy(d1) 0.5
@test_approx_eq pdf(d1, 1) 0
@test_approx_eq pdf(d1, 3.1) 0.2
@test_approx_eq pdf(d1, 3.5) 1
@test_approx_eq pdf(d1, 4.0) 2/3
@test_approx_eq pdf(d1, 5.0) 0
@test_approx_eq pdf(d1, 5.1) 0
# Not sure these make sense
@test_approx_eq mgf(d1, 2) 3141.7032569285
# Not sure these make sense
@test_approx_eq cf(d1, 1) -0.706340963216717-0.578379106849962im
@test_approx_eq cf(d1, 3) 0.117846739742880-0.400401950797709im
@test_approx_eq cf(d1, 20) 0.00754699480217164+0.00752726705107994im
@test isnan(quantile(d1,-1))
@test isnan(quantile(d1,6))
@test_approx_eq quantile(d1, 0) a1
@test_approx_eq quantile(d1, 1) b1
@test_approx_eq quantile(d1, 0.1) 3+1/sqrt(10)
@test_approx_eq quantile(d1, 0.4) 5-3/sqrt(5)
@test_approx_eq quantile(d1, 0.5) 5-sqrt(3/2)
@test_approx_eq quantile(d1, 0.9) 5-sqrt(3/10)
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