# Custom multivariate distribution in Turing.jl

**URL:** <https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [October 7, 2022, 6:12pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414 "2022-10-07T18:12:06Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![marouane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marouane/32/16267_2.png) [@marouane](https://discourse.julialang.org/u/marouane)\
**Post date:** [October 7, 2022, 6:12pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/1 "2022-10-07T18:12:06Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 7, 2022, 6:26pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/2 "2022-10-07T18:26:26Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![marouane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marouane/32/16267_2.png) [@marouane](https://discourse.julialang.org/u/marouane)\
**Post date:** [October 8, 2022, 1:55pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/3 "2022-10-08T13:55:37Z")

</div>

i wanted to try it, but then it shows me:

```julia
 ERROR: syntax: too few type parameters specified in "new{...}" around REPL[25]:1
Stacktrace:
 [1] top-level scope
   @ REPL[25]:1

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 8, 2022, 4:22pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/4 "2022-10-08T16:22:24Z")

</div>

Can you please provide the code that generated that error?

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 8, 2022, 4:30pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/5 "2022-10-08T16:30:02Z")

</div>

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.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 8, 2022, 4:35pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/6 "2022-10-08T16:35:19Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![marouane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marouane/32/16267_2.png) [@marouane](https://discourse.julialang.org/u/marouane)\
**Post date:** [October 10, 2022, 11:44am UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/7 "2022-10-10T11:44:34Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [October 10, 2022, 11:59am UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/8 "2022-10-10T11:59:56Z")

</div>

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**
>
> ```julia
> 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![marouane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marouane/32/16267_2.png) [@marouane](https://discourse.julialang.org/u/marouane)\
**Post date:** [October 10, 2022, 1:41pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/9 "2022-10-10T13:41:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [October 10, 2022, 9:15pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/10 "2022-10-10T21:15:46Z")

</div>

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](https://turing.ml/v0.21/docs/using-turing/advanced#how-to-define-a-customized-distribution)  
for more information regarding these things.

---

<div class="post-metadata">

**Author:** ![marouane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marouane/32/16267_2.png) [@marouane](https://discourse.julialang.org/u/marouane)\
**Post date:** [October 13, 2022, 9:31am UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/11 "2022-10-13T09:31:37Z")

</div>

i’m just trying to use a metropolis-hastings algorithm to begin with to see how it goes.

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [October 13, 2022, 1:31pm UTC](https://discourse.julialang.org/t/custom-multivariate-distribution-in-turing-jl/88414/12 "2022-10-13T13:31:35Z")

</div>

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.