# Choosing elements of an array based on given probabilities for each element

I'm not well-versed in statistics at all. I'm currently working on an implementing an optimization algorithm, and am having a bit of trouble on this part.

I have an array of $n$ elements $A = \left\{a_{i}\right\}$ where $a_{i}\in\mathbb{R}$ and $1 \le i \le n$. Furthermore, I have an array of probabilities given by $P = \left\{p_{i}\right\}$ with $p_{i}\in[0,1]$. (Also, $\sum_{i} p_{i} = 1$ of course) If it matters, the probability distribution is triangular with $i=1$ having the highest probability and $i=n$ having the smallest probability.

Now I want to choose $q < n$ distinct elements from $A$ using the probabilities of selection from $P$. I'm currently coding the algorithm in MATLAB, but I'm not sure how to go about choosing the $q$ elements.

• Sampling with replacement, or without? When you say the "probability is triangular", the result is underdetermined - what is the ratio $p_2/p_1$? – Glen_b Jun 1 '14 at 0:18
• If you're sampling with replacement than the randsample function can do exactly what you ask (the exact line will be something like $samples=randsample(A,q,true,P)$). – Pat Jun 1 '14 at 1:30
• If without replacement then you'll need to specify exactly how the probabilities change when you remove a sample from the pool. The obvious way is to just renormalise the probabilites $P$ and repeat, until you've drawn $q$ samples (so you'll probably end up building a $for$ loop over $q$ iterations, each time calling randsample to draw a single weighted sample from the remaining available ones). But be sure this renormalisation makes physical sense for the system - it's entirely concievable that some systems may need a more involved re-calculation of sample probabilities. – Pat Jun 1 '14 at 1:36

## 2 Answers

To get a weighted random sample without replacement of size $m<n$ draw $n$ independent $u_i$'s with $\mathrm{Uniform}[0,1]$ distribution (using rand()), compute the keys $k_i=u_i^{1/p_i}$, and pick the $m$ elements with largest $k_i$'s. The $p_i$'s don't need to be normalized.

This amazingly simple algorithm is due to:

Efraimidis, P.S. and Spirakis, P.G. Information Processing Letters, 97, 181–185 (2006).

• Fabulous!!! I'll thank you and that paper in my code :) – Yan King Yin Jul 10 '18 at 21:03

In Matlab you can call this line q times

r(ii) = 1 + sum( rand() > cumsum(P) );

• could you explain this? I'm proficient in MATLAB but I'm not sure what this line of code is supposed to do. I used my P array but it just returns 1. And what is ii supposed to be? – Justin May 31 '14 at 21:45