Designing an experiment before there is any data · step 6 of 6
Running it
The counts come in. What does the experiment actually say?
Every station is a Poisson draw around what the apparatus would really have produced. Then the five parameters are profiled over the one the experiment is about, and the bound is read where the likelihood has climbed far enough above its minimum.
It comes back one-sided, and that is the honest shape: there is no lower limit on the radius here. Everything below the bound fits the counts equally well, because the beam cannot resolve past D/2 whatever the counting time. Rutherford published 34 fm in 1911.
the single-scattering tail, cut off by a charge of size λ
u = Apply(fn="sin", arg=Mul(factors=(
Const(value=0.5, dimension=dimensionless()), theta)))
tail = Mul(factors=(
_exp(log_a),
Pow(base=u, exponent=Fraction(-4)),
_exp(Mul(factors=(
Const(value=-1.0, dimension=dimensionless()),
Pow(base=_ratio(u, lam), exponent=Fraction(2)),
))),
))
the multiple-scattering core
core = Mul(factors=(
_exp(log_c),
_exp(Mul(factors=(
Const(value=-0.5, dimension=dimensionless()),
Pow(base=_ratio(theta, log_w), exponent=Fraction(2)),
))),
))
expected counts, and the variance-stabilising transform
scattered = Mul(factors=(
Data(name="foil", dimension=dimensionless()),
Add(terms=(tail, core)),
))
expected = Mul(factors=(
Mul(factors=(RATE_UNIT, Add(terms=(scattered, _exp(log_b))))),
Data(name="omega", dimension=dimensionless()),
Data(name="exposure", dimension=BASES.time),
))
mean = Mul(factors=(
Const(value=2.0, dimension=dimensionless()),
Pow(base=expected, exponent=Fraction(1, 2)),
))
# 2*sqrt(Poisson) has variance 1, so a normal
# likelihood on y is the right one
import warnings
from axiom.design import profile_likelihood
counts = np.random.default_rng(1911).poisson(surface.counts(plan, s.HARD_CENTRE))
measured = dict(plan)
measured["root_count"] = 2.0 * np.sqrt(counts)
inn = [bool(st.foil) for st in s.PLAN]
if_diffuse = surface.counts(plan, s.DIFFUSE)
print(f"at 150 degrees: {int(counts[inn][-1]):,} counts observed")
print(f" a diffuse atom predicts {float(if_diffuse[inn][-1]):.0f}")
print(f" that is a Poisson excess of "
f"{(counts[inn][-1] - if_diffuse[inn][-1]) / math.sqrt(counts[inn][-1]):,.0f} sd")
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
values, drops = profile_likelihood(
s.model(), measured, {**s.HARD_CENTRE, "lam": 0.0},
"lam", grid=np.linspace(-4, 4, 41),
)
drops = np.asarray(drops)
lowest = int(np.argmin(drops))
last = int(np.where(drops[:lowest] > 1.921)[0][-1])
crossing = float(np.interp(
1.921, [drops[last + 1], drops[last]], [values[last + 1], values[last]]
))
print()
print(f"R < {s.radius_of(crossing) * 1e15:.1f} fm (95%, one-sided)")
print(f"the beam could not have resolved below {s.D_CLOSEST / 2 * 1e15:.1f} fm")
print("Rutherford published 34 fm in 1911.")
at 150 degrees: 94,117 counts observed a diffuse atom predicts 119 that is a Poisson excess of 306 sd R < 33.5 fm (95%, one-sided) the beam could not have resolved below 29.6 fm Rutherford published 34 fm in 1911.
Every line above was captured by running this code at build time. The steps share one namespace, so the code on this page is the code you would type after the steps before it.