title: Sampling with replacement from streams
body: <p>I'm looking for efficient algorithm for sampling data with replacement from very large files (or streams). Specifically:</p>
  <ul>
  <li>There is a file large enough so that I cannot load it into memory, but can iterate through rows and write selected rows to stream;</li>
  <li>I need to sample <em>k</em> rows with replacement from it, total number of rows <em>n</em> is known in advance (<em>k</em> might be greater than <em>n</em>);</li>
  <li>Sampling has to be uniform (meaning that probability of getting any possible <em>k</em>-sized sample is same).</li>
  </ul>
  <p>I know the solution for case without replacement - keep track of <em>k_i</em> and <em>n_i</em> (number of sampled and iterated rows so far up to row <em>i</em>) and take <em>i</em>-th row with probability <em>(k - k_i) / (n - n_i)</em>; but applying same logic for sampling with replacement doesn't work - I get very complicated formulas for probability to take <em>i</em>-th row (or I'm doing it wrong).</p>
  <p>Is there an efficient way to do this?</p>
  <p><strong>UPDATED</strong></p>
  <p>Okay, seems that I found a semi-decent way to do this.</p>
  <p>Let <em>binom(n, k)</em> be binomial coefficient. Number of <em>k</em>-sized samples with replacement from <em>n</em> elements is <em>binom(n + k - 1, n - 1)</em>. Hence</p>
  <ul>
  <li>Probability to select 1st element 0 times is <em>p0 = binom(n + k - 2, n - 2) / binom(n + k - 1, n - 1) = (n - 1) / (n + k - 1)</em>;</li>
  <li>Probability to select 1st element 1 time is <em>p1 = binom(n + k - 3, n - 2) / binom(n + k - 1, n - 1) = p0 * k / (n + k - 2)</em>;</li>
  <li>Probability to select 1st element 2 times is <em>p2 = binom(n + k - 4, n - 2) / binom(n + k - 1, n - 1) = p1 * (k - 1) / (n + k - 3)</em>;</li>
  <li>etc.</li>
  </ul>
  <p>With this nice recursive formulas we can iterate through rows, roll how many times to take <em>i</em>-th row while keeping track of how many was already taken. I simulated this process with python function below:</p>
  <pre class="lang-py prettyprint-override"><code>import random
  
  
  def sample_wr(xs, k0):
      res = []
      n, k = len(xs), k0
      ix = 0
      
      while len(res) &lt; k0:
          if n == 1 and k &gt; 0:
              res.extend([xs[-1]] * k)
          else:
              m = 0
              u = random.random()
              p = (n - 1.0) / (n + k - 1.0)
              sump = p
  
              while sump &lt; u:
                  m += 1
                  p *= (k - m + 1.0) / (n + k - m - 1.0)
                  sump += p
                  
              if m:
                  res.extend([xs[ix]] * m)
                  k -= m
              
              n -= 1
              ix += 1
          
      return tuple(res)
  </code></pre>
  <p>After 10M simulations on list <em>[1, 2, 3]</em> all 4-element samples seem to come out uniformly.</p>
  <p>Thanks everyone for insights.</p>