A step's chance: the period's chance, one independent trial a step¶
The chance of an event in a step, from its chance in a period, each step an independent trial: p_dt = 1 - (1 - p_T)^(dt / T), which keeps the chance of at least once in the period, 1 - (1 - p_dt)^(T / dt) = p_T.
A chance is dimensionless; the period is the game's (Bourse's minute, 60 s), never a literal. At dt = T it returns the period's chance itself, by the branch: without it the chance would not come back exactly, since 1 - (1 - p) loses a p below 1/2 whose bits are finer than 1 - p keeps, and binary32's pow(x, 1) is x but for a few x in a thousand, one ULP off. A chance near 0 loses its low bits to 1 - p the same way (binary32 holds 1 - p to 2^-25): a step's chance is good to a few 2^-24 absolute, not relative to a small chance. Source: Ross, A First Course in Probability, 3.4 (independent trials) and 4.6 (the geometric law). Checked against a reference in numpy, which has its own tests.
- Source
- Ross, A First Course in Probability, 3.4 and 4.6.
- Proof
- curve, invariant, limit.
step_chance.trait
namespace lusine.stochastic;
float step_chance(float per_period, float<s> period, float<s> dt)
{
return dt == period ? per_period : 1.0 - pow(1.0 - per_period, dt / period);
}
The chance of an event in a step dt, from its chance in a period T, each step an independent trial, p_dt = 1 − (1 − p_T)^(dt / T), so the chance of at least once in the period, 1 − (1 − p_dt)^(T / dt), is p_T. A chance is dimensionless; the period is the game's, passed in.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
U = EPS / 2 # binary32's unit roundoff, 2^-24
f32 = lambda a: np.float32(a).astype(float)
The kernel¶
The law as the library states it.
print(trait_run.source('stochastic/step_chance.trait'))
// lusine.stochastic.step_chance: the chance of an event in a step, from its chance in a period, each step an
// independent trial: p_dt = 1 - (1 - p_T)^(dt / T), which keeps the chance of at least once in the period, 1 - (1 -
// p_dt)^(T / dt) = p_T. A chance is dimensionless; the period is the game's (Bourse's minute, 60 s), never a literal.
// At dt = T it returns the period's chance itself, by the branch: without it the chance would not come back exactly,
// since 1 - (1 - p) loses a p below 1/2 whose bits are finer than 1 - p keeps, and binary32's pow(x, 1) is x but for
// a few x in a thousand, one ULP off. A chance near 0 loses its low bits to 1 - p the same way (binary32 holds 1 - p to
// 2^-25): a step's chance is good to a few 2^-24 absolute, not relative to a small chance.
// Source: Ross, A First Course in Probability, 3.4 (independent trials) and 4.6 (the geometric law).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.stochastic;
float step_chance(float per_period, float<s> period, float<s> dt)
{
return dt == period ? per_period : 1.0 - pow(1.0 - per_period, dt / period);
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets the chance, the period and the step, and main() binds the kernel over the items.
print(trait_run.law_program('stochastic/step_chance.trait', 'step_chance'))
def chance(per_period, period, dt):
return trait_run.law('stochastic/step_chance.trait', 'step_chance', per_period=per_period, period=period, dt=dt)
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float per_period;
host float<s> period;
host float<s> dt;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.stochastic.step_chance(items.per_period, items.period, items.dt); }
Curve: 0.35 a minute is 0.35 a minute, and 0.1938 a half¶
The textbook's cases: a chance of 0.35 a minute is 0.35 at dt = 60 s, and 1 − 0.65^½ ≈ 0.1938 at dt = 30 s. Against the step the law bends below the planted fault, the chance scaled linearly (p_T dt / T), which it meets at dt = T and which passes 1 beyond.
c = chance(0.35, 60.0, [60.0, 30.0])
assert c[0] == f32(0.35) and abs(c[1] - (1 - np.sqrt(0.65))) <= 4 * U and round(c[1], 4) == 0.1938, c
print(f'dt = 60 s: {c[0]:.6f}; dt = 30 s: {c[1]:.6f}')
dt = np.linspace(0.0, 240.0, 97)
fig, ax = plt.subplots(figsize=(7, 3.2))
for p in (0.1, 0.35, 0.7):
c = chance(p, 60.0, dt)
assert np.all(np.diff(c) > 0) and np.all((c >= 0) & (c < 1)), p
line, = ax.plot(dt, c, label=f'p_T = {p}')
ax.plot(dt, p * dt / 60.0, '--', color=line.get_color(), lw=0.8)
ax.axvline(60.0, color='grey', lw=0.5)
ax.set_ylim(0, 1.05)
ax.set_xlabel('dt (s)')
ax.set_ylabel('chance in the step')
ax.legend(title='solid: the law; dashed: the planted fault')
plt.show()
dt = 60 s: 0.350000; dt = 30 s: 0.193774
Invariant: n steps compound back to the period's chance¶
Over n steps of dt = T / n, the chance of at least once is 1 − (1 − p_dt)^n, which must be p_T at any n; over an hour, 1 − (1 − p_T)^60 at any dt, as Bourse's clock asks. The page compounds the emitted binary32 chances in float64: each step's chance is good to about 2 × 2⁻²⁴ absolute (the last section), and n steps compound that to at most n times as much, so the check is held to n × 2.5 × 2⁻²⁴.
p = f32([0.02, 0.35, 0.9])
for n in (1, 2, 3, 10, 60, 600, 3600):
c = chance(p, 60.0, f32(60.0 / n))
back = 1 - (1 - c) ** (60.0 / f32(60.0 / n))
assert np.all(np.abs(back - p) <= n * 2.5 * U), (n, back - p)
hour = 1 - (1 - f32(0.02)) ** 60
for dt in f32([0.25, 1.0, 10.0, 60.0, 600.0]):
once = float(1 - (1 - chance(0.02, 60.0, dt)[0]) ** (3600.0 / dt))
assert abs(once - hour) <= (3600.0 / dt) * 2.5 * U, (dt, once, hour)
print(f'dt = {dt:7.2f} s: at least once in an hour {once:.7f}, want {hour:.7f}')
dt = 0.25 s: at least once in an hour 0.7023927, want 0.7024468 dt = 1.00 s: at least once in an hour 0.7024382, want 0.7024468 dt = 10.00 s: at least once in an hour 0.7024477, want 0.7024468 dt = 60.00 s: at least once in an hour 0.7024468, want 0.7024468 dt = 600.00 s: at least once in an hour 0.7024465, want 0.7024468
Limit: none, certain, a step of the period, and of nothing¶
A chance of nothing stays nothing and a certainty stays certain, at any step; a step of nothing has no chance (pow(x, 0) is 1). A step of the period gives the period's chance itself, bit for bit, by the kernel's branch. Without it the period's chance would not come back: the program below writes the law without the branch, and over two thousand chances it misses a good share of those below ½, since 1 − (1 − p) drops a p's bits finer than 1 − p keeps, and binary32's pow(x, 1) is not always x.
p = f32([0.0, 1.0])
for dt in (0.0, 0.5, 60.0, 240.0):
c = chance(p, 60.0, dt)
assert c[0] == 0.0 and c[1] == (0.0 if dt == 0.0 else 1.0), (dt, c)
assert np.all(chance(f32([0.01, 0.35, 0.99]), 60.0, 0.0) == 0.0)
rng = np.random.default_rng(302)
p = f32(np.concatenate([rng.uniform(0, 1, 1000), 10 ** rng.uniform(-8, 0, 1000)]))
assert np.array_equal(chance(p, 60.0, 60.0), p) and np.array_equal(chance(p, 0.5, 0.5), p)
UNBRANCHED = '''namespace probe;
lusine.version version = { { number = 1 } };
struct item { host float x; float a; float b; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);
(float, float) unbranched(float x) { return (pow(x, 1.0), 1.0 - pow(1.0 - x, 1.0)); }
void main() { (items.a, items.b) = unbranched(items.x); }
'''
run = trait_run.Program(['stochastic/step_chance.trait'], UNBRANCHED, {'probe.items': [len(p)]})
run.store('probe.items', 'x', p)
a, b = (o[0] for o in run.steps('probe.main', 1, [('probe.items', 'a'), ('probe.items', 'b')]))
assert np.any(b[p < 0.5] != p[p < 0.5]) and np.any(a != p)
print(f'with the branch: every p back exactly; without: {np.mean(b != p):.1%} missed, '
f'pow(x, 1) != x for {np.mean(a != p):.2%}, at most {np.max(np.abs(a - p) / np.spacing(np.float32(p))):.0f} ULP')
with the branch: every p back exactly; without: 62.8% missed, pow(x, 1) != x for 0.10%, at most 1 ULP
Agreement: binary32 against the float64 reference¶
Three thousand cases, chances from 10⁻⁸ to 1 and steps from a sixtieth of a second to four periods, each input a binary32: the emitted kernel against the reference in float64 on the same inputs (written as −expm1((dt / T) log1p(−p)), which keeps a small chance's digits). The kernel rounds 1 − p to binary32 first, so its error is absolute, not relative: at most 2.5 × 2⁻²⁴ (about 1.5 × 10⁻⁷) here, which is a large share of a chance of 10⁻⁷ and a small one of a chance of 0.01. Relative to the step's chance, it is held to 2.5 × 2⁻²⁴ / p_dt.
from lusine.kernels.stochastic.step_chance import step_chance as reference
rng = np.random.default_rng(302)
p = f32(np.concatenate([10 ** rng.uniform(-8, 0, 1500), rng.uniform(0, 1, 1500)]))
dt = f32(60 * 10 ** rng.uniform(-3.5, 0.6, 3000))
want = reference(p, 60.0, dt)
err = np.abs(chance(p, 60.0, dt) - want)
assert np.max(err) <= 2.5 * U, np.max(err) / U
big = want >= 0.01
rel = err[big] / want[big]
assert np.all(rel <= 2.5 * U / want[big])
print(f'absolute error at most {np.max(err):.3g} = {np.max(err) / U:.2f} x 2^-24; '
f'relative, where the step\'s chance is at least 0.01, at most {np.max(rel):.3g}')
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.loglog(want[want > 0], err[want > 0] / U, '.', ms=2)
ax.set_xlabel('the step\'s chance (float64 reference)')
ax.set_ylabel('|binary32 − float64| / 2⁻²⁴')
plt.show()
absolute error at most 1.13e-07 = 1.90 x 2^-24; relative, where the step's chance is at least 0.01, at most 3.61e-06
