I want to use a custom multivariate distribution that actually generates arrays, and I want to know if there’s anything I need to modify in order to be adapted to be used in Turing
using Distributions
struct OrthoNNDist <: DiscreteMultivariateDistribution
x0::Vector{Float64}
oc::Array{Float64,2}
x1s::Array
prob::Float64
#return a new uniform distribution with all vectors in x1s orthogonal to oc
function OrthoNNDist(x0::Vector{Float64}, oc::Array{Float64,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) = rand(d.x1s)
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))
I don’t know if it’s fit to be used in Turing or no.
I don’t think your code will work with automatic differentiation. The problem is that your types need to be a sub-type of real.
This might get you started
struct OrthoNNDist{T<:Real} <: DiscreteMultivariateDistribution
x0::Vector{T}
oc::Array{T,2}
x1s::Array{T}
prob::T
#return a new uniform distribution with all vectors in x1s orthogonal to oc
function OrthoNNDist(x0::Vector{T}, oc::Array{T,2}) where {T<:Real}
x1s = T[]
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
2 Likes
i wanted to try it, but then it shows me:
ERROR: syntax: too few type parameters specified in "new{...}" around REPL[25]:1
Stacktrace:
[1] top-level scope
@ REPL[25]:1
Can you please provide the code that generated that error?
Ok. I see that you were trying to run the code for the constructor. You have two options. 1) you can pass the type information to new as new{T}, 2) you can pull your constructor outside of the struct definition. I personally prefer the second option because its easier to read.
using Distributions
struct OrthoNNDist{T<:Real} <: DiscreteMultivariateDistribution
x0::Vector{T}
oc::Array{T,2}
x1s::Array{T}
prob::T
end
#return a new uniform distribution with all vectors in x1s orthogonal to oc
function OrthoNNDist(x0::Vector{T}, oc::Array{T,2}) where {T<:Real}
x1s = T[]
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
return OrthoNNDist(x0, oc, x1s, 1.0/length(x1s)::T)
end
1 Like
A few notes that might be helpful. In the return statement above, I annotated the last argument as 1.0/length(x1)::T, to ensure that prob assumes type T.
The code assumes that all fields will be an array with elements of type T, or type T itself. If that is too restrictive, you can define T1,T2, … to make allow fields to have different subtypes of Real.
1 Like
I tried the last code you put in the comment and I took these arguments to run the code: oc=[3.0 2.0 1.0; 0.0 1.0 1.0], x0=[1.0;2.0]
and i got this:
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float64
Closest candidates are:
convert(::Type{T}, ::Base.TwicePrecision) where T<:Number at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/twiceprecision.jl:262
convert(::Type{T}, ::AbstractChar) where T<:Number at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/char.jl:185
convert(::Type{T}, ::CartesianIndex{1}) where T<:Number at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/multidimensional.jl:136
...
Stacktrace:
[1] push!(a::Vector{Float64}, item::Vector{Float64})
@ Base ./array.jl:994
[2] OrthoNNDist(x0::Vector{Float64}, oc::Matrix{Float64})
@ Main ./REPL[19]:7
[3] top-level scope
@ REPL[20]:1
There were a few issues that I fixed. First, upon closer inspection, it is clear that x1 is a vector of vectors. I updated the code to reflect that. Also, I introduced an error in defining probs, which is now fixed. The function nonneg was not defined, so I added a definition. Finally, size(oc)[2] makes an unnecessary memory allocation. The changing it to size(oc, 2) fixes that minor performance issue.
The code below runs without an error. It should work with Turing as long as all elements should be of type T.
Summary
using Distributions
nonneg(x) = x >= zero(x)
struct OrthoNNDist{T<:Real} <: DiscreteMultivariateDistribution
x0::Vector{T}
oc::Array{T,2}
x1s::Vector{Vector{T}}
prob::T
end
#return a new uniform distribution with all vectors in x1s orthogonal to oc
function OrthoNNDist(x0::Vector{T}, oc::Array{T,2}) where {T<:Real}
x1s = Vector{Vector{T}}(undef,0)
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
prob::T = 1.0 / length(x1s)
return OrthoNNDist(x0, oc, x1s, prob)
end
oc = [3.0 2.0 1.0; 0.0 1.0 1.0]
x0 = [1.0;2.0]
dist = OrthoNNDist(x0, oc)
1 Like
yeah concerning the nonneg function i forgot to put it down with the code it’s defined like this: nonneg(v) = all(v.>=0) ? true : false where v is a vector
If you’re going to use it with samplers such as HMC, you also need to
define bijector for your distribution.
I recommend checking out
https://turing.ml/v0.21/docs/using-turing/advanced#how-to-define-a-customized-distribution
for more information regarding these things.
1 Like
i’m just trying to use a metropolis-hastings algorithm to begin with to see how it goes.
Depending on which MH approach you’re using, this might still be
required. For example, Random walk MH, i.e. MH with MvNormal as
proposal assumes the space is unconstrained and so requires impl of
bijector.
If you’re specifying your own proposal which statisfies the constraints
of the variables, then AFAIK you should be fine even without bijector.
But I’d still recomemend just implementing it; shouldn’t be too difficult.