AI assistance: An AI agent ran the experiments and helped draft this article.
Try this in Python:
0.5 ** 1200
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.
Imagine two fixed models of independent tosses:
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.
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.
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?
The sweep also computed each likelihood by repeatedly multiplying its previous value by p.
At n = 1,200, that route produced:
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.
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.
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)
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.
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.
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.