Wrong results for weighted quantile

[andreasnoack]

From https://discourse.julialang.org/t/median-vs-50th-quantile-giving-different-answers/18414/4

julia> x = [1; 4; 3; 2; 2.5; 7];

julia> w = [0.1;0.3;0.05;0.05;0.2;0.3];

julia> quantile(x, weights(w), 0.5)
3.5

but

julia> sum(w[x .<= 3.5]) >= 0.5
false

so 0.5 is not a median.

cc @matthieugomez

[matthieugomez]

There are multiple ways to extend the base function quantile to the case of weights, and I am happy to discuss how to do it properly.

But it is not obvious to me that the existing function has the wrong output. You reject that 3.5 is the median because you say: if there is only one value v in the vector that satisfies P(x <=v) >= p and P(x >= v ) >= 1 - p, then v should be the p quantile.

But this is not true with the definition of quantile in Julia. For instance:

x = [1, 2, 3, 4]
df= DataFrame(plow = cumsum([0.25, 0.25, 0.25, 0.25]), phigh = reverse(cumsum(reverse([0.25, 0.25, 0.25, 0.25]))), x = x)
#>│ Row │ plow    │ phigh   │ x     │
#>│     │ Float64 │ Float64 │ Int64 │
#>├─────┼─────────┼─────────┼───────┤
#>│ 1   │ 0.25    │ 1.0     │ 1     │
#>│ 2   │ 0.5     │ 0.75    │ 2     │
#>│ 3   │ 0.75    │ 0.5     │ 3     │
#>│ 4   │ 1.0     │ 0.25    │ 4     │

Julia says that the 0.49 quantile is 2.47 but, according to your definition, it should be 2.
So even though the property happens to be true for the unweighted 0.5 quantile, I am not sure it should be true for the weighted 0.5 quantile.

That being said, I think the thread shows two problems with the existing implementation:

• Whatever the definition of weighted quantile we choose in the end, wmedian should return the same output as wquantile(., 0.5). For now they return different things because they use different definitions of weighted medians.
• with frequency weigths, wquantile fails. This is because I wrote this function assuming that fweights were restricted to integer. How to interpret a frequency weight that is not an integer?

[matthieugomez]

Btw, note that the current definition of weighted median in StatsBase is wrong because the median with frequency weights does not give the same thing as the median of the vector with repeated values.

median([1, 2, 3, 4], fweights([2, 1, 3, 2]))
# 2.5
median([1, 1, 2, 3, 3, 3, 4, 4])
# 3

[matthieugomez]

At the end of the day, quantile type 7 is just not straighforward to extend to weights. I would be happier to have quantile defined as type 1 or 3 in base. It would make it simpler to extend it to weight. It would also make it possible to use quantile with types that are only ordinal. JuliaLang/julia#27367

[nalimilan]

Using the same default as R and NumPy is nice, though. For ordinal data, it sounds OK to require people to choose a different method.

• Whatever the definition of weighted quantile we choose in the end, wmedian should return the same output as wquantile(., 0.5). For now they return different things because they use different definitions of weighted medians.

Can't we just have it call quantile(x, 0.5)? Or would it be less efficient?

• with frequency weigths, wquantile fails. This is because I wrote this function assuming that fweights were restricted to integer. How to interpret a frequency weight that is not an integer?

Cf. preivous discussion here: #313 (comment). I guess several generalizations are possible, we just need one which is reasonable (for some definition of that term). Does the current answer (7.0) fit that bill?

[matthieugomez]

• quantile(x, 0.5) seems to be as (in)efficient as wmedian so there is no performance reason to have a separate implementation.

x = rand(10_000_000)
w = rand(10_000_000)
@time wquantile(x, w, 0.5)
  1.908000 seconds (44 allocations: 382.676 MiB, 2.88% gc time)
0.5001186970922967

julia> @time wmedian(x, w)
  1.911641 seconds (27 allocations: 231.275 MiB, 2.02% gc time)
0.5001187598955092

• For the second part, I just don't understand what non integer frequency weight mean. In Stata, it is impossible to have frequency weights that are not integer. Why is it possible in StatsBase? The previous discussion you linked to is about non-integer probability weights, which is fine.

[nalimilan]

For many algorithms there's an obvious generalization of frequency weights to non-integer values. For example, to compute the mean it doesn't matter whether the weight is integer or not. Generally speaking, it sounds simple to give a meaning to non-integer weights: having two observations with weight w/2 should be equivalent to having one observation with weight w. That said I'm not sure it's very useful (even if I remember somebody asked for them in another context), so if we're unsure we could throw an error in the presence of non-integer weights for now.

Regarding the weighted median, let's remove the separate implementation then.

[matthieugomez]

But, in your example, it should then be interpreted as a probability weight.I am having a hard time thinking of a case where someone wants non frequency weights to be different than probability weights. From Stata:

fweights, or frequency weights, are weights that indicate the number of duplicated observations.
pweights, or sampling weights, are weights that denote the inverse of the probability that the
        observation is included because of the sampling design.

[nalimilan]

I'm not saying they have to be different when computing quantiles. That's the case e.g. for inference of when computing the Bessel correction for variance, but not necessarily for descriptive statistics.

[matthieugomez]

While I am doing a pull request, why is wmean deprecated but not wsum or wmedian?

[nalimilan]

I guess that's just an oversight.

Do we have any idea why Julia uses a separate method for median rather than calling quantile? Just because it's simple? I'd just like to be sure we're not missing something.

[nalimilan]

So apart from #436, does something need to be addressed here @matthieugomez? We could reject non-integer frequency weights, but do you have an example where they give a clearly incorrect result?

[matthieugomez]

I don't really know what correct or incorrect result would mean in this situation, so I prefer to return an error. This is what Stata does (not sure about whether other softwares have a concept of frequency weights).

…

On Fri, Feb 1, 2019 at 10:33 AM Milan Bouchet-Valat < ***@***.***> wrote: So apart from #436 <#436>, does something need to be addressed here @matthieugomez <https://github.com/matthieugomez>? We could reject non-integer frequency weights, but do you have an example where they give a clearly incorrect result? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#435 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AF733f6O6AECWPKCokhZ_mcVzUsTIiasks5vJF5JgaJpZM4ZIayv> .

[nalimilan]

I don't really know what correct or incorrect result would mean in this situation,

For example, wouldn't the rule I describe above make sense? "having two observations with weight w/2 should be equivalent to having one observation with weight w".

[matthieugomez]

I'm not sure this rule would work. Remember the example:

quantile([1, 2], fweights([2, 2]), [0.25]) = 1
#but
quantile([1, 2], fweights([1, 1]), [0.25]) > 1

which violates your rule.