CODEWITHSUNDEEP
julia
Ben S.
Ben S.

Reputation: 3545

Make a probability distribution from two distributions in Julia

using Distributions
d1 = Exponential(0.2)
d2 = Exponential(0.5)
p = 0.7

Is there any easy way I construct a distribution in Julia, that behaves like a distribution in that I can call rand() and rand!, that behaves as follows: draw from distribution d1 with probability p and draw from distribution d2 with probability 1-p. Thank you.

Upvotes: 3

Views: 622

Answers (2)

mbauman
mbauman

Reputation: 31342

You can just use a MixtureModel:

julia> m = MixtureModel([d1,d2],[p,1-p])
MixtureModel{Distributions.Exponential{Float64}}(K = 2)
components[1] (prior = 0.7000): Distributions.Exponential{Float64}(θ=0.2)
components[2] (prior = 0.3000): Distributions.Exponential{Float64}(θ=0.5)

julia> mean(m)
0.29000000000000004

julia> pdf(m, 0)
4.1

julia> rand(m)
0.2574516697519676

julia> rand!(m, zeros(1,5))
1×5 Array{Float64,2}:
 0.0704624  0.264519  0.636179  0.11479  0.41158

Upvotes: 4

Dan Getz
Dan Getz

Reputation: 18217

Distributions.jl basically prepared all the tools to define new distributions. So, in this case, my attempt looks like:

using Distributions

struct CompoundBernoulli{T<:Distributions.VariateForm,
                         S<:Distributions.ValueSupport} <:
         Distributions.Sampleable{T, S}
    p::Bernoulli
    d1::Distributions.Sampleable{T,S}
    d2::Distributions.Sampleable{T,S}
end

# outer constructor  
CompoundBernoulli(p::Real, 
                  d1::Distributions.Sampleable{S, T},
                  d2::Distributions.Sampleable{S, T}) where 
  {S<:Distributions.VariateForm, T<:Distributions.ValueSupport} = 
  CompoundBernoulli{S,T}(Bernoulli(p),d1,d2)

Base.rand(cb::CompoundBernoulli) = rand(cb.p)==0 ? rand(cb.d1) : rand(cb.d2)

With these definitions:

julia> cb = CompoundBernoulli(0.7,Exponential(0.2),Exponential(0.5))
CompoundBernoulli{Distributions.Univariate,Distributions.Continuous}
(Distributions.Bernoulli{Float64}(p=0.7), 
Distributions.Exponential{Float64}(θ=0.2), 
Distributions.Exponential{Float64}(θ=0.5))

julia> rand(cb)
0.3247418465183849

julia> rand(cb,3,3)
3×3 Array{Float64,2}:
 0.33105   0.231418  0.271571
 0.413905  0.662144  1.42725 
 0.20196   0.091628  0.194761

More functions can be defined and functions can be specialized for this specific type as the application requires.

Upvotes: 2

Related Questions