CODEWITHSUNDEEP
randomjuliadistribution
Marouane1994
Marouane1994

Reputation: 534

Generating random values from a customized distribution

I have this distribution below:

using Distributions
struct OrthoNNDist <: DiscreteMultivariateDistribution
    x0::Vector{Int64}
    oc::Array{Int64,2}
    x1s::Array
    prob::Float64
    #return a new uniform distribution with all vectors in x1s orthogonal to oc
    function OrthoNNDist(x0::Vector{Int}, oc::Array{Int,2})
        x1s = []
        for i = 1:size(oc)[2]
            x1 = x0 + oc[:, i]
            if nonneg(x1)
                push!(x1s, x1)
            end
            x1 = x0 - oc[:, i]
            if nonneg(x1)
                push!(x1s, x1)
            end
        end
        new(x0, oc, x1s, 1.0/length(x1s))
    end
end

Base.length(d::OrthoNNDist) = length(d.x0)

Distributions.rand(d::OrthoNNDist, N::Integer=1) = rand(d.x1s, 1)

Distributions.pdf(d::OrthoNNDist, x::Vector) = x in d.x1s ? D.prob : 0.0
Distributions.pdf(d::OrthoNNDist) = fill(d.prob, size(d.x1s))
Distributions.logpdf(d::OrthoNNDist, x::Vector) = log(PDF(d, x))

and I want to generate random values from it, I don't know how I've tried:rand(OrthoNNDist,1000)and it didn't work, I'm kinda new on probabilistic programming, I don't know how I can do that.

Upvotes: 3

Views: 109

Answers (1)

Jeffrey Sarnoff
Jeffrey Sarnoff

Reputation: 1757

The function nonneg is no longer provided, that is easy enough to remedy:

nonneg(x::Real) = zero(x) <= x
nonneg(x::Vector{<:Real}) = all(nonneg, x)

Almost all of what you wrote is fine:

using Distributions

struct OrthoNNDist <: DiscreteMultivariateDistribution
    x0::Vector{Int64}
    oc::Array{Int64,2}
    x1s::Array
    prob::Float64
    #return a new uniform distribution with all vectors in x1s orthogonal to oc
    function OrthoNNDist(x0::Vector{Int}, oc::Array{Int,2})
        x1s = []
        for i = 1:size(oc)[2]
            x1 = x0 + oc[:, i]
            if nonneg(x1)
                push!(x1s, x1)
            end
            x1 = x0 - oc[:, i]
            if nonneg(x1)
                push!(x1s, x1)
            end
        end
        new(x0, oc, x1s, 1.0/length(x1s))
    end
end

Base.length(d::OrthoNNDist) = length(d.x0)

Distributions.pdf(d::OrthoNNDist, x::Vector) = x in d.x1s ? D.prob : 0.0
Distributions.pdf(d::OrthoNNDist) = fill(d.prob, size(d.x1s))
Distributions.logpdf(d::OrthoNNDist, x::Vector) = log(PDF(d, x))

To get rand working .. you are almost there

using Distributions: rand

Distributions.rand(d::OrthoNNDist, n::Int=1) = rand(d.x1s, n)

now with some data

julia> x0 = rand(1:1_000_000,5);
julia> oc = reshape(rand(1:1_000_000,25), (5,5));
julia> dist = OrthoNNDist(x0,oc);
julia> Distributions.rand(dist, 4)
4-element Array{Any,1}:
 [1330729, 656190, 927615, 470782, 1435138]
 [1382946, 1058057, 778316, 488440, 1304526]
 [1330729, 656190, 927615, 470782, 1435138]
 [1409093, 353679, 454229, 698320, 1271674]

(thanks to Mohamed Tarek for his direction)

Upvotes: 4

Related Questions