Commit e200dc93 by Dahua Lin

### use StatsBase.make_alias_table! method

parent 837e2ee6
 julia 0.3- ArrayViews 0.4.3- PDMats 0.1.3- StatsBase 0.5- StatsBase 0.5.1-
 ... ... @@ -238,8 +238,12 @@ immutable BinomialAliasSampler <: Sampleable{Univariate,Discrete} table::AliasTable end BinomialAliasSampler(n::Int, p::Float64) = BinomialAliasSampler(make_alias_table!(binompvec(n, p))) function BinomialAliasSampler(n::Int, p::Float64) pv = binompvec(n, p) alias = Array(Int, n+1) StatsBase.make_alias_table!(pv, 1.0, pv, alias) BinomialAliasSampler(AliasTable(pv, alias, RandIntSampler(n+1))) end rand(s::BinomialAliasSampler) = rand(s.table) - 1 ... ...
 ... ... @@ -32,55 +32,13 @@ immutable AliasTable <: Sampleable{Univariate,Discrete} end ncategories(s::AliasTable) = length(s.accept) function make_alias_table!(a::AbstractVector{Float64}) # input probabilities via a, which is then # overriden by acceptance probabilities # and used internally by the AliasTable instance # n = length(a) for i=1:n @inbounds a[i] *= n end alias = Array(Int,n) larges = Int[] smalls = Int[] for i = 1:n acci = a[i] if acci > 1.0 push!(larges,i) elseif acci < 1.0 push!(smalls,i) end end while !isempty(larges) && !isempty(smalls) s = pop!(smalls) l = pop!(larges) @inbounds alias[s] = l @inbounds a[l] = (a[l] - 1.0) + a[s] if a[l] > 1 push!(larges,l) else push!(smalls,l) end end # this loop should be redundant, except for rounding for s in smalls @inbounds a[s] = 1.0 end AliasTable(a, alias, RandIntSampler(n)) end function AliasTable{T<:Real}(probs::AbstractVector{T}) n = length(probs) n > 0 || error("The input probability vector is empty.") a = Array(Float64, n) for i = 1:n @inbounds a[i] = probs[i] end make_alias_table!(a) accp = Array(Float64, n) alias = Array(Int, n) StatsBase.make_alias_table!(probs, 1.0, accp, alias) AliasTable(accp, alias, RandIntSampler(n)) end rand(s::AliasTable) = ... ...
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!