# The sampling without replacement

I have an EmpiricalDistribution

  d = EmpiricalDistribution[{1, 1, 1, 2, 2, 2, 3, 3, 3, 3}]


now, I want to sample twice without replacement $x_1,x_2$.

afterward, compute probability distribution of $|x_1-x_2|$ to get:

  X=|x_1-x_2|   0        1       2
P        4/15    7/15     4/15


What command should be used?

• d = {1,1,1,2,2,2,3,3,3,3}; samples=Subsets[d, {2}]; Map[{First[#], Length[#]/Length[samples]} &, Split[Sort[Map[Abs[First[#]-Last[#]]&, samples]]]] which gives you {{0, 4/15}, {1, 7/15}, {2, 4/15}}
– Bill
Commented Nov 30, 2016 at 1:00
• This is a brute force algorithm. Is there a little more simple method? Commented Nov 30, 2016 at 1:13
• be aware if you want to work with EmpiricalDistribution the second draw comes from a different distribution. Likely more complicated than Bills method Commented Nov 30, 2016 at 4:01
• Hi tiankonghewo If your problem has more than 3 types of numbers or you need more than 3 sums then the solution that Bill gave is going to be a lot more attractive to you than the one I posted below. Commented Nov 30, 2016 at 4:45
• Hi tiankonghewo; Always wait about 24 hours or so before accepting an answer. The first answer is not always the best and never is when it is mine. It is just to get the ball rolling and if you wait without accepting you will encourage others to answer. Commented Nov 30, 2016 at 6:15

You can do it this way but it is a bit clumsy:

For $\left | x2-x1 \right |=0$

Probability[x == 2 || y == 2 || z == 2, {x, y, z} \[Distributed]
MultivariateHypergeometricDistribution[2, {3, 3, 4}]]


For $\left | x2-x1 \right |=1$

Probability[(x == 1 && y == 1) || (y == 1 && z == 1), {x, y, z}
\[Distributed]MultivariateHypergeometricDistribution[2, {3, 3, 4}]]


For $\left | x2-x1 \right |=2$

Probability[(x == 1 && z == 1), {x, y, z}
\[Distributed]MultivariateHypergeometricDistribution[2, {3, 3, 4}]]

• Great! this is just what I need. I have learned something new. Thank you very much. Commented Nov 30, 2016 at 5:19
• You can write this as p[x_] = Piecewise[{Probability[#[[1]], {x, y, z} \[Distributed] MultivariateHypergeometricDistribution[2, {3, 3, 4}]], x == #[[2]]} & /@ {{x == 2 || y == 2 || z == 2, 0}, {(x == 1 && y == 1) || (y == 1 && z == 1), 1}, {x == 1 && z == 1, 2}}] Commented Nov 30, 2016 at 6:07
• Neater and cleaner, thanks +1 Commented Nov 30, 2016 at 6:16
• OK, Piecewise is a new thing to me. I learned something new. Commented Nov 30, 2016 at 7:18

I am a bit late but another way to use MultivariateHypergeometricDistribution:

md = MultivariateHypergeometricDistribution[2, {3, 3, 4}];
f[{___, 2, ___}] := 0
f[{1, 0, 1}] := 2
f[{___, 1, 1, ___}] := 1
td = TransformedDistribution[
f[{x, y, z}], {x, y, z} \[Distributed] md];
res = Probability[x == #, x \[Distributed] td] & /@ Range[0, 2]
Histogram[RandomVariate[td, 10000], Automatic, "PDF",
Epilog -> {Red, PointSize[0.02],
Point[MapIndexed[{#2[[1]] - 1/2, #1} &, res]]}]


• That is nice, I played with the Histogram command but could not make it work.+1 Commented Nov 30, 2016 at 6:21
• @bobbym thank you, your answer got to the guts of the matter (reciprocal +1). I just wanted to show another way to use MultivariateHypergeometricDistribution :) Commented Nov 30, 2016 at 6:23
• This method is more simple and convenient. Commented Nov 30, 2016 at 7:30
• @tiankonghewo please do not feel compelled to accept my answer. I think MultivariateHypergeometricDistribution was the key insight and bobbym and all our subsequent answers. I do thank you for the accept. Commented Nov 30, 2016 at 7:34

p[x_] = Piecewise[{Probability[#[[1]], {x, y, z} \[Distributed]
MultivariateHypergeometricDistribution[2, {3, 3, 4}]],
x == #[[2]]} & /@
{{x == 2 || y == 2 || z == 2, 0},
{(x == 1 && y == 1) || (y == 1 && z == 1), 1},
{x == 1 && z == 1, 2}}]


The distribution for |x1-x2| is then

dist = ProbabilityDistribution[p[x], {x, 0, 2, 1}];


This distribution can be used like any other distribution

PDF[dist, x] // Simplify


CDF[dist, x]


Mean[dist]

(*  1  *)

Variance[dist]

(*  8/15  *)

SeedRandom[1]

RandomVariate[dist, 10]

(*  {2, 0, 2, 0, 0, 0, 1, 0, 1, 1}  *)

• instructive amplification +1 (TransformedDistribution also works nicely in this case) Commented Nov 30, 2016 at 6:45

in case you wanted to work with EmpericalDistribution :

list = {1, 1, 1, 2, 2, 2, 3, 3, 3, 3};
outcomes = Tuples[ConstantArray[Union[list], 2]]


{{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 2}, {2, 3}, {3, 1}, {3, 2}, {3, 3}}

{Abs[Subtract @@ #], PDF[EmpiricalDistribution[list], #[[1]]]
PDF[EmpiricalDistribution[
Drop[list, First@Position[list, #[[1]]]]], #[[2]]] } & /@
outcomes;
{#[[1, 1]], Total[#[[All, 2]]]} & /@ GatherBy[%, First]


{{0, 4/15}, {1, 7/15}, {2, 4/15}}

of course , PDF[EmpiricalDistribution[list], i] just gives you Count[list,i]/Length@list