# Custom multivariate distribution in Turing.jl

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
required. For example, Random walk MH, i.e. MH with `MvNormal` as
`bijector`.
of the variables, then AFAIK you should be fine even without `bijector`.