How to get median based on probability distribution?

From some calculation, I have a distribution of discrete data $$(i,P_i)$$ and then want to get the median based on this distribution. Naive way is to create a list

list = Table[Sum[p[[j]], {j, i}], {i, n}]


then find the first element greater than $$0$$:

pos = Position[list, _?(# >= 0.5 &)][[1, 1]]


but this way is pretty slow, especially when I have a very large amount of data. The reason is I have to do Sum everytime. I did several attempts to speed up, such as

list = Table[0, {i, n}]
For[i = 1, i <= n, i ++, list[[i]] = If[i == 1, p[[i]], list[[i - 1]] + p[[i]]]]


or even smarter by using Accumulate

Accumulate[p]


and then do same Position operation. This made everything much faster and I'm pretty happy with it. I'm wondering whether some similar function is already in Mathematica, so we don't have to manually implement this. However after lookup the Median, I have no results relate to this. Do you guys have any idea?

• Have you seen EmpiricalDistribution[]? – J. M. is away Oct 15 '18 at 17:44
• @J.M.iscomputer-less thank you! – RoderickLee Oct 15 '18 at 18:04

To elaborate a bit on what J.M. hinted at, this is one way of achieving what you want with EmpiricalDistribution.

First let's get an example table of pairs of {value,probability} like you showed in your question

list = Transpose[{Range, #/Total[#] &[RandomReal[{0, 1}, 10]]}]


{{1, 0.0538923}, {2, 0.00538521}, {3, 0.158895}, {4, 0.10697}, {5, 0.0799713}, {6, 0.17624}, {7, 0.112601}, {8, 0.128191}, {9, 0.156779}, {10, 0.0210756}}

Then we make this into a EmpiricalDistribution:

dist = EmpiricalDistribution[list[[All, 2]] -> list[[All, 1]]]
Plot[CDF[dist, x], {x, 0, 11}, Filling -> Axis, Exclusions -> None] Here we use the syntax where we have used the probabilities as weights to single sample examples to get the right distribution. If you have your original data sample before binning that's even better as an input and EmpiricalDistribution will do the binning for you.

Now we can easily get the median by calling Median on our distribution:

Median[dist]


6

• Thanks for following through. ;) One thing: you could have used Normalize[RandomReal[{0, 1}, 10], Total] as well. – J. M. is away Oct 15 '18 at 18:00
• @J.M. Ah, cool idea! Just thought about #/Norm[#,1]&, but haven't thought about using Normalize in that way! – Thies Heidecke Oct 15 '18 at 18:11

You can also use WeightedData as follows:

SeedRandom
list = Transpose[{Range, #/Total[#] &[RandomReal[{0, 1}, 10]]}];
wd = WeightedData @@ Transpose[list];


You can use Median or Quantile with wd

Median[wd]


5

Quantile[wd, 1/2]


5

Alternatively,

Median[EmpiricalDistribution[wd]]
Quantile[EmpiricalDistribution[wd], 1/2]
InverseCDF[EmpiricalDistribution[wd], 1/2]


5

If it is really about performance, here is some inspiration.

Some data...

n = 10000000;
i = Range[n];
SeedRandom;
p = Normalize[RandomReal[{0, 1}, n], Total];


Thies' solution (with timing):

Median@EmpiricalDistribution[p -> i] // RepeatedTiming


{0.332, 4999860}

klgr's solution:

{0.18, 4999860}

Median[WeightedData @@ {i, p}] // RepeatedTiming


Computing the CDF and using Nearest to determine the elements closed to 0.5:

With[{P = Accumulate[p]},
If[P[[#]] < 0.5, # + 1, #] &@Last@Nearest[P -> i, 0.5]
] // RepeatedTiming


{0.046, 4999860}

Computing the CDF and using good ol' BinarySearch to determine the element closest to 0.5 (breaking ties in favor of the largest element below 0.5):

Needs["Combinatorica"]
With[{P = Accumulate[p]},
Ceiling@BinarySearch[P, 0.5]
] // RepeatedTiming


{0.023, 4999860}

Using OP's idea in a CompiledFunction (breaking ties in favor of the largest element below 0.5):

First compile

cf = Compile[{{p, _Real, 1}},
Block[{P = 0., i = 0},
While[P < 0.5, i++; P += CompileGetElement[p, i];];
i
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];


and then compute

i[[cf[p]]] // RepeatedTiming


{0.00568, 4999860}

• I don't think that the method with Nearest always gives the correct value for the median with a discrete distribution. Consider i = Range; p = {0.1, 0.39, 0.3, 0.01, 0.2}; Nearest[Accumulate[p] -> i, 0.5]. The result is 2 when the median is 3. – JimB Oct 16 '18 at 4:03
• @JimB Thanks for the hint. I was in the wrong assumption that the median were defined as the element in the list with cummulated probability closest to 0.5. Instead, we have to look for the smallest element with cummulated probability above to 0.5, right? I fixed the code, also for the binary searches. I hope that it is correct now. – Henrik Schumacher Oct 16 '18 at 6:50
• A median (rather than "the" median) for a univariate probability distribution is any value such that at least half the probability is below that value and at least half of the probability is above that value. For some discrete distributions that could result in 2 possible values that satisfy that definition. – JimB Oct 16 '18 at 13:06