mean([1±1, 2±1]) should give a scalar

[mauro3]

Currently:

julia> mean([1±1, 2±1])
2-element Array{Float64,1}:
 1.0000000000000002
 2.0000000000000004

shouldn't this be what mean.([1±1, 2±1]) returns? I would expect it to return:

julia> ((1±1) + (2±1))/2
Particles{Float64,2000}
 1.5 ± 0.706

[baggepinnen]

The problem here is described in #88 (comment)
With the unfortunate fact that Distributions.mean === Statistics.mean, there is no way of distinguishing between mean of a vector and the mean of a multivariate distribution. A vector of particles is supposed to behave like both, but in some cases, it fails. I'm leaning towards introducing new functions dmean, dstd etc. that take the responsibility of what should have been Distributions.mean , but perhaps the issue will vanish if https://github.com/cscherrer/MeasureTheory.jl is used as backend instead of Distributions.jl? @cscherrer how does MeasureTheory treat functions like mean,var,std? Are they the same as those in Statistics or MeasureTheory defines its own?

[cscherrer]

Yeah this is not good:

julia> using Distributions
[ Info: Precompiling Distributions [31c24e10-a181-5473-b8eb-7969acd0382f]

julia> d = MvNormal(ones(3))
ZeroMeanDiagNormal(
dim: 3
μ: 3-element FillArrays.Zeros{Float64,1,Tuple{Base.OneTo{Int64}}} = 0.0
Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
)


julia> mean(d)
3-element FillArrays.Zeros{Float64,1,Tuple{Base.OneTo{Int64}}} = 0.0

julia> mean(rand(d,1000))
-0.018331363866233977

MeasureTheory doesn't have anything for this yet, but we should make sure we don't end up with strange behavior like this.

I think the core of the problem is that a sample from MvNormal produces a matrix. The semantics work better if this would instead be an array of arrays:

julia> mean(eachcol(rand(d,1000)))
3-element Array{Float64,1}:
 -0.007034913250783585
 -0.02964041764216485
 -0.028973082765793752

Distributions.jl made this choice for efficiency, but now it could make more sense to use something like ArraysOfArrays.jl.

@baggepinnen if MeasureTheory were to use Statistics.mean, but also have the semantics that for any measure d and n::Int,

typeof(mean(d)) == typeof(mean(rand(d, n)))

Would that fix things?

[mauro3]

More:

julia> mean(Any[0.411 ± 0.12, 0.342 ± 0.099])
Particles{Float64,2000}
 0.3765 ± 0.0763

julia> mean([0.411 ± 0.12, 0.342 ± 0.099])
2-element Array{Float64,1}:
 0.41099999999999987
 0.342

and

julia> std([1.0±1, 10±1])
Particles{Float64,2000}
 10.1 ± 0.994

[baggepinnen]

@baggepinnen if MeasureTheory were to use Statistics.mean, but also have the semantics that for any measure d and n::Int,
typeof(mean(d)) == typeof(mean(rand(d, n)))
Would that fix things?

I have a feeling that it's a step in the right direction, but might not completely solve the problem experienced here. The problem in MCM can be explained as p::Particles trying to act like two different types, a scalar and a distribution. Functions that are defined for both makes this hard to accomplish. In @mauro3 's first example, the vector is really a vector of vectors, and the expected result of mean(vec_of_vec) is the mean as if the element type of the outer vector was a scalar. Or put another way, the mean acted on the outer array and returned an array similar to the inner array type. In @cscherrer's proposed approach, a sample from a multivariate distributions is also a vector of vectors, with the mean acting on the outer array producing an element similar to the inner arrays.

There is still one problem though, the data storage layout between the two examples is transposed, and so their meaning is also transposed. @mauro3 is expecting the sample mean of the vector of particles, viewing the vector of particles not as a multivariate distributions but as a sample from a univariate distribution (which happens to be particle valued). The only way I can see out of his is to have two different functions, one that treats the data as coming from a multivariate distribution, and one that considers the data coming from a vector-valued univariate distribution (as weird as that sounds). I.e., one samplemean and one distributionmean. I hope that I am incorrect one someone finds a more elegant solution to this problem.

[cscherrer]

I haven't thought about this problem in a while in the MCM context, but I've been having some good luck lately with ArraysOfArrays.jl. Maybe that can help? Example:

julia> using Distributions

julia> using ArraysOfArrays

julia> d = MvNormal(ones(3))
ZeroMeanDiagNormal(
dim: 3
μ: 3-element FillArrays.Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}} = 0.0
Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
)


julia> mean(d)
3-element FillArrays.Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}} = 0.0

help?> nestedview
search: nestedview

  nestedview(A::AbstractArray{T,M+N}, M::Integer)
  nestedview(A::AbstractArray{T,2})

  AbstractArray{<:AbstractArray{T,M},N}

  View array A in as an M-dimensional array of N-dimensional arrays by wrapping it into an ArrayOfSimilarArrays.

  It's also possible to use a StaticVector of length S as the type of the inner arrays via

  nestedview(A::AbstractArray{T}, ::Type{StaticArrays.SVector{S}})
  nestedview(A::AbstractArray{T}, ::Type{StaticArrays.SVector{S,T}})

julia> mean(nestedview(rand(d,1000),1))
3-element Vector{Float64}:
  0.02048182849608377
 -0.03202973353056356
  0.01099568630469374

[mauro3]

But however this is resolved, my example above shows that mean and std behave inconsistently. Shouldn't they, in the mean time until this is resolved, both behave the same? Also, again in above example, that the result is dependent on the array-eltype is very un-julian. So, there are two issues here:

1. how to solve the bigger picture question: is a vector of Particles a sample from an n-dim multivariate distribution or n-samples from a Univariate distribution
2. how to make the current implementation consistent, sans solving 1

[baggepinnen]

They definitely should, and I believe the only way of making that happen is to let them behave like a vector of scalars rather than a MV distribution if there is a conflict. One can either supply new functions that treat them as an MV distribution, e.g, dmean, dstd, or provide a wrapper type ParticleDistribution <: Distribution

[mauro3]

But wouldn't it then make sense to wrap dimensions of MV distribution into something other than a Vector? Then you could keep the same function names. Also:

[...] to let them behave like a vector of scalars rather than a MV distribution if there is a conflict.

That seems like the "if there is a conflict" would be very confusing, no? Me, a stats-noob, I have no clear picture in my head which function makes sense for only MV. So likely I'd try around stuff with my vector of univariate distributions and then end up doing something which is totally inappropriate. But I guess, it's also ok for you to say that you not going to do any hand-holding with MCM.jl.

[baggepinnen]

Many of these problems should be fixed on master now, where Particles always behave like a scalar. Wrappers to make them behave like distributions are coming. To get the old behavior of mean, std, maximum etc., the functions pmean, pstd, pmaximum have been introduced

julia> mean([1±1, 2±1])
1.5 ± 0.71

julia> mean(Any[0.411 ± 0.12, 0.342 ± 0.099])
0.3765 ± 0.0792

julia> mean([0.411 ± 0.12, 0.342 ± 0.099])
0.3765 ± 0.078

julia> std([1.0±1, 10±1])
6.36396 ± 0.997

[baggepinnen]

A wrapper type around Particles that give them the semantics of a distributions is planned. I think all problems in this issue are resolved by today's changes.
Feel free to comment if they're not and I'll reopen :)