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.

  • $\begingroup$ Sampling with replacement, or without? When you say the "probability is triangular", the result is underdetermined - what is the ratio $p_2/p_1$? $\endgroup$ – Glen_b Jun 1 '14 at 0:18
  • $\begingroup$ 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)$). $\endgroup$ – Pat Jun 1 '14 at 1:30
  • $\begingroup$ 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. $\endgroup$ – Pat Jun 1 '14 at 1:36

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).

| cite | improve this answer | |
  • 1
    $\begingroup$ Fabulous!!! I'll thank you and that paper in my code :) $\endgroup$ – 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) );
| cite | improve this answer | |
  • $\begingroup$ 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? $\endgroup$ – Justin May 31 '14 at 21:45

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.