Commit 084559e6 authored by Dahua Lin's avatar Dahua Lin

samplers all moved to Sampling.jl

parent 6e99d1b8
......@@ -277,6 +277,8 @@ include(joinpath("univariate", "symtriangular.jl"))
include(joinpath("univariate", "uniform.jl"))
include(joinpath("univariate", "weibull.jl"))
include("samplers.jl")
# Multivariate distributions
include(joinpath("multivariate", "dirichlet.jl"))
include(joinpath("multivariate", "multinomial.jl"))
......@@ -296,15 +298,6 @@ include(joinpath("univariate", "truncated", "normal.jl"))
# Mixture distributions
include("mixturemodel.jl")
# Samplers
include(joinpath("samplers", "gamma.jl"))
include(joinpath("samplers", "beta.jl"))
include(joinpath("samplers", "binomial.jl"))
include(joinpath("samplers", "discreteinversetransform.jl"))
include(joinpath("samplers", "poisson.jl"))
# REPL representations
include("show.jl")
......
sampler(d::Binomial) = BinomialPolySampler(d.size, d.prob)
# delegation of samplers
rand(d::Binomial) = rand(sampler(d))
rand!(d::Binomial,a::Array) = rand!(sampler(d),a)
sampler(d::Binomial) = BinomialPolySampler(d.size, d.prob)
sampler(d::Gamma) = GammaMTSampler(d.shape, d.scale)
sampler(d::Beta) = BetaGamSampler(d::Beta)
rand(d::Beta) = rand(sampler(d))
rand!(d::Beta,a::Array) = rand!(sampler(d),a)
immutable BetaGamSampler{Ga,Gb} <: Sampler{Univariate,Continuous}
alpha::Ga
beta::Gb
end
BetaGamSampler(d::Beta) = BetaGamSampler(sampler(Gamma(d.alpha)),sampler(Gamma(d.beta)))
function rand(s::BetaGamSampler)
u = rand(s.alpha)
v = rand(s.beta)
u/(u+v)
end
# inverse transform sampler for discrete distributions
# efficient for small values, right-skewed distributions
immutable DiscreteITSampler{T<:AbstractVector} <: Sampler{Univariate,Discrete}
values::T
cdf::Vector{Float64}
end
function rand(s::DiscreteITSampler)
u = rand()
i = 0
for i = 1:length(s.cdf)
@inbounds if u <= s.cdf[i]
return s.values[i]
end
end
s.values[i+1]
end
# Sampler for drawing random number of a Gamma distribution
function sampler(d::Gamma)
d.shape > 1.0 ? GammaMTSampler(d) :
d.shape == 1.0 ? Exponential(d.scale) :
d.shape == 0.5 ? NormalSqSampler(d) :
GammaMTPowerSampler(d)
end
rand(d::Gamma) = rand(sampler(d))
rand!(d::Gamma,a::Array) = rand!(sampler(d),a)
sampler(d::Chisq) = sampler(Gamma(0.5*d.df,2.0))
rand(d::Chisq) = rand(sampler(d))
rand!(d::Chisq,a::Array) = rand!(sampler(d),a)
immutable NormalSqSampler <: Sampler{Univariate,Continuous}
scale::Float64
end
NormalSqSampler(d::Gamma) = NormalSqSampler(0.5*d.scale)
rand(s::NormalSqSampler) = s.scale*(randn())^2
# A simple method for generating gamma variables - Marsaglia and Tsang (2000)
# http://www.cparity.com/projects/AcmClassification/samples/358414.pdf
# Page 369
# basic simulation loop for pre-computed d and c
immutable GammaMTSampler <: Sampler{Univariate,Continuous}
d::Float64
c::Float64
scale::Float64
end
function GammaMTSampler(g::Gamma)
d = g.shape - 1.0/3.0
GammaMTSampler(d, 1.0 / sqrt(9.0 * d), g.scale * d)
end
# for shape < 1.0: see note on page 371.
immutable GammaMTPowerSampler <: Sampler{Univariate,Continuous}
gmt::GammaMTSampler
::Float64
end
function GammaMTPowerSampler(g::Gamma)
d = g.shape + 2.0/3.0
GammaMTPowerSampler(GammaMTSampler(d, 1.0 / sqrt(9.0 * d), g.scale * d), 1.0/g.shape)
end
function rand(s::GammaMTSampler)
v = 0.0
while true
x = randn()
v = 1.0 + s.c * x
while v <= 0.0
x = randn()
v = 1.0 + s.c * x
end
v *= (v * v)
u = rand()
x2 = x^2
if u < 1.0 - 0.331 * x2^2
break
end
if log(u) < 0.5 * x2 + s.d * (1.0 - v + log(v))
break
end
end
v * s.scale
end
rand(s::GammaMTPowerSampler) = rand(s.gmt) * rand()^s.
sampler(d::Poisson) = d.lambda <= 10.0 ? DiscreteITSampler(d) : PoissonADSampler(d)
rand(d::Poisson) = rand(sampler(d))
rand!(d::Poisson,a::Array) = rand!(sampler(d),a)
function DiscreteITSampler(d::Poisson)
d.lambda <= 10.0 || error("lambda too large")
values = 0:44
cdf = Array(Float64,44)
p = exp(-d.lambda)
c = p
cdf[1] = p
for i = 1:43
p *= d.lambda/i
c += p
cdf[i+1] = c
end
DiscreteITSampler(values,cdf)
end
# algorithm from:
# J.H. Ahrens, U. Dieter (1982)
# "Computer Generation of Poisson Deviates from Modified Normal Distributions"
# ACM Transactions on Mathematical Software, 8(2):163-179
# μ >= 10.0
immutable PoissonADSampler <: Sampler{Univariate,Discrete}
μ::Float64
s::Float64
d::Float64
L::Int
end
function PoissonADSampler(d::Poisson)
μ = d.lambda
PoissonADSampler(μ,sqrt(μ),6.0*μ^2,ifloor(μ-1.1484))
end
function rand(s::PoissonADSampler)
# Step N
G = s.μ + s.s*randn()
if G >= 0.0
K = ifloor(G)
# Step I
if K >= s.L
return K
end
# Step S
U = rand()
if s.d*U >= (s.μ-K)^3
return K
end
# Step P
px,py,fx,fy = procf(s.μ,K,s.s)
# Step Q
if fy*(1-U) <= py*exp(px-fx)
return K
end
end
while true
# Step E
E = Base.Random.randmtzig_exprnd()
U = 2.0*rand()-1.0
T = 1.8+copysign(E,U)
if T <= -0.6744
continue
end
K = ifloor(s.μ + s.s*T)
px,py,fx,fy = procf(s.μ,K,s.s)
c = 0.1069/s.μ
# Step H
if c*abs(U) <= py*exp(px+E)-fy*exp(fx+E)
return K
end
end
end
# Procedure F
function procf(μ,K,s)
# can be pre-computed, but does not seem to affect performance
ω = 0.3989422804014327/s
b1 = 0.041666666666666664/μ
b2 = 0.3*b1*b1
c3 = 0.14285714285714285*b1*b2
c2 = b2 - 15.0*c3
c1 = b1 - 6.0*b2 + 45.0*c3
c0 = 1.0 - b1 + 3.0*b2 - 15.0*c3
if K < 10
px = -μ
py = μ^K/factorial(K)
else
δ = 0.08333333333333333/K
δ -= 4.8*δ^3
V = (μ-K)/K
px = K*log1pmx(V) - δ # avoids need for table
py = 0.3989422804014327/sqrt(K)
end
X = (K-μ+0.5)/s
X2 = X^2
fx = -0.5*X2 # missing negation in pseudo-algorithm, but appears in fortran code.
fy = ω*(((c3*X2+c2)*X2+c1)*X2+c0)
return px,py,fx,fy
end
......@@ -24,63 +24,6 @@ logexpm1(x::Integer) = logexpm1(float(x))
# log(1+x^2)
log1psq(x::FloatingPoint) = (ax = abs(x); ax < maxintfloat(x) ? log1p(ax*ax) : 2*log(ax))
# log(1+x)-x
# accurate ~2ulps for -0.227 < x < 0.315
function log1pmx_kernel(x::Float64)
r = x/(x+2.0)
t = r*r
w = @horner(t,
6.66666666666666667e-1, # 2/3
4.00000000000000000e-1, # 2/5
2.85714285714285714e-1, # 2/7
2.22222222222222222e-1, # 2/9
1.81818181818181818e-1, # 2/11
1.53846153846153846e-1, # 2/13
1.33333333333333333e-1, # 2/15
1.17647058823529412e-1) # 2/17
hxsq = 0.5*x*x
r*(hxsq+w*t)-hxsq
end
# use naive calculation or range reduction outside kernel range.
# accurate ~2ulps for all x
function log1pmx(x::Float64)
if !(-0.7 < x < 0.9)
return log1p(x) - x
elseif x > 0.315
u = (x-0.5)/1.5
return log1pmx_kernel(u) - 9.45348918918356180e-2 - 0.5*u
elseif x > -0.227
return log1pmx_kernel(x)
elseif x > -0.4
u = (x+0.25)/0.75
return log1pmx_kernel(u) - 3.76820724517809274e-2 + 0.25*u
elseif x > -0.6
u = (x+0.5)*2.0
return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5*u
else
u = (x+0.625)/0.375
return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625*u
end
end
# log(x) - x + 1
function logmxp1(x::Float64)
if x <= 0.3
return (log(x) + 1.0) - x
elseif x <= 0.4
u = (x-0.375)/0.375
return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625*u
elseif x <= 0.6
u = 2.0*(x-0.5)
return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5*u
else
return log1pmx(x-1.0)
end
end
φ(z::Real) = exp(-0.5*z*z)/√2π
logφ(z::Real) = -0.5*(z*z + log2π)
......
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