In the below, I'd have expected either rand(d, 2) to give a vector of vectors or rand(d, 2, 4) to give a 3x2x4 array, but currently it seems that the behaviour of these methods is inconsistent. I think this is true of all multivariate distributions (I've tested a handful).
For matrix distributions it always returns array-of-matrix regardless of how many extra dimensions are added.
julia> using Distributions; d = Dirichlet(ones(3));
julia> rand(d)
3-element Vector{Float64}:
0.07995751567688866
0.7342415416771756
0.18580094264593575
julia> rand(d, 2)
3×2 Matrix{Float64}:
0.154931 0.544246
0.323862 0.29087
0.521207 0.164884
julia> rand(d, 2, 4)
2×4 Matrix{Vector{Float64}}:
[0.1947, 0.197046, 0.608253] [0.155361, 0.826544, 0.0180956] … [0.157334, 0.764753, 0.0779129]
[0.305394, 0.605478, 0.089128] [0.524005, 0.280374, 0.195621] [0.17374, 0.612635, 0.213625]
In the below, I'd have expected either
rand(d, 2)to give a vector of vectors orrand(d, 2, 4)to give a 3x2x4 array, but currently it seems that the behaviour of these methods is inconsistent. I think this is true of all multivariate distributions (I've tested a handful).For matrix distributions it always returns array-of-matrix regardless of how many extra dimensions are added.