# Improving Distributions.jl type consistency

**URL:** https://discourse.julialang.org/t/improving-distributions-jl-type-consistency/30082
**Category:** Statistics
**Created:** [October 19, 2019, 2:24pm UTC](https://discourse.julialang.org/t/improving-distributions-jl-type-consistency/30082 "2019-10-19T14:24:50Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)
#### Post date: [October 19, 2019, 2:24pm UTC](https://discourse.julialang.org/t/improving-distributions-jl-type-consistency/30082/1 "2019-10-19T14:24:50Z")

</div>

The `rand(dist, n::Int)` method is really handy, but there’s a lot of irregularity in the type that’s returned. Ideally we’d have some type `A` with

```julia
typeof(rand(dist, n::Int)) == A{typeof(rand(dist)), n}

```

Having an identity like this at the type level can allow us to build very general code, rather than relying no special cases. The obvious choice here is `A == Array`, but this loses some efficiency; for example, a regular array-of-vectors is more efficiently stored as a Matrix.

Currently, the relation is not a type constructor, but a _function_, something like

```julia
typeof(rand(dist, n::Int)) == f(typeof(rand(dist)), n)

f(Array{T, n}) = Array{T, n+1}
f(T) = Vector{T}

```

[RecursiveArrayTools.jl](https://github.com/JuliaDiffEq/RecursiveArrayTools.jl) gives a way to treat, for example, a vector-of-vectors as a matrix. It seems in the above case we need something to go the other direction, allowing internal representation as a regular array, but indexing as nested arrays.

Or maybe this is entirely wrong, and we need a different approach. Any ideas?

---

<div class="post-metadata">

### Author: ![johnczito](https://avatars.discourse-cdn.com/v4/letter/j/53a042/32.png) [@johnczito](https://discourse.julialang.org/u/johnczito)
#### Post date: [October 21, 2019, 5:50am UTC](https://discourse.julialang.org/t/improving-distributions-jl-type-consistency/30082/2 "2019-10-21T05:50:44Z")

</div>

> there’s a lot of irregularity in the type that’s returned

I can think of the fact that `rand(d::UnivariateDistribution, n::Int)` returns a vector and `rand(d::MultivariateDistribution, n::Int)` returns a matrix, but `rand(d::MatrixDistribution, n::Int)` returns a vector of matrices. What else do you have in mind?

(This came up [here](https://github.com/JuliaStats/Distributions.jl/issues/981), at least.)

---

<div class="post-metadata">

### Author: ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)
#### Post date: [October 21, 2019, 1:49pm UTC](https://discourse.julialang.org/t/improving-distributions-jl-type-consistency/30082/3 "2019-10-21T13:49:13Z")

</div>

Hi @johnczito, welcome!

The problem is the lack of a consistent relationship between `typeof(rand(dist))` and `typeof(rand(dist, n::Int))`. Starting with `dist = Normal()`, we have

```julia
julia> rand(Normal())
1.0200255085620293

julia> rand(Normal()) |> typeof
Float64

julia> rand(Normal(),2)
2-element Array{Float64,1}:
 -0.5708464683263881
 -0.5539728783723908

```

From this, one might expect that an extra `Int` argument to `rand` always changes the type from `T` to `Array{T,1}`. But that’s not it at all:

```julia
julia> rand(MvNormal(zeros(3)))
3-element Array{Float64,1}:
 -0.0
  0.0
  0.0

julia> rand(MvNormal(zeros(3)),2)
3×2 Array{Float64,2}:
 -0.0 0.0
 -0.0 0.0
  0.0 0.0

```

This might not seem like a big deal. After all, it’s easy enough to remember to treat these cases separately. The pain point comes when trying to write code that deals with an arbitrary distribution. As it is, generic functions on distributions behave strangely:

```julia
julia> samplemean(dist,n) = mean(rand(dist,n))
samplemean (generic function with 1 method)

julia> samplemean(Normal(),100)
0.13491200393928363

julia> samplemean(MvNormal(zeros(3)),100)
0.0

```
