Python returned zero. The coin model still assigned a positive probability.
AI assistance: An AI agent ran the experiments and helped draft this article.
Try this in Python:
0.5 ** 1200
# 0.0
That is a strange result if we read zero as “this cannot happen.” Under an independent fair-coin model, 1,200 heads in a row has a positive probability. A very small one, yes. But positive.
In the cookie-jar example, some outcomes really were impossible: we could not draw more plain cookies than the jar contained. Here the formula permits the outcome. The number system runs out of room.
The interesting part is what happens when we compare two models. The calculation can lose the difference between them, or even turn a rounding artifact into apparent certainty.
Start with two almost identical coins
Imagine two fixed models of independent tosses:
- Model A assigns a heads probability of 0.5.
- Model B assigns a heads probability of 0.5001.
For a specified sequence of n heads, each model gives a likelihood of p raised to n. No combination factor is needed: the sequence is all heads.
These are synthetic arithmetic inputs. Nobody observed 1,200 heads, trained a model, or collected evidence of a real biased coin. The models were chosen to make the numerical comparison easy to inspect.
At n = 1,200, high-precision Decimal arithmetic gives approximately:
| Quantity | Model A: p = 0.5 | Model B: p = 0.5001 |
| --- | ---: | ---: |
| Likelihood | 5.807714 × 10⁻³⁶² | 7.382874 × 10⁻³⁶² |
| Natural log likelihood | −831.776616672 | −831.536640669 |
| Normalized likelihood weight | 44.02923% | 55.97077% |
Both ordinary Python powers return 0.0. Dividing either by their sum would mean 0/0, so the direct calculation cannot produce those weights.
Here, “normalized likelihood weight” means L_A / (L_A + L_B), and similarly for B. It compares these two specified models on this sequence. It is not the chance of heads on the next toss. With a separate assumption of equal prior probabilities over exactly these two models, the same numbers could be interpreted as posterior model probabilities. That extra assumption matters.
The little gap between zero and the smallest float
On the tested machine, Python floats use binary64 arithmetic. The smallest positive value prints as 5e-324. Values near zero are spaced apart; a sufficiently small positive result rounds to zero. A value just below the smallest positive float can still round up to it. The halfway boundary and rounding rule matter.
There is also a region of subnormal values below the smallest positive normal value. They extend the available range, with less relative precision. NumPy's floating-point limits documentation explains this distinction.
The local sweep checked n = 1 through 1,300. Direct powers first became subnormal at n = 1,023 for both models. Model A first rounded to zero at n = 1,075; Model B at n = 1,076. These thresholds describe the tested binary64 routes and inputs, rather than every language, data type or device.
!Saved likelihood measurements and normalized weights across synthetic sequence lengths. Direct powers disappear below binary64 range; repeated multiplication approaches zero and one while stable weights remain near 44 and 56 percent.
Agent-run arithmetic sweep, 1,300 sequence lengths. The upper chart uses log10 likelihood so the scale remains readable. The Decimal reference continues below binary64 range. The lower chart zooms into the failure region. A missing direct-power curve means normalization is undefined after both values become zero.
Okay, so perhaps we should multiply one toss at a time instead?
Multiplying slowly can give a different wrong answer
The sweep also computed each likelihood by repeatedly multiplying its previous value by p.
At n = 1,200, that route produced:
- Model A: 0.0.
- Model B: 5e-324.
Normalizing those two values gives 0% and 100%.
The second model has become stuck. Multiplying the smallest positive float by 0.5001 gives an exact intermediate value just above half of that smallest float. Rounding brings it back to the same smallest value. Repeating the operation repeats the rounding. With p = 0.5, half of the smallest float is an exact tie between zero and that smallest value; binary64 round-to-nearest, ties-to-even selects zero.
The direct power and repeated product represent the same mathematical expression, but they round at different stages. The surviving nonzero value does not establish that Model B is certain.
!At 1,200 synthetic heads, Decimal and log-space comparison give model weights of 44.03 and 55.97 percent. Repeated multiplication gives zero and 100 percent. Direct-power normalization is undefined.
Comparison on the same saved n = 1,200 inputs. The apparent 100% weight comes from rounding in repeated multiplication. Both direct powers are zero, so there is no valid direct-power bar to display.
This is a useful reason to inspect a calculation even when its final weights add up to one.
Keep the logarithm before the probability disappears
Products become sums in log space:
log L = log(p × p × … × p) = n × log(p)
The log likelihoods around −832 fit comfortably in binary64, even though the corresponding probabilities do not.
For comparison, we also need the logarithm of the sum of the likelihoods. Computing the tiny likelihoods first would lose them again. The standard log-sum-exp identity handles this:
m = max(a, b)
log(exp(a) + exp(b)) = m + log(exp(a − m) + exp(b − m))
For finite a and b, the shifted exponents are nonpositive and at least one is zero. The shared scale is removed before exponentiation. SciPy's logsumexp implements a stable calculation, including cases beyond this two-number example.
Here is the tested calculation:
import math
import numpy as np
from scipy.special import logsumexp
n = 1200
log_likelihoods = np.array([
n * math.log(0.5),
n * math.log(0.5001),
])
weights = np.exp(log_likelihoods - logsumexp(log_likelihoods))
print(weights)
# [0.44029226 0.55970774]
The log-likelihood difference is approximately 0.2399760032. For just two models, we can check Model B's weight directly:
weight_B = 1 / (1 + exp(−0.2399760032)) ≈ 0.5597077356
This is the same comparison after the common scale cancels. The calculation retains useful relative information with small rounding error; the logs are not exact.
The weights are ordinary-sized numbers because the subtraction removes the shared tiny scale. Exponentiating the original log likelihoods still returns zero on this machine. Taking a logarithm after a value has already underflowed cannot recover what was lost.
Check both extremes and a real zero
The checks also included inputs outside the coin example.
For log values −1,000 and −1,001, the naive exponential sum underflowed and its logarithm became minus infinity. SciPy returned approximately −999.686738312.
For 1,000 and 1,001, the naive exponentials overflowed. SciPy returned approximately 1,001.313261688.
Adding a common shift of −1,000 or +1,000 to a separate two-score example preserved the normalized weights within about 3.45 × 10⁻¹⁵ in this run. That is the shared-scale cancellation in action.
A genuinely impossible event still has probability zero and log probability minus infinity. Two minus-infinity inputs give a log sum of minus infinity. They do not define useful normalized model weights. A stable function cannot supply evidence that the model does not contain.
There is another boundary before multiplication: rounding the input itself. With p = 10⁻²⁰, ordinary 1 - p rounded to 1.0, and math.log(1 - p) returned 0.0. math.log1p(-p) retained approximately −10⁻²⁰. Python documents log1p for accurate logarithms of one plus a small input.
What was actually checked
The local run used Python 3.13.9, NumPy 2.5.4 and SciPy 1.18.1 on macOS arm64. The inputs and tolerances were fixed before execution.
For each of the 1,300 sequence lengths, Decimal references used exact decimal inputs at both 100 and 160 digits of precision. The two reference calculations agreed to better than 10⁻⁹⁰ under the recorded relative-probability and absolute-weight checks. Small cases at n = 1, 2 and 10 also matched exact fractions.
Across the grid, the largest absolute log-likelihood error was about 1.27 × 10⁻¹³. The largest absolute normalized-weight error was about 7.01 × 10⁻¹⁴. Both were below the preset tolerances of 10⁻¹⁰ and 10⁻¹². Errors were compared against high-precision Decimal values before reducing the reported discrepancy to a float.
The comparison includes binary64 input rounding: 0.5001 has an exact decimal reference but is not exactly representable as a binary float. This small input difference is part of the measured discrepancy.
The installed SciPy source matched its 1.18.1 tagged source. Reading the 1.14.1 and 1.16.2 source also showed maximum-shift calculations; those older versions were not executed. Numerical stability here is an established technique, not a new library improvement. For a deeper treatment, see Blanchard, D. J. Higham and N. J. Higham's analysis of log-sum-exp and softmax.
There is also a model-choice catch. The weights compare only the candidates supplied. A tidy 56/44 split does not establish that either model describes the data well. Tiny likelihoods of exact long sequences alone do not establish that a model is wrong, because many possible sequences each have tiny probability. An all-heads count of 1,200 would be striking under either nearly fair model; checking that statistic against the model is a separate question from normalizing two likelihoods.
The experiment covers two small synthetic models and a few numerical controls. It does not establish accuracy for every distribution, backend, precision, cancellation pattern or trained model. No SciPy defect was identified.
A displayed zero can therefore need two very different explanations: the model rules out the event, or the chosen arithmetic has rounded away a positive value. Keeping likelihoods in log space preserved the comparison in this example. Before interpreting a zero, ask whether it comes from what the model permits or from how the computer stores the result.
10 reactions · 2 comments
Probabilities can sum to one and still describe the wrong experiment
AI assistance: An AI agent ran the experiments and helped draft this article.
Imagine a jar with seven chocolate cookies and three plain cookies. You are
packing a box of six. Each time, you pick blindly, so every cookie still in the
jar has the same chance of being selected. Once a cookie goes into the box, it
stays there.
Before working out any probabilities, try to picture a box containing only two
chocolate cookies. The other four would have to be plain.
But the jar only had three plain cookies to begin with.
Oh, that is the catch. This is already impossible, before a probability formula
has entered the picture.
Start with what can actually happen
Use all three plain cookies and there are still three places left in the box.
Those places must contain chocolate cookies. So the smallest possible chocolate
count is six minus three, or three.
The largest possible count is six. There are seven chocolate cookies available,
but the box has only six places.
That leaves just four possible counts: 3, 4, 5 and 6.
!An illustrated cookie jar has seven chocolate and three plain cookies. A six-cookie box with two chocolate would need an unavailable fourth plain cookie. Valid boxes contain at least three and at most six chocolate cookies.
*An original illustrated sequence of possible and impossible boxes—not a
recorded random draw. C means chocolate; P means plain. The final scene shows
one of the possible four-chocolate boxes. [Static storyboard of the starting
jar and the two boundaries](panditshivam.github.io).*
This is a counting constraint, or a pigeonhole-style argument: there are not
enough plain cookies to fill more than three of the six places. The constraint
comes from what is available, not from how lucky the draw is. Wikipedia's
pigeonhole principle is a
useful introduction to this kind of reasoning.
For a finite count like this, the set of values with nonzero probability is
called the support. Here that set is {3, 4, 5, 6}. The
hypergeometric distribution
describes this sort of count: selecting uniformly from a finite population
without putting the selected objects back.
Okay, now we know what can happen. That still does not tell us which of those
four counts is most likely.
Four possible counts does not mean four equal chances
It is tempting to give 3, 4, 5 and 6 a probability of one quarter each. But those
counts can be made in different numbers of ways.
Treat the ten cookies as ten distinct objects. There are 210 different sets of
six cookies that could go into the box. Under the drawing rule above, those
sets are equally likely.
Here is what exact enumeration of all 210 sets gives:
| Chocolate cookies | Boxes out of 210 | Exact probability | Percentage |
| --- | ---: | ---: | ---: |
| 2 | 0 | 0 | 0% |
| 3 | 35 | 1/6 | 16.6667% |
| 4 | 105 | 1/2 | 50% |
| 5 | 63 | 3/10 | 30% |
| 6 | 7 | 1/30 | 3.3333% |
| 7 | 0 | 0 | 0% |
So a box with four chocolate cookies has probability 50%, while six
chocolate cookies has probability about 3.33%. Both are possible; they are
far from equally likely.
The difference becomes easier to see by counting a four-chocolate box. We need
to choose four of the seven chocolate cookies and two of the three plain ones.
C(7, 4) × C(3, 2) = 35 × 3 = 105 boxes
Here C(a, b) means the number of ways to choose b distinct objects from a,
ignoring their order. The total number of six-cookie sets is C(10, 6) = 210.
Dividing gives:
P(four chocolate) = 105 / 210 = 1/2
That is the idea behind the hypergeometric formula. It counts the selections
that fit the requested chocolate count, then divides by all possible selections.
SciPy's hypergeometric documentation
gives the general formula and the bounds used by its implementation.
Putting the cookie back changes the question
Now imagine a slightly different activity. Pick a cookie, record whether it is
chocolate, put it back, and mix the jar before the next pick. Do that six
times. This counts chocolate results in six draws; it does not pack six distinct
cookies into a box. The same physical cookie can appear again.
Every draw now has a chocolate probability of 7/10 = 0.7. Returning the cookie
keeps the jar's composition unchanged, and independent uniform picks give six
independent trials with that fixed probability. That is a
binomial model, written as
Binomial(6, 0.7).
In the original box-packing activity, the next probability depends on what was
already taken:
- After taking a chocolate cookie, 6 of the remaining 9 are chocolate.
- After taking a plain cookie, 7 of the remaining 9 are chocolate.
So the binomial formula is answering a valid question. It is just a different
question from packing six distinct cookies without replacement.
!Two probability bar charts use the same axes. Without replacement, only chocolate counts 3–6 have positive probability. With replacement, counts 0–2 also receive probability, totaling 7.047%.
*Top: exact counts from all 210 distinct six-cookie boxes. Bottom: exact
Binomial(6, 0.7) arithmetic for six independent draws with replacement. The
shaded counts are impossible in the original packing activity. Both charts use
the saved experiment values; the cookie story changes the explanation, not the
data.*
Adding up to one is only one check
The binomial model assigns probabilities of 0.0729%, 1.0206% and 5.9535% to
zero, one and two chocolate results. Together, those add to 7.047%.
Those results can happen when cookies are returned to the jar. They cannot
happen in the original six-cookie box, because there are only three plain
cookies available.
And yet both distributions add up to one. Both also give an expected chocolate
count of 4.2. That is an average over the possible results, not a claim that
a single box contains 4.2 intact cookies.
This is why checking the total and the mean is not enough. A distribution can
pass both checks while representing the wrong process.
The weights on possible outcomes matter too. For four chocolate results, the
binomial model gives 32.4135%; the original without-replacement calculation
gives 50%. Simply knowing that four is possible does not tell us which weight
belongs there.
There is one small distinction worth keeping straight. Probability assigned to
impossible outcomes means the possible outcomes collectively receive too little
mass. It does not, in general, prove that every individual possible outcome has
the wrong weight. In this particular comparison, all four feasible weights do
differ.
Connecting the jar back to the symbols
Once the example makes sense, the general bounds are less intimidating. Using
SciPy's names, let M be the total number of cookies, n the chocolate count in
the jar, and N the number selected. For our jar, M = 10, n = 7, N = 6.
The chocolate count k must be an integer with:
max(0, N - (M - n)) <= k <= min(n, N)
On the left, M - n is the plain-cookie supply. Any places left after using that
supply must be chocolate. The max keeps the lower bound from going below zero.
On the right, the chocolate count cannot exceed either the chocolate supply or
the number of places in the box.
One easy source of confusion: letters vary across references. Wikipedia writes
this same setup with population size N = 10, chocolate supply K = 7, and draw
size n = 6. Match the meanings before substituting the numbers.
What the small experiment checked
This is an illustration of established probability, not a new result. The local
experiment enumerated the 210 sets directly, compared them with a separate
combination calculation, and checked SciPy 1.18.1. The exact fractions sum to one.
It also covered every valid integer marked supply and draw size for populations
of 1–16 objects: 1,784 parameter combinations and 16,164 point-probability
comparisons. The largest difference from exact enumeration was about
3.33e-16, within the predefined absolute tolerance of 1e-12.
Five boundary checks included counts 2, 3, 6, 7 and 3.5. SciPy returned zero for
impossible counts under valid parameters, including the noninteger count 3.5.
Trying to draw eleven cookies from ten is a different problem: the parameters
are invalid, and that check returned NaN.
A first-run check also required exact floating-point equality to one for an
empty draw. SciPy returned 0.9999999999999997. Using the same predefined
tolerance as the other comparisons fixed that check; the exact answer is one.
These are small synthetic cases, not a large-parameter stability study. Six out
of ten is deliberately a large sampling fraction that makes the distinction
visible. This example does not say that a binomial approximation is always
unsuitable; an approximation needs to be judged against the errors its intended
use can tolerate.
The habit to keep is simple: picture how the selection happens, write down what
it permits, and then choose a probability model. The formula should follow the
story.
References for going a little deeper
- Wikipedia: Pigeonhole principle — an introduction to counting constraints.
- Wikipedia: Hypergeometric distribution — finite-population sampling without replacement, its formula and examples.
- Wikipedia: Binomial distribution — independent trials with a fixed probability.
- SciPy: hypergeom and binom — parameter definitions and the implemented probability functions.
- SciPy 1.18.1 source — the version inspected for the numerical comparison.
Follow-up: The cookie jar behind A/B testing: Fisher’s test to churn connects this example to randomized tests, Simpson’s paradox, churn and shared-stock interference.
Related experiment: When a positive likelihood becomes zero follows the numerical side: underflow, a misleading 100% weight, and stable log-space comparison.
9 reactions
Mentions
Irodori-TTS on macOS 27.0.1でIrodori-TTSのモデル読み込みに失敗するの、SciPyの読み込みに失敗してるからっぽい scipy#25635
Python returned zero. The coin model still assigned a positive probability. AI assistance: An AI agent ran the experiments and helped draft this article. Try this in Python: 0.5 ** 1200 # 0.0 That is a strange result if we read zero as “this cannot happen.” Under an independent fair-coin model, 1,200 heads in a row has a positive probability. A very small one, yes. But positive. In the cookie-jar example, some outcomes really were impossible: we could not draw more plain cookies than the jar contained. Here the formula permits the outcome. The number system runs out of room. The interesting part is what happens when we compare two models. The calculation can lose the difference between them, or even turn a rounding artifact into apparent certainty. Start with two almost identical coins Imagine two fixed models of independent tosses: - Model A assigns a heads probability of 0.5. - Model B assigns a heads probability of 0.5001. For a specified sequence of n heads, each model gives a likelihood of p raised to n. No combination factor is needed: the sequence is all heads. These are synthetic arithmetic inputs. Nobody observed 1,200 heads, trained a model, or collected evidence of a real biased coin. The models were chosen to make the numerical comparison easy to inspect. At n = 1,200, high-precision Decimal arithmetic gives approximately: | Quantity | Model A: p = 0.5 | Model B: p = 0.5001 | | --- | ---: | ---: | | Likelihood | 5.807714 × 10⁻³⁶² | 7.382874 × 10⁻³⁶² | | Natural log likelihood | −831.776616672 | −831.536640669 | | Normalized likelihood weight | 44.02923% | 55.97077% | Both ordinary Python powers return 0.0. Dividing either by their sum would mean 0/0, so the direct calculation cannot produce those weights. Here, “normalized likelihood weight” means L_A / (L_A + L_B), and similarly for B. It compares these two specified models on this sequence. It is not the chance of heads on the next toss. With a separate assumption of equal prior probabilities over exactly these two models, the same numbers could be interpreted as posterior model probabilities. That extra assumption matters. The little gap between zero and the smallest float On the tested machine, Python floats use binary64 arithmetic. The smallest positive value prints as 5e-324. Values near zero are spaced apart; a sufficiently small positive result rounds to zero. A value just below the smallest positive float can still round up to it. The halfway boundary and rounding rule matter. There is also a region of subnormal values below the smallest positive normal value. They extend the available range, with less relative precision. NumPy's floating-point limits documentation explains this distinction. The local sweep checked n = 1 through 1,300. Direct powers first became subnormal at n = 1,023 for both models. Model A first rounded to zero at n = 1,075; Model B at n = 1,076. These thresholds describe the tested binary64 routes and inputs, rather than every language, data type or device. !Saved likelihood measurements and normalized weights across synthetic sequence lengths. Direct powers disappear below binary64 range; repeated multiplication approaches zero and one while stable weights remain near 44 and 56 percent. Agent-run arithmetic sweep, 1,300 sequence lengths. The upper chart uses log10 likelihood so the scale remains readable. The Decimal reference continues below binary64 range. The lower chart zooms into the failure region. A missing direct-power curve means normalization is undefined after both values become zero. Okay, so perhaps we should multiply one toss at a time instead? Multiplying slowly can give a different wrong answer The sweep also computed each likelihood by repeatedly multiplying its previous value by p. At n = 1,200, that route produced: - Model A: 0.0. - Model B: 5e-324. Normalizing those two values gives 0% and 100%. The second model has become stuck. Multiplying the smallest positive float by 0.5001 gives an exact intermediate value just above half of that smallest float. Rounding brings it back to the same smallest value. Repeating the operation repeats the rounding. With p = 0.5, half of the smallest float is an exact tie between zero and that smallest value; binary64 round-to-nearest, ties-to-even selects zero. The direct power and repeated product represent the same mathematical expression, but they round at different stages. The surviving nonzero value does not establish that Model B is certain. !At 1,200 synthetic heads, Decimal and log-space comparison give model weights of 44.03 and 55.97 percent. Repeated multiplication gives zero and 100 percent. Direct-power normalization is undefined. Comparison on the same saved n = 1,200 inputs. The apparent 100% weight comes from rounding in repeated multiplication. Both direct powers are zero, so there is no valid direct-power bar to display. This is a useful reason to inspect a calculation even when its final weights add up to one. Keep the logarithm before the probability disappears Products become sums in log space: log L = log(p × p × … × p) = n × log(p) The log likelihoods around −832 fit comfortably in binary64, even though the corresponding probabilities do not. For comparison, we also need the logarithm of the sum of the likelihoods. Computing the tiny likelihoods first would lose them again. The standard log-sum-exp identity handles this: m = max(a, b) log(exp(a) + exp(b)) = m + log(exp(a − m) + exp(b − m)) For finite a and b, the shifted exponents are nonpositive and at least one is zero. The shared scale is removed before exponentiation. SciPy's logsumexp implements a stable calculation, including cases beyond this two-number example. Here is the tested calculation: import math import numpy as np from scipy.special import logsumexp n = 1200 log_likelihoods = np.array([ n * math.log(0.5), n * math.log(0.5001), ]) weights = np.exp(log_likelihoods - logsumexp(log_likelihoods)) print(weights) # [0.44029226 0.55970774] The log-likelihood difference is approximately 0.2399760032. For just two models, we can check Model B's weight directly: weight_B = 1 / (1 + exp(−0.2399760032)) ≈ 0.5597077356 This is the same comparison after the common scale cancels. The calculation retains useful relative information with small rounding error; the logs are not exact. The weights are ordinary-sized numbers because the subtraction removes the shared tiny scale. Exponentiating the original log likelihoods still returns zero on this machine. Taking a logarithm after a value has already underflowed cannot recover what was lost. Check both extremes and a real zero The checks also included inputs outside the coin example. For log values −1,000 and −1,001, the naive exponential sum underflowed and its logarithm became minus infinity. SciPy returned approximately −999.686738312. For 1,000 and 1,001, the naive exponentials overflowed. SciPy returned approximately 1,001.313261688. Adding a common shift of −1,000 or +1,000 to a separate two-score example preserved the normalized weights within about 3.45 × 10⁻¹⁵ in this run. That is the shared-scale cancellation in action. A genuinely impossible event still has probability zero and log probability minus infinity. Two minus-infinity inputs give a log sum of minus infinity. They do not define useful normalized model weights. A stable function cannot supply evidence that the model does not contain. There is another boundary before multiplication: rounding the input itself. With p = 10⁻²⁰, ordinary 1 - p rounded to 1.0, and math.log(1 - p) returned 0.0. math.log1p(-p) retained approximately −10⁻²⁰. Python documents log1p for accurate logarithms of one plus a small input. What was actually checked The local run used Python 3.13.9, NumPy 2.5.4 and SciPy 1.18.1 on macOS arm64. The inputs and tolerances were fixed before execution. For each of the 1,300 sequence lengths, Decimal references used exact decimal inputs at both 100 and 160 digits of precision. The two reference calculations agreed to better than 10⁻⁹⁰ under the recorded relative-probability and absolute-weight checks. Small cases at n = 1, 2 and 10 also matched exact fractions. Across the grid, the largest absolute log-likelihood error was about 1.27 × 10⁻¹³. The largest absolute normalized-weight error was about 7.01 × 10⁻¹⁴. Both were below the preset tolerances of 10⁻¹⁰ and 10⁻¹². Errors were compared against high-precision Decimal values before reducing the reported discrepancy to a float. The comparison includes binary64 input rounding: 0.5001 has an exact decimal reference but is not exactly representable as a binary float. This small input difference is part of the measured discrepancy. The installed SciPy source matched its 1.18.1 tagged source. Reading the 1.14.1 and 1.16.2 source also showed maximum-shift calculations; those older versions were not executed. Numerical stability here is an established technique, not a new library improvement. For a deeper treatment, see Blanchard, D. J. Higham and N. J. Higham's analysis of log-sum-exp and softmax. There is also a model-choice catch. The weights compare only the candidates supplied. A tidy 56/44 split does not establish that either model describes the data well. Tiny likelihoods of exact long sequences alone do not establish that a model is wrong, because many possible sequences each have tiny probability. An all-heads count of 1,200 would be striking under either nearly fair model; checking that statistic against the model is a separate question from normalizing two likelihoods. The experiment covers two small synthetic models and a few numerical controls. It does not establish accuracy for every distribution, backend, precision, cancellation pattern or trained model. No SciPy defect was identified. A displayed zero can therefore need two very different explanations: the model rules out the event, or the chosen arithmetic has rounded away a positive value. Keeping likelihoods in log space preserved the comparison in this example. Before interpreting a zero, ask whether it comes from what the model permits or from how the computer stores the result.
Probabilities can sum to one and still describe the wrong experiment AI assistance: An AI agent ran the experiments and helped draft this article. Imagine a jar with seven chocolate cookies and three plain cookies. You are packing a box of six. Each time, you pick blindly, so every cookie still in the jar has the same chance of being selected. Once a cookie goes into the box, it stays there. Before working out any probabilities, try to picture a box containing only two chocolate cookies. The other four would have to be plain. But the jar only had three plain cookies to begin with. Oh, that is the catch. This is already impossible, before a probability formula has entered the picture. Start with what can actually happen Use all three plain cookies and there are still three places left in the box. Those places must contain chocolate cookies. So the smallest possible chocolate count is six minus three, or three. The largest possible count is six. There are seven chocolate cookies available, but the box has only six places. That leaves just four possible counts: 3, 4, 5 and 6. !An illustrated cookie jar has seven chocolate and three plain cookies. A six-cookie box with two chocolate would need an unavailable fourth plain cookie. Valid boxes contain at least three and at most six chocolate cookies. *An original illustrated sequence of possible and impossible boxes—not a recorded random draw. C means chocolate; P means plain. The final scene shows one of the possible four-chocolate boxes. [Static storyboard of the starting jar and the two boundaries](panditshivam.github.io).* This is a counting constraint, or a pigeonhole-style argument: there are not enough plain cookies to fill more than three of the six places. The constraint comes from what is available, not from how lucky the draw is. Wikipedia's pigeonhole principle is a useful introduction to this kind of reasoning. For a finite count like this, the set of values with nonzero probability is called the support. Here that set is {3, 4, 5, 6}. The hypergeometric distribution describes this sort of count: selecting uniformly from a finite population without putting the selected objects back. Okay, now we know what can happen. That still does not tell us which of those four counts is most likely. Four possible counts does not mean four equal chances It is tempting to give 3, 4, 5 and 6 a probability of one quarter each. But those counts can be made in different numbers of ways. Treat the ten cookies as ten distinct objects. There are 210 different sets of six cookies that could go into the box. Under the drawing rule above, those sets are equally likely. Here is what exact enumeration of all 210 sets gives: | Chocolate cookies | Boxes out of 210 | Exact probability | Percentage | | --- | ---: | ---: | ---: | | 2 | 0 | 0 | 0% | | 3 | 35 | 1/6 | 16.6667% | | 4 | 105 | 1/2 | 50% | | 5 | 63 | 3/10 | 30% | | 6 | 7 | 1/30 | 3.3333% | | 7 | 0 | 0 | 0% | So a box with four chocolate cookies has probability 50%, while six chocolate cookies has probability about 3.33%. Both are possible; they are far from equally likely. The difference becomes easier to see by counting a four-chocolate box. We need to choose four of the seven chocolate cookies and two of the three plain ones. C(7, 4) × C(3, 2) = 35 × 3 = 105 boxes Here C(a, b) means the number of ways to choose b distinct objects from a, ignoring their order. The total number of six-cookie sets is C(10, 6) = 210. Dividing gives: P(four chocolate) = 105 / 210 = 1/2 That is the idea behind the hypergeometric formula. It counts the selections that fit the requested chocolate count, then divides by all possible selections. SciPy's hypergeometric documentation gives the general formula and the bounds used by its implementation. Putting the cookie back changes the question Now imagine a slightly different activity. Pick a cookie, record whether it is chocolate, put it back, and mix the jar before the next pick. Do that six times. This counts chocolate results in six draws; it does not pack six distinct cookies into a box. The same physical cookie can appear again. Every draw now has a chocolate probability of 7/10 = 0.7. Returning the cookie keeps the jar's composition unchanged, and independent uniform picks give six independent trials with that fixed probability. That is a binomial model, written as Binomial(6, 0.7). In the original box-packing activity, the next probability depends on what was already taken: - After taking a chocolate cookie, 6 of the remaining 9 are chocolate. - After taking a plain cookie, 7 of the remaining 9 are chocolate. So the binomial formula is answering a valid question. It is just a different question from packing six distinct cookies without replacement. !Two probability bar charts use the same axes. Without replacement, only chocolate counts 3–6 have positive probability. With replacement, counts 0–2 also receive probability, totaling 7.047%. *Top: exact counts from all 210 distinct six-cookie boxes. Bottom: exact Binomial(6, 0.7) arithmetic for six independent draws with replacement. The shaded counts are impossible in the original packing activity. Both charts use the saved experiment values; the cookie story changes the explanation, not the data.* Adding up to one is only one check The binomial model assigns probabilities of 0.0729%, 1.0206% and 5.9535% to zero, one and two chocolate results. Together, those add to 7.047%. Those results can happen when cookies are returned to the jar. They cannot happen in the original six-cookie box, because there are only three plain cookies available. And yet both distributions add up to one. Both also give an expected chocolate count of 4.2. That is an average over the possible results, not a claim that a single box contains 4.2 intact cookies. This is why checking the total and the mean is not enough. A distribution can pass both checks while representing the wrong process. The weights on possible outcomes matter too. For four chocolate results, the binomial model gives 32.4135%; the original without-replacement calculation gives 50%. Simply knowing that four is possible does not tell us which weight belongs there. There is one small distinction worth keeping straight. Probability assigned to impossible outcomes means the possible outcomes collectively receive too little mass. It does not, in general, prove that every individual possible outcome has the wrong weight. In this particular comparison, all four feasible weights do differ. Connecting the jar back to the symbols Once the example makes sense, the general bounds are less intimidating. Using SciPy's names, let M be the total number of cookies, n the chocolate count in the jar, and N the number selected. For our jar, M = 10, n = 7, N = 6. The chocolate count k must be an integer with: max(0, N - (M - n)) <= k <= min(n, N) On the left, M - n is the plain-cookie supply. Any places left after using that supply must be chocolate. The max keeps the lower bound from going below zero. On the right, the chocolate count cannot exceed either the chocolate supply or the number of places in the box. One easy source of confusion: letters vary across references. Wikipedia writes this same setup with population size N = 10, chocolate supply K = 7, and draw size n = 6. Match the meanings before substituting the numbers. What the small experiment checked This is an illustration of established probability, not a new result. The local experiment enumerated the 210 sets directly, compared them with a separate combination calculation, and checked SciPy 1.18.1. The exact fractions sum to one. It also covered every valid integer marked supply and draw size for populations of 1–16 objects: 1,784 parameter combinations and 16,164 point-probability comparisons. The largest difference from exact enumeration was about 3.33e-16, within the predefined absolute tolerance of 1e-12. Five boundary checks included counts 2, 3, 6, 7 and 3.5. SciPy returned zero for impossible counts under valid parameters, including the noninteger count 3.5. Trying to draw eleven cookies from ten is a different problem: the parameters are invalid, and that check returned NaN. A first-run check also required exact floating-point equality to one for an empty draw. SciPy returned 0.9999999999999997. Using the same predefined tolerance as the other comparisons fixed that check; the exact answer is one. These are small synthetic cases, not a large-parameter stability study. Six out of ten is deliberately a large sampling fraction that makes the distinction visible. This example does not say that a binomial approximation is always unsuitable; an approximation needs to be judged against the errors its intended use can tolerate. The habit to keep is simple: picture how the selection happens, write down what it permits, and then choose a probability model. The formula should follow the story. References for going a little deeper - Wikipedia: Pigeonhole principle — an introduction to counting constraints. - Wikipedia: Hypergeometric distribution — finite-population sampling without replacement, its formula and examples. - Wikipedia: Binomial distribution — independent trials with a fixed probability. - SciPy: hypergeom and binom — parameter definitions and the implemented probability functions. - SciPy 1.18.1 source — the version inspected for the numerical comparison. Follow-up: The cookie jar behind A/B testing: Fisher’s test to churn connects this example to randomized tests, Simpson’s paradox, churn and shared-stock interference. Related experiment: When a positive likelihood becomes zero follows the numerical side: underflow, a misleading 100% weight, and stable log-space comparison.
Second one today 🫡 scipy#26092