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
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
typeof(rand(dist, n::Int)) == f(typeof(rand(dist)), n)
f(Array{T, n}) = Array{T, n+1}
f(T) = Vector{T}
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?
1 Like
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, at least.)
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> 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> 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> 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