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.