Inconsistency in MvNormal constructors, covariance or deviation

[dehann]

Hi,

I'm using Distributions.jl in IncrementalInference.jl and found what looks to be inconsistent behavior when constructing MvNormals. One dimensional cases take deviation while higher dimensions almost always take covariance.

using Distributions
julia> Base.std(rand(Normal(0.0,100.0), 1000))
96.7377564234093

julia> Base.std(rand(MvNormal([0.0],[100.0]), 1000), 2)
1×1 Array{Float64,2}:
 103.241

julia> Base.std(rand(MvNormal([0.0;0.0],[100.0;10.0]), 1000), 2)
2×1 Array{Float64,2}:
 101.026 
  10.0854

julia> Base.std(rand(MvNormal([0.0;0.0],[[100.0;0.0]';[0.0;10.0]']), 1000), 2)
2×1 Array{Float64,2}:
 10.1583 
  3.18971

Notice the last call expects covariance where all others take standard deviation. Do we want to persist this behavior? I don't really mind(think it is worth fixing), but it's a bit confusing and cost me some time debugging.

Thanks for putting this package out there!

[felixrehren]

This inconsistency also got me.

Mixing variance and volatility/deviation is in some sense type-instable for dimensional quantities (e.g. for market prices measured in USD, the deviation is in USD, the variance in USD^2).

It feels surprising that

• you can either provide a diagonal matrix, or the diagonal of the same matrix, but these two will construct different distributions
• univariate MvNormal is different from multivariate MvNormal (if the correlation is specified -- which is presumably the main reason to use multivariate normal distributions), which makes d=1 a very special case indeed

This is the number 2 result for google > "julia distributions", so it seems to strike a chord. But I have no idea how this could be changed without breaking people's code. Perhaps the documentation of the arguments should be put into all-caps?

[simonbyrne]

I agree it is a bit of a mess.

Currently we have univariate methods:

Normal() # N(0,1)
Normal(mean::Real) # N(mean,1)
Normal(mean::Real, stddev::Real)

for multivariate methods:

MvNormal(len::Int, stddev::Real) # zero mean
MvNormal(stddevvec::Vector) # zero mean
MvNormal(covar::Matrix) # zero mean
MvNormal(meanvec::Vector, stddev::Real)
MvNormal(meanvec::Vector, stddevvec::Vector)
MvNormal(meanvec::Vector, covar::Matrix)

Maybe keyword args are the way to go?

[dehann]

I would vote for not having special cases, and if this is going to be fixed it should be before Julia 1.0... Agreed that this would have deep affects in the package code base -- the only way would be clear documentation, announcements, and maximum lead-time warning upon using Distributions?

Warning:  Distributions.Normal and Distributions.MvNormal constructors are being standardized to deviation/covariance/?.  Julia 1.0 versions will use the new standard, i.e. Distributions v???+.

Maybe the Julia 0.7 cycle is long enough? It's not great, but rather now than waiting until Julia 2.0 or never at all. I would also make JuliaComputing aware of this and let them make the final call?

i was recently caught off guard by the DiffEqs and NLsolve API change that switched the order of function arguments -- but wasn't that bad once the fixes were in.

Lastly, I think it is okay to change, as long as the package tag numbers clearly indicate the transition without any other changes.

[dehann]

Keyword arguments are a good idea too. Maybe some combination?

[simonbyrne]

Having thought about this a bit more, I think all distributions should have:

• a single "canonical" constructor (e.g. Normal(mean, stddev))
• a keyword arg constructor, which may accept different ways to specify the distribution, and possible default values. e.g. Normal(mean=2) == Normal(2,1), Normal(var=9) == Normal(sigma=3) == Normal(0, 3).

[Affie]

The difference in these two constructors leads to potential errors, I had to test it to make sure of the usage:

julia> MvNormal([1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 3.0])
ZeroMeanFullNormal(
dim: 3
μ: [0.0, 0.0, 0.0]
Σ: [1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 3.0]
)


julia> MvNormal([1.0,2.0,3.0])
ZeroMeanDiagNormal(
dim: 3
μ: [0.0, 0.0, 0.0]
Σ: [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 9.0]
)

Perhaps the documentation can be a bit more explicit on this part:

3. vector of type Vector{T}: indicating a diagonal covariance as diagm(abs2(sig)),

Its easy to miss and perhaps more emphasis should be on the fact that it is the vector of the standard deviations for the diagonal, and not the diagonal of the covariance itself. So the abs2(sig) part.

or perhaps by including #584 (comment):

MvNormal(len::Int, stddev::Real) # zero mean
MvNormal(stddevvec::Vector) # zero mean
MvNormal(covar::Matrix) # zero mean
MvNormal(meanvec::Vector, stddev::Real)
MvNormal(meanvec::Vector, stddevvec::Vector)
MvNormal(meanvec::Vector, covar::Matrix)

[andreasnoack]

Recently, we've had issues with the difference between

julia> Normal(0.1)
Normal{Float64}(μ=0.1, σ=1.0)

julia> MvNormal([0.1;;])
ZeroMeanFullNormal(
dim: 1
μ: Zeros(1)
Σ: [0.1;;]
)

Within our application, the one-argument MvNormal constructor makes sense while the one-argument Normal constructor` does not. Generally, I don't really see how that constructor could ever be useful but before deprecating it, I'd like to ask if anybody here consider it useful?

[devmotion]

Generally, I wonder if keyword arguments are the better approach for distributions with many scalar parameters. An (now outdated) draft was #1405.