<p>Since NumPy 1.17, the reason is largely backward compatibility. See also <a href="https://stackoverflow.com/questions/40914862/why-is-random-sample-faster-than-numpys-random-choice/62951059#62951059">this question</a> and <a href="https://stackoverflow.com/questions/62309424/does-numpy-random-seed-always-give-the-same-random-number-every-time/62310046#62310046">this question</a>.</p>
<p>As of NumPy 1.17, <code>numpy.random.*</code> functions, including <code>numpy.random.choice</code>, are legacy functions and &quot;SHALL remain the same as they currently are&quot;, according to <a href="https://github.com/numpy/numpy/blob/master/doc/neps/nep-0019-rng-policy.rst" rel="nofollow noreferrer">NumPy's new RNG policy</a>, which also introduced a <a href="https://numpy.org/doc/stable/reference/random/index.html" rel="nofollow noreferrer">new random generation system for NumPy</a>. The reasons for making them legacy functions include the recommendation to avoid global state. Even so, however, NumPy did not deprecate any <code>numpy.random.*</code> functions in version 1.17, although a future version of NumPy might.</p>
<p>Recall that in your examples, <code>numpy.random.choice</code> takes an array of <code>float</code>s as weights. An array of integer weights would lead to more exact random number generation. And although any <code>float</code> could be converted to a rational number (leading to rational-valued weights and thus integer weights), the legacy NumPy version appears not to do this. These and other implementation decisions in <code>numpy.random.choice</code> can't be changed without breaking backward compatibility.</p>
<p>By the way, arithmetic coding is not the only algorithm that seeks to avoid wasting bits. Perhaps the canonical algorithm for sampling for a discrete distribution is the Knuth and Yao algorithm (1976), which exactly chooses a random integer based on the binary expansion of the probabilities involved, and treats the problem as a random walk on a binary tree. (This algorithm uses, on average, up to 2 bits away from the theoretical lower bound.) Any other integer generating algorithm can be ultimately described in the same way, namely as a random walk on a binary tree. For example, the <a href="https://arxiv.org/abs/2003.03830v2" rel="nofollow noreferrer">Fast Loaded Dice Roller</a> is a recent algorithm that has a guaranteed bound on the average number of bits it uses (in this case, no more than 6 bits away from the theoretical lower bound). The Han and Hoshi algorithm (from 1997) is another of this kind, but uses cumulative probabilities. See also my section, &quot;<a href="https://peteroupc.github.io/randomfunc.html#Weighted_Choice_With_Replacement" rel="nofollow noreferrer">Weighted Choice With Replacement</a>&quot;.</p>
Why does numpy.random.choice not use arithmetic coding?
<p>If I evaluate something like:</p>
<p><code>numpy.random.choice(2, size=100000, p=[0.01, 0.99])</code></p>
<p>using one uniformly-distributed random <code>float</code>, say <code>r</code>, and deciding if <code>r &lt; 0.01</code> will presumably waste many of the random bits (entropy) generated. I've heard (second-hand) that generating psuedo-random numbers is computationally expensive, so I assumed that <code>numpy</code> would not be doing that, and rather would use a scheme like <a href="https://en.wikipedia.org/wiki/Arithmetic_coding" rel="nofollow noreferrer">arithmetic coding</a> in this case.</p>
<p>However, at first <a href="https://github.com/numpy/numpy/blob/master/numpy/random/mtrand.pyx#L949" rel="nofollow noreferrer">glance</a> it appears that <code>choice</code> does indeed generate a <code>float</code> for every sample it is asked for. Further, a quick <code>timeit</code> experiment shows that generating <code>n</code> uniform floats is actually quicker than <code>n</code> samples from <code>p=[0.01, 0.99]</code>.</p>
<pre><code>&gt;&gt;&gt; timeit.timeit(lambda : numpy.random.choice(2, size=100000, p=[0.01, 0.99]), number=1000)
1.74494537999999
&gt;&gt;&gt; timeit.timeit(lambda : numpy.random.random(size=100000), number=1000)
0.8165735180009506
</code></pre>
<p>Does <code>choice</code> really generate a <code>float</code> for each sample, as it would appear? Would it not significantly improve performance to use some compression algorithm in some cases (particularly if <code>size</code> is large and <code>p</code> is distributed unevenly)? If not, why not?</p>