CategoricalMonteCarlo.algorithm2_1!Method
algorithm2_1!(p::Vector{T}, I::Vector{Int}, w::Vector{<:Real}) where {T<:Real}

Fill p with the probabilities that result from normalizing the weights w[I]. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_1

CategoricalMonteCarlo.algorithm2_1Method
algorithm2_1(I::Vector{Int}, w::Vector{<:Real})

Create a vector of probabilities by normalizing the weights selected by I from w. It is assumed that 0 ≤ wᵢ < Inf, i ∈ I.

Mathematically, given:

I ∈ ℕᴺ, 𝐰 ∈ ℝᴰ; I ⊆ {1,…,D}

The iᵗʰ term will be computed as: pᵢ = 𝐰ᵢ / ∑ⱼ 𝐰ⱼ; j ∈ I

Examples

julia> I = [1, 5, 2]; w = [5, 4, 3, 2, 1];

julia> algorithm2_1(I, w)
3-element Vector{Float64}:
0.5
0.1
0.4

julia> algorithm2_1([1, 1, 2], Rational.(w))
3-element Vector{Rational{Int64}}:
5//14
5//14
2//7

julia> w[2] = -w[2];

julia> algorithm2_1(I, w)                # Nonsense results if wᵢ constraints violated
3-element Vector{Float64}:
2.5
0.5
-2.0

julia> algorithm2_1(I, [5, 4, 3, 2, Inf])
3-element Vector{Float64}:
0.0
NaN
0.0

julia> algorithm2_1(I, [5, NaN, 3, 2, 1])
3-element Vector{Float64}:
NaN
NaN
NaN
CategoricalMonteCarlo.algorithm2_1_algorithm3!Method
algorithm2_1_algorithm3!(p::Vector{T}, I::Vector{Int}, w::Vector{<:Real}, u::Real) where {T<:Real}

Normalize w[I] to probabilities, storing the result in p, then spreading probability mass 0 ≤ u ≤ 1 across the 0 or more elements of w[I] which are equal to zero. If all values of w[I] are zero and u ≠ 0, p will be filled with uniform probability mass. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_1_algorithm3

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; u = 0.5;

julia> algorithm2_1_algorithm3!(similar(I, Rational{Int}), I, w, u)
4-element Vector{Rational{Int64}}:
1//6
1//4
1//4
1//3
CategoricalMonteCarlo.algorithm2_1_algorithm3Method
algorithm2_1_algorithm3(I::Vector{Int}, w::Vector{<:Real}, u::Real)

Return a vector of probabilities, normalizing the components selected from w by the index set I, then spreading the probability mass 0 ≤ u ≤ 1 across the 0 or more elements which are equal to zero. If all values of w[I] are zero and u ≠ 0, a vector of uniform probability mass is returned. Equivalent to algorithm3!(algorithm2_1(I, w), u) but more efficient.

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; u = 0.5;

julia> algorithm2_1_algorithm3(I, w, u)
4-element Vector{Float64}:
0.16666666666666666
0.25
0.25
0.3333333333333333

julia> algorithm3!(algorithm2_1(I, w), u)
4-element Vector{Float64}:
0.16666666666666666
0.25
0.25
0.3333333333333333
CategoricalMonteCarlo.algorithm2_1_algorithm3_ratio!Method
algorithm2_1_algorithm3!(p::Vector{T}, I::Vector{Int}, w::Vector{<:Real}, u::Real) where {T<:Real}

Normalize w[I] to probabilities, storing the result in p, then spreading probability mass u = r / (1 + r) across the 0 or more elements of w[I] which are equal to zero such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. Note that T must be a type which is able to hold the result of inv(one(T)).

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; r = 1.0;

julia> algorithm2_1_algorithm3_ratio!(similar(I, Rational{Int}), I, w, r)
4-element Vector{Rational{Int64}}:
1//6
1//4
1//4
1//3
CategoricalMonteCarlo.algorithm2_1_algorithm3_ratioMethod
algorithm2_1_algorithm3_ratio(I::Vector{Int}, w::Vector{<:Real}, r::Real)

Return a vector of probabilities, normalizing the components selected from w by the index set I, then spreading the probability mass u = r / (1 + r) across the 0 or more elements of w[I] which are equal to zero such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. If all values of w[I] are zero and r ≠ 0, a vector of uniform probability mass is returned. Equivalent to algorithm3_ratio!(algorithm2_1(I, w), u) but more efficient.

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; r = 1.0;

julia> algorithm2_1_algorithm3_ratio(I, w, r)
4-element Vector{Float64}:
0.16666666666666666
0.25
0.25
0.3333333333333333

julia> algorithm3_ratio!(algorithm2_1(I, w), r)
4-element Vector{Float64}:
0.16666666666666666
0.25
0.25
0.3333333333333333
CategoricalMonteCarlo.algorithm2_2!Method
algorithm2_2!(p::Vector{T}, Is::NTuple{M, Vector{Int}}, ws::NTuple{M, Vector{<:Real}}) where {T<:Real, M}

Fill p with the probabilities that result from normalizing the element-wise product of weights selected by the index set, Is[m], respective to each weight vector, ws[m]. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_2

Examples

julia> Is = ([1,2,3], [4,5,6], [7,8,9]); ws = ([1,2,3], fill(1/6, 6), fill(1//10, 9));

julia> algorithm2_2!(zeros(3), Is, ws)
3-element Vector{Float64}:
0.16666666666666666
0.3333333333333333
0.5

julia> algorithm2_2!(zeros(Rational{Int}, 3), Is, (ws[1], fill(1//6, 6), ws[3]))
3-element Vector{Rational{Int64}}:
1//6
1//3
1//2
CategoricalMonteCarlo.algorithm2_2Method
algorithm2_2(Is::NTuple{M, Vector{Int}}, ws::NTuple{M, Vector{<:Real}}) where {M}

Compute the product of weights selected by the respective index sets Is, then normalize the resultant weight vector to probabilities.

Mathematically, given:

I₁ ∈ ℕᴺ , 𝐰₁ ∈ ℝᴰ¹, 0 ≤ 𝐰₁ᵢ < Inf, i ∈ I₁

I₂ ∈ ℕᴺ , 𝐰₂ ∈ ℝᴰ², 0 ≤ 𝐰₂ᵢ < Inf, i ∈ I₂

⋮ , ⋮

Iₘ ∈ ℕᴺ , 𝐰ₘ ∈ ℝᴰᵐ, 0 ≤ 𝐰ₘᵢ < Inf, i ∈ Iₘ

The iᵗʰ term will be computed as: pᵢ = ∏ₘ₌₁ᴹ 𝐰ₘ[Iₘ[i]] / ∑ⱼ₌₁ᴺ ∏ₘ₌₁ᴹ 𝐰ₘ[Iₘ[j]]

Examples

julia> Is = ([1,2,3], [4,5,6], [7,8,9]); ws = ([1.0, 2.0, 3.0], fill(0.5, 6), fill(0.1, 9));

julia> algorithm2_2(Is, ws)
3-element Vector{Float64}:
0.16666666666666666
0.3333333333333333
0.5

julia> w = ws[1][Is[1]] .* ws[2][Is[2]] .* ws[3][Is[3]]    # unnormalized
3-element Vector{Float64}:
0.05
0.1
0.15000000000000002

julia> w ./= sum(w)                                        # normalized
3-element Vector{Float64}:
0.16666666666666666
0.3333333333333333
0.5
CategoricalMonteCarlo.algorithm3!Method
algorithm3!(p::Vector{T}, w::Vector{<:Real}, u::Real) where {T<:Real}

Normalize w to probabilities, storing the result in p, spreading probability mass 0 ≤ u ≤ 1 across the 0 or more elements of w which are equal to zero. If all values of w are zero and u ≠ 0, p will be filled with uniform probability mass. Note that T must be a type which is able to hold the result of inv(one(T)).

Examples

julia> w = [0, 10, 5, 1]; u = 0.5;

julia> algorithm3!(similar(w, Float64), w, u)
4-element Vector{Float64}:
0.5
0.3125
0.15625
0.03125
CategoricalMonteCarlo.algorithm3!Method
algorithm3!(p::Vector{T}, u::Real) where {T<:Real}

Normalize p to probabilities, spreading probability mass u across the 0 or more elements of p which are equal to zero. If all values of w are zero and u ≠ 0, p will be filled with uniform probability mass. Note that T must be a type which is able to hold the result of inv(one(T)).

Examples

julia> algorithm3!(Rational{Int}[0, 10, 5, 0], 0.5)
4-element Vector{Rational{Int64}}:
1//4
1//3
1//6
1//4
CategoricalMonteCarlo.algorithm3Method
algorithm3(w::Vector{<:Real}, u::Real)

Return the vector of probabilities created by normalizing w to probabilities, then spreading the probability mass 0 ≤ u ≤ 1 across the 0 or more elements of w which are equal to zero. If all values of w are zero and u ≠ 0, a vector of uniform probability mass is returned.

Mathematically, given:

𝐰 ∈ ℝᴺ, 0 ≤ wᵢ < ∞, u ∈ ℝ, 0 ≤ u ≤ 1, J = {i : 𝐰ᵢ = 0}

pᵢ =
Case 1: if J ≠ ∅
u / |J|                     if i ∈ J
(1 - u) * 𝐰ᵢ / ∑ᵢ₌₁ᴺ 𝐰ᵢ     otherwise
Case 2: if J = {1,…,N}
u / (u * N)                 Equivalent to 1/N if 𝐰 ≠ ̲0

Examples

julia> algorithm3([0, 10, 5, 0], 0.5)
4-element Vector{Float64}:
0.25
0.3333333333333333
0.16666666666666666
0.25

julia> algorithm3(Rational{Int}[0, 0, 0], 0.25)    # 𝐰 = ̲0
3-element Vector{Rational{Int64}}:
1//3
1//3
1//3

julia> algorithm3([0, 0], 0.0)                     # 𝐰 = ̲0 and u = 0
2-element Vector{Float64}:
NaN
NaN

julia> algorithm3([1, 2, 3], 0.9)                  # in absence of 0's, just normalize
3-element Vector{Float64}:
0.16666666666666666
0.3333333333333333
0.5
CategoricalMonteCarlo.algorithm3_ratio!Method
algorithm3_ratio!(p::Vector{T}, w::Vector{<:Real}, r::Real) where {T<:Real}

Normalize w to probabilities, storing the result in p, then spread the probability mass u = r / (1 + r) across the 0 or more elements of w such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. Note that T must be a type which is able to hold the result of inv(one(T)).

Examples

julia> w = [0, 10, 5, 1]; r = 1.0;

julia> algorithm3_ratio!(similar(w, Float64), w, r)
4-element Vector{Float64}:
0.5
0.3125
0.15625
0.03125
CategoricalMonteCarlo.algorithm3_ratio!Method
algorithm3_ratio!(p::Vector{T}, r::Real) where {T<:Real}

Normalize p to probabilities, then spread the probability mass u = r / (1 + r) across the 0 or more elements of p such that the ratio of the sum of (inititally) zero elements to the sum of the non-zero elements is equal to r. Note that T must be a type which is able to hold the result of inv(one(T)).

CategoricalMonteCarlo.algorithm3_ratioMethod
algorithm3_ratio(w::Vector{<:Real}, r::Real)

Return a vector of probabilities created by normalizing w to probabilities, then spread the probability mass u = r / (1 + r) across the 0 or more elements of w which are equal to zero such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. If all values of w are zero and r ≠ 0, a vector of uniform probability mass is returned.

Mathematically, given:

𝐰 ∈ ℝᴺ, 0 ≤ wᵢ < ∞, r ∈ ℝ, 0 ≤ r ≤ Inf, J = {i : 𝐰ᵢ = 0}

pᵢ =
Case 1: if J ≠ ∅
(r / (1+r)) / |J|                     if i ∈ J
(1 / (1+r)) * 𝐰ᵢ / ∑ᵢ₌₁ᴺ 𝐰ᵢ           otherwise
Case 2: if J = {1,…,N}
r / (r * N)                           Equivalent to 1/N if 𝐰 ≠ ̲0

Examples

julia> w = Rational{Int}[1, 0, 3, 0, 5]; r = 3;

julia> p = algorithm3_ratio(w, r)
5-element Vector{Rational{Int64}}:
1//36
3//8
1//12
3//8
5//36

julia> r′ = sum(p[findall(iszero, w)]) / sum(p[findall(!iszero, w)]); (r′, r′ == r)
(3//1, true)

julia> algorithm3(w, r / (1 + r))              # Note equivalence
5-element Vector{Rational{Int64}}:
1//36
3//8
1//12
3//8
5//36

julia> algorithm3_ratio(w, Inf)                # r = Inf ⟹ u = 1
5-element Vector{Rational{Int64}}:
0//1
1//2
0//1
1//2
0//1
CategoricalMonteCarlo.algorithm4!Method
algorithm4!(p::Vector{T}, 𝐰₁::Vector{<:Real}, 𝐰₂::Vector{<:Real}) where {T<:Real}

Fill p with the probabilities which result from algorithm4(𝐰₁, 𝐰₂). Note that T must be a type which is able to hold the result of inv(one(T)).

CategoricalMonteCarlo.algorithm4!Method
algorithm4!(𝐰₁::Vector{T}, 𝐰₂::Vector{<:Real}) where {T<:Real}

Fill 𝐰₁ with the probabilities which result from algorithm4(𝐰₁, 𝐰₂). Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm4

CategoricalMonteCarlo.algorithm4Method
algorithm4(𝐰₁::Vector{<:Real}, 𝐰₂::Vector{<:Real})

Return a vector of probabilities constructed according to the following algorithm:

Define:

I = {1,…,N}

J₁ = {i: 𝐰₁ᵢ = 0}, I₁′ = {i: 𝐰₁ᵢ ≠ 0} = I ∖ J₁

J₂ = {i: 𝐰₂ᵢ = 0}, I₂′ = {i: 𝐰₂ᵢ ≠ 0} = I ∖ J₂

𝐰₁ ∈ ℝᴺ : initial weights, 0 ≤ 𝐰₁ᵢ < Inf

𝐰₂ ∈ ℝᴺ : augment weights, 0 ≤ 𝐰₂ᵢ < Inf; a value of zero indicates no re-weight

Then:

pᵢ = 𝐰₁ᵢ / ∑ₗ₌₁ᴺ 𝐰₁ₗ, i ∈ I ∖ I₂′

pᵢ = (𝐰₂ᵢ * ∑ₗ 𝐰₁ₗ, l ∈ I₂′) / (∑ₗ₌₁ᴺ 𝐰₂ₗ * ∑ₗ₌₁ᴺ 𝐰₁ₗ), i ∈ I₂′

This algorithm can produce a wide variety of probability vectors as the result of the various combinations of intersections which can be formed from J₁, J₂, I₁′, and I₂′. However, complexity of outputs aside, the motivating concept is quite simple: take a vector of weights, 𝐰₁ and re-weight some subset (I₂′) of those weights using a second set of weights, 𝐰₂, while preserving the proportion of probability mass derived from 𝐰₁. That is, given p = algorithm4(𝐰₁, 𝐰₂), the following relationships are preserved: sum(p[J₂]) ≈ sum(𝐰₁[J₂]) / sum(𝐰₁[I₁′]), sum(w₁[J₂]) / sum(w₁[I₂′]) ≈ sum(p[J₂]) / sum(p[I₂′]).

See also: algorithm4!

Examples

julia> w₁ = [1, 1, 1, 1, 0];

julia> algorithm4(w₁, [2, 1, 3, 4, 0])    # J₁ ∩ I₂′ = ∅
5-element Vector{Float64}:
0.2
0.1
0.30000000000000004
0.4
0.0

julia> algorithm4(w₁, [2, 1, 3, 0, 5])    # J₂ = [4] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
0.13636363636363635
0.06818181818181818
0.20454545454545453
0.25
0.3409090909090909

julia> w₁ = [1, 1, 1, 0, 0];

julia> algorithm4(w₁, [2, 1, 3, 4, 0])    # J₂ = [5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
0.2
0.1
0.30000000000000004
0.4
0.0

julia> w₁ = [1, 1, 0, 1, 0];

julia> algorithm4(w₁, [0, 1, 0, 4, 0])    # J₂ = [1,3,5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
0.3333333333333333
0.13333333333333333
0.0
0.5333333333333333
0.0

julia> algorithm4(w₁, [0, 0, 3, 4, 0])    # J₂ = [1,2,5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
0.3333333333333333
0.3333333333333333
0.14285714285714285
0.19047619047619047
0.0

julia> algorithm4(w₁, [2, 0, 3, 0, 0])    # J₂ = [2,4,5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
0.13333333333333333
0.3333333333333333
0.2
0.3333333333333333
0.0
CategoricalMonteCarlo.normalize1!Method
normalize1!(A::AbstractArray{<:Real})

Normalize the values in A such that sum(A) ≈ 1 and 0 ≤ A[i] ≤ 1 ∀i. This is not quite the L¹-norm, which would require that abs(A[i]) be used. It is assumed that 0 ≤ A[i] < Inf ∀i. Inf values are not handled and will result in NaN's.

See also: normalize1

julia> normalize1!([1.0, 2.0, 3.0])
3-element Vector{Float64}:
0.16666666666666666
0.3333333333333333
0.5

julia> normalize1!([1.0, 2.0, Inf])
3-element Vector{Float64}:
0.0
0.0
NaN

julia> normalize1!([1.0, 2.0, NaN])     # NaN propagates, as expected
3-element Vector{Float64}:
NaN
NaN
NaN

julia> normalize1!([1.0, -2.0, 3.0])    # not the L¹-norm
3-element Vector{Float64}:
0.5
-1.0
1.5
CategoricalMonteCarlo.normalize1!Method
normalize1!(B::AbstractArray{<:Real}, A::AbstractArray{<:Real})

Normalize the values in A such that sum(B) ≈ 1 and 0 ≤ B[i] ≤ 1 ∀i, storing the result in B. It is assumed that A[i] ≥ 0 ∀i.

CategoricalMonteCarlo.sample!Method
sample!(B::AbstractArray, A::AbstractArray)

Draw samples, summing (and potentially reducing) in-place into B. The shape of B determines the extent of reduction performed.

CategoricalMonteCarlo.sampleMethod
sample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing])

Draw n_sim samples from the distribution which corresponds to the sum of the independent categorical distributions defined by the probability mass vector(s), A, storing the result in an array of eltype T of sufficient size and dimension as determined by the input, A, the number of categories, n_cat, and the (potential) reduction dimensions, dims.

dims is an optional keyword argument used to specify an in-place sum on the indices of A (if A is an array of arrays). n_cat may be optionally specified, or inferred from the available data; validity indices will be checked regardless.

In the simplest case, A may be a single probability mass vector, for which compatible types are AbstractVector{<:Real}, SparseVector{<:Real, <:Integer}, and AbstractVector{Tuple{<:Integer, <:Real}}, the last of which is simply another representation of a sparse vector.

In the second case, A may be an AbstractArray{V, N} where {V, N} in which V is any of the types specified above for a single probability mass vector. That is, an array of vectors.

In the third case, A may be an AbstractArray{W, M} where {W<:AbstractArray{V, N}, M} where {V, N} in which W is any of the types specified in the second case. In other words, an array of arrays of vectors.

CategoricalMonteCarlo.splitrangesMethod
splitranges(start, stop, chunksize)

Divide the range start:stop into segments, each of size chunksize. The last segment will contain the remainder, (start - stop + 1) % chunksize, if it exists.

CategoricalMonteCarlo.tsample!Method
tsample!(B::AbstractArray, A::AbstractArray; [chunksize=5000])

Identical to sample! except that thread-based parallelism is used to accelerate the computation.

CategoricalMonteCarlo.tsampleMethod
tsample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing], [chunksize=5000])

See sample for full documentation; identical in behavior except that thread-based parallelism is used to accelerate the computation.

The optional chunksize argument provides a fixed upper bound for the number of simulations to be passed to each thread. Depending on number of independent categorical distributions to be sample and number of simulations to be performed, a chunksize of 5000 is likely too conservative. The user is encouraged to try smaller chunk sizes.

CategoricalMonteCarlo.vsample!Method
vsample!(B::AbstractArray, A::AbstractArray)

Identical to sample! except that except that the underlying Marsaglia sampler uses LoopVectorization.

CategoricalMonteCarlo.vsampleMethod
vsample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing])

See sample for full documentation; identical in behavior except that the underlying Marsaglia sampler uses LoopVectorization. Sometimes, but not always, confers speedup.

CategoricalMonteCarlo.vtsample!Method
tsample!(B::AbstractArray, A::AbstractArray; [chunksize=5000])

Identical to sample! except that thread-based parallelism is used to accelerate the computation. also, the underlying Marsaglia sampler uses LoopVectorization.

CategoricalMonteCarlo.vtsampleMethod
tsample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing], [chunksize=5000])

See sample for full documentation; identical in behavior except that thread-based parallelism is used to accelerate the computation; also, the underlying Marsaglia sampler uses LoopVectorization.

The optional chunksize argument provides a fixed upper bound for the number of simulations to be passed to each thread. Depending on number of independent categorical distributions to be sample and number of simulations to be performed, a chunksize of 5000 is likely too conservative. The user is encouraged to try smaller chunk sizes.