## Reservoir Sampling

### January 31, 2014

Our code is very concise:

`(define (sample k xs)`

(let ((resv (list->vector (take k xs))))

(do ((xs (drop k xs) (cdr xs)) (i k (+ i 1))) ((null? xs) resv)

(let ((j (randint i)))

(when (< j k) (vector-set! resv j (car xs)))))))

And here are some examples:

`> (define xs '(a b c d e f g h i j k l m n o p q r s t u v w x y z))`

> (sample 3 xs)

#(l z r)

> (sample 3 xs)

#(y g h)

> (sample 3 xs)

#(d k i)

We used `take`

and `drop`

from the Standard Prelude, and `rand`

and `randint`

from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/VgxH2mwk.

In Python. Note that it is important to choose a random number inclusive i. Otherwise the first k items in the sequence will be undersampled. What is interesting about this form of reservir sampling is that you do not have to know the number of items in advance. On the downside, it is slow, as you need a random sample for every item.

Fun! Python or something. (Same as Paul’s but different.)

Some test results:

100-sampling [0, 1, ..., 999] 1000 times:

population mean: 499.5

mean deviation from population mean: 0.027230000000000528

`2-sampling [3, 1, 4] 999 times:`

335 [4, 1]

328 [3, 1]

337 [3, 4]

In this paper various methods for reservoir sampling are given. The simplest is algorithm R and is the same as suggested above. Also 2 more algorithms are given, named X and Z. I give here implementations. Timings for k=10 (or 100) and N=10000000 are R: 13 sec, X: 2.4 sec and Z: 0.17 sec. (I left out the code for R, as that is the same as my last post.

There is an error in algorithm_Z of my last post. I corrected it and the full code including test code can be found at here. Timing for Z is now 0.35 sec.

In Guile Scheme.

Here’s a solution in Julia.