I have a population of 12 items, let's call them A through L. I'm going to be drawing 6 of those items out of the population without replacement. So you can't choose the same item twice.

What I know is that I want A to be selected as any one of the 6 items approximately 99% of the time, F to be selected as any one of the 6 items approximately 57% of the time, etc.. I'll post the full values below.

I'm using Python's NumPy library to sample without replacement, but I believe the methodology is similar to R. As each item is chosen from the population, its weight is distributed among the remaining items and then set to 0.

Kudos to @Glen_b for the example:

"If we're selecting 2 items from a pool of 3 items with first draw probabilities (0.1,0.4,0.5) and we select the third item first, the probabilities for the second item will be (0.2,0.8,0)"

What I'm trying to deduce is how to reverse engineer what weights would get sent to a sampling function to end up with the results below.

# Probabilities that each item is one of the 6 selected out of 12  

itemSelected = {
A = 0.99
B = 0.97
C = 0.95
D = 0.82
E = 0.68
F = 0.57
G = 0.48
H = 0.33
I = 0.11
J = 0.09
K = 0.01
L = 0.00
  • 1
    $\begingroup$ Since draws are without replacement, the probabilities will be different on each draw. So, could you explain more specifically what you mean by "the individual probabilities"? $\endgroup$ – user20160 Oct 14 '19 at 3:02
  • $\begingroup$ Presumably the intent is that there's an implicit set of weights, $w_i$ such that $p_i=w_i/\sum_i w_i$ and when an item is selected it can't be drawn again (its weight is set to 0). This is exactly how R's sample works when sampling without replacement. e.g. if we're selecting 2 items from a pool of 3 items with first draw probabilities $(0.1,0.4,0.5)$ and we select the third item first, the probabilities for the second item will be $(0.2,0.8,0)$ (this works for every set of weights that give the same $p_i$). However, if that's the intent the OP should make it explicit in the question. $\endgroup$ – Glen_b -Reinstate Monica Oct 14 '19 at 3:14
  • $\begingroup$ You should probably indicate how you know that a solution to your set of conditions is possible. Is it because you have observed what you're trying to fit? $\endgroup$ – Glen_b -Reinstate Monica Oct 14 '19 at 3:16
  • $\begingroup$ @user20160 My terminology may be poor, I meant that I'm seeking the probability that each item started with in order to end up with the observed results in the first list of numbers. The second list of numbers is my best guess at what they all started at, by simple trial and error. But I'd like a more repeatable method. Thank you! $\endgroup$ – Justin Geoghegan Oct 14 '19 at 4:12
  • $\begingroup$ @Glen_b Yes exactly, I'm trying to find a method to deduce what the set of weights are that are leading to the observed results in my first list, when selecting 6 items from the population of 12. I'm using python to sample without replacement but it sounds like the methodology is the same as R. I will edit my post to try to find some clarity. $\endgroup$ – Justin Geoghegan Oct 14 '19 at 4:17

At first, I thought a solution might require summing over all possible sequences of draws, with a resulting combinatorial explosion. But, thinking about the problem from a slightly different angle, there's actually an efficient way to approximate the solution. It only requires solving a simple system of linear equations.

Assuming the target probabilities given in the question and 6 draws without replacement, I calculated the item weights using the method described below. The weights for the first 11 items are:

$$w = [.2864, .2181, .1863, .1066, .0709, .0525, .0407, .0249, .0072, .0059, .0006]$$

The 12th item has zero weight and can be ignored, since it has zero probability of being drawn. Then, I simulated random draws using these weights, and calculated the empirical probabilities. They match the target probabilities fairly well:

enter image description here

Formulating the problem

Suppose there are $n$ unique items (which I'll refer to by number instead of letters). Let $w = [w_1, \dots, w_n]$ be weights for each item, which are assumed to be nonnegative and sum to one. Suppose we draw $k$ items without replacement. On the first draw, the probability of each item is given by its weight. On subsequent draws, the probability of each item is given by its weight divided by the summed weights of all remaining items (or zero if the the item has already been drawn).

Let $X = [X_1, \dots, X_n]$ be a binary random vector indicating whether each item has been drawn after $k$ draws; $X_i = 1$ if item $i$ has been drawn, otherwise 0. Let $\mu = [\mu_1, \dots, \mu_n]$ denote the expected value of $X$. Then, $\mu_i$ can be interpreted as the probability that item $i$ is included in the final set of $k$ selected items. Therefore, the problem can be framed as follows: Find the weights $w$ that give rise to a specified set of values in $\mu$ (target probabilities).

The relationship between $w$ and $\mu$

The distribution of $X$ can be identified by drawing an analogy to a paricular ball-and-urn problem. Let each item correspond to a weighted, colored ball in an urn. Each ball is uniquely colored, so there are $n$ colors with 1 ball each. We draw $k$ balls from the urn without replacement, according to their weights $w$ (as in the data generating process described above). $X$ can be interpreted as the number of times we observe each color. Then, the distribution of $X$ is a special case of the multivariate version of Wallenius' noncentral hypergeometric distribution. This is a relatively obscure, but fascinating distribution; see the references below for details. The general form allows for multiple balls of each color, but things simplify slightly in our case, since balls are uniquely colored.

$\mu$ is the mean of this distribution. Unfortunately, it doesn't have a closed form expression. But, Manly (1974) proposed an approximation (which is also described clearly in Fog 2008). Under this approximation (in our simplified case), $\mu$ and $w$ are related as:

$$(1-\mu_1)^{\frac{1}{w_1}} = (1-\mu_2)^{\frac{1}{w_2}} = \cdots = (1-\mu_n)^{\frac{1}{w_n}}$$

Solving for the weights

Take the reciprocal of the log of all sides:

$$\frac{w_1}{\log (1-\mu_1)} = \frac{w_2}{\log (1-\mu_2)} = \cdots = \frac{w_n}{\log (1-\mu_n)}$$

This is a system of $n (n-1) / 2$ linear equations. That is, for each unique $(i,j)$ pair:

$$\frac{w_i}{\log (1-\mu_i)} - \frac{w_j}{\log (1-\mu_j)} = 0$$

We also know that $w$ must sum to one, which can be expressed as an additional linear equation:

$$\sum_i w_i = 1$$

The task is to solve these equations for $w$, given $\mu$. This can be done by writing the entire system in the form $A w = b$, where $A$ is a matrix (constructed using the given $\mu$) and $b$ is a vector. $A$ and $b$ can then be passed to a standard linear solver to obtain a numerical solution for $w$.

Note that elements of the given $\mu$ must sum to $k$ by the definition of the problem. Furthermore, they must be strictly positive; otherwise $\log(0)$ would appear in the above equations. This isn't a big deal, since $\mu_i=0$ would imply that item $i$ has zero probability of being drawn, and we might as well drop it from the problem.


Manly (1974). A Model for Certain Types of Selection Experiments.

Fog (2008). Calculation Methods for Wallenius' Noncentral Hypergeometric Distribution.

  • $\begingroup$ This is above and beyond the answer I expected, thank you so much. As a novice in this field of study, I understand the description of the process all the way up to the solution. In this scenario, what does the matrix A and the vector b look like? $\endgroup$ – Justin Geoghegan Oct 17 '19 at 3:40
  • 1
    $\begingroup$ @JustinGeoghegan Glad to help. This was an interesting problem to think about; I'm happy it makes sense to you. Each row of $A$ and the corresponding element of $b$ represent an equation. To represent $\frac{w_i}{\log(1-\mu_i)} - \frac{w_j}{\log(1-\mu_j)} = 0$, set the $i$th column of the current row of $A$ to $\frac{1}{\log(1-\mu_i)}$, set the $j$th column to $-\frac{1}{\log(1-\mu_j)}$, and set all other columns to 0. Then, set the corresponding element of $b$ to 0. To represent $\sum_i w_i=1$, set all columns of the current row of $A$ to 1, and set the corresponding element of $b$ to 1. $\endgroup$ – user20160 Oct 17 '19 at 4:47
  • $\begingroup$ For the question example, how many rows of A are there? I think the answer is 66 because of 12(12-1)/2. Or I suppose maybe 55 if we throw out the one with zero probability and it's 11(11-1)/2. And then I build the matrix so it populates the first two columns of every row to be a unique i,j pair. Am I tracking so far? $\endgroup$ – Justin Geoghegan Oct 17 '19 at 18:59
  • 1
    $\begingroup$ @JustinGeoghegan For $n=11$ items, $A$ has 11 columns and 56 rows (55 rows for the $(i,j)$ pairs, plus an extra row for the sum-to-one equation). For each $(i,j)$ pair, you would use columns $i$ and $j$, not the first two columns. $\endgroup$ – user20160 Oct 18 '19 at 6:16
  • 1
    $\begingroup$ @JustinGeoghegan The calculations in your last comment look correct (assuming base 10 logs; any base will give the same final answer, although $A$ will contain different values). Keep in mind that the solution is an approximation. The quality of the approximation will be lower for small values of $n$. Looking through this answer again, I think there's actually a closed form solution; solving numerically is not necessary. I'll edit soon to include details. $\endgroup$ – user20160 Oct 23 '19 at 7:30

Initially, and even after some of the exchanges in Comments, I'm not sure I understand the question. Is this anything like what you're looking for? If not, maybe there are ideas here you can use. Everything below is in R:

pct = c(99,97,95,82, 68,57,48,33, 11,9,1,0)
prob = pct/sum(pct); prob
[1] 0.165000000 0.161666667 0.158333333 0.136666667 0.113333333 0.095000000
[7] 0.080000000 0.055000000 0.018333333 0.015000000 0.001666667 0.000000000
[1] 1
prob.4 = round(prob, 4); prob.4
[1] 0.1650 0.1617 0.1583 0.1367 0.1133 0.0950
[7] 0.0800 0.0550 0.0183 0.0150 0.0017 0.0000
[1] 1

Then for sampling:

 [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L"
sort(sample(LETTERS[1:12], 6, p=prob))
[1] "A" "B" "D" "E" "F" "G"
sort(sample(LETTERS[1:12], 6, p=prob))
[1] "A" "B" "C" "D" "E" "G"
sort(sample(LETTERS[1:12], 6, p=prob))
[1] "A" "B" "C" "D" "E" "F"
sort(sample(LETTERS[1:12], 6, p=prob))
[1] "A" "B" "C" "D" "I" "K"
sort(sample(LETTERS[1:12], 6, p=prob))
[1] "A" "B" "C" "E" "F" "H"

Below pick is a $6 \times 100$ matrix; t(pick) is its $100 \times 6$ transpose, and head(t(pick)) shows the first six rows.

pick = replicate(100, sort(sample(LETTERS[1:12], 6, p=prob)))
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] "A"  "B"  "C"  "F"  "G"  "H" 
[2,] "A"  "C"  "D"  "E"  "G"  "J" 
[3,] "A"  "B"  "C"  "D"  "H"  "J" 
[4,] "A"  "B"  "C"  "D"  "F"  "G" 
[5,] "A"  "C"  "D"  "F"  "G"  "H" 
[6,] "A"  "B"  "C"  "D"  "F"  "J" 

Note: Actually, if you're going to use the sample procedure in R, there is no need for its parameter p to be a probability vector. If it's not, R scales p as in the first few steps I showed.

When some letters have already selected, the selection probabilties among the remaining ones are chosen proportionately.


pp= cumsum(pct[12:1]/sum(pct))[12:1]
[1] 0.2512562814 0.2097989950 0.1691792295 0.1293969849
[5] 0.0950586265 0.0665829146 0.0427135678 0.0226130653
[9] 0.0087939698 0.0041876047 0.0004187605 0.0000000000

pick = replicate(10^5, sample(LETTERS[1:12], 6, p=pp))
head(t(pick))     # shows order in which letters chosen
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] "D"  "A"  "C"  "E"  "B"  "G" 
[2,] "C"  "G"  "A"  "D"  "E"  "B" 
[3,] "H"  "A"  "B"  "F"  "C"  "D" 
[4,] "B"  "C"  "A"  "F"  "E"  "G" 
[5,] "D"  "B"  "C"  "F"  "E"  "A" 
[6,] "B"  "G"  "D"  "C"  "A"  "F" 

round(table(as.vector(pick))/10^5,2) # aprx prob each letter chosen

   A    B    C    D    E    F    G    H    I    J    K 
0.97 0.95 0.91 0.86 0.77 0.64 0.47 0.27 0.11 0.05 0.01 
  • $\begingroup$ I'm so sorry I'm not communicating well. Basically, I'm trying to reverse engineer what the weights are. What weights would get sent to a sampling function such that it would lead to my results? $\endgroup$ – Justin Geoghegan Oct 14 '19 at 5:07
  • $\begingroup$ Fussing around a little more, I think I understand your goal. I think an exact reverse-engineering may be a combinatorial nightmare trying to account for step-wise changing of probabilities as various 'letters' get chosen. Assuming that letters get chosen roughly in alphabetical order, I made a 'reverse-cumulative' probability vector that seems to work a little better. See Addendum. With that I'll stop for now. $\endgroup$ – BruceET Oct 14 '19 at 5:50
  • $\begingroup$ Hey Bruce, what language is this, R? $\endgroup$ – Justin Geoghegan Oct 15 '19 at 5:24
  • $\begingroup$ Yes, sorry for saying only once in passing. Everything is in R. Editing now. $\endgroup$ – BruceET Oct 15 '19 at 5:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.