A muon-decay experiment records decay times in the window

PHYSICS · Multi-part · Fermi University Physics · Scientific Practice & Experimental Methods · Computational Physics · Data analysis and model fitting in the laboratory

Error, Chance and the Machine, Ch Markov-chain Monte Carlo

Problem

A muon-decay experiment records N=1000N = 1000 decay times in the window 0≤t≤T=10 μs0 \le t \le T = 10\ \mu\mathrm{s}. A fraction bb of the recorded times are stray signals, uniform on the window, so one recorded time has density f(t∣τ,b)=(1−b) e−t/ττ (1−e−T/τ)+bT,f(t\mid\tau, b) = (1-b)\,\frac{e^{-t/\tau}}{\tau\,(1 - e^{-T/\tau})} + \frac{b}{T}, where τ\tau is the mean life. Generate the data with NumPy's default generator (PCG64), seed 11, τ=2.197 μs\tau = 2.197\ \mu\mathrm{s} and b=0.05b = 0.05: draw u = rng.random(1000) and then v = rng.random(1000). Event ii is a stray signal if ui<bu_i < b, and then ti=Tvit_i = Tv_i. Otherwise ti=−τln⁡[1−vi(1−e−T/τ)]t_i = -\tau\ln[1 - v_i(1 - e^{-T/\tau})]. Take the prior π(τ,b)\pi(\tau, b) flat on 0<τ≤20 μs0 < \tau \le 20\ \mu\mathrm{s} and 0≤b≤10 \le b \le 1. ((a)) Write ln⁡L(τ,b)\ln L(\tau, b) and the target p(τ,b∣t)p(\tau, b\mid\mathbf{t}). Explain why the prior needs the upper limit on τ\tau. ((b)) Sample the target by random-walk Metropolis with separate Gaussian steps sτs_\tau and sbs_b, rejecting any proposal outside the box. Tune the steps with pilot runs of 50005000 steps from (2.2 μs,0.05)(2.2\ \mu\mathrm{s}, 0.05), seed 77, until the acceptance lies between 0.30.3 and 0.50.5. Report the acceptance for the steps (0.15 μs,0.03)(0.15\ \mu\mathrm{s}, 0.03). ((c)) With those steps, run four chains from (τ0,b0)=(1.0 μs,0.30)(\tau_0, b_0) = (1.0\ \mu\mathrm{s}, 0.30), (6.0 μs,0.02)(6.0\ \mu\mathrm{s}, 0.02), (4.0 μs,0.60)(4.0\ \mu\mathrm{s}, 0.60) and (2.0 μs,0.90)(2.0\ \mu\mathrm{s}, 0.90), PCG64 seeds 100100 to 103103. Each chain discards 20002000 steps and keeps 20 00020\,000; draw xi = rng.standard_normal((n_total, 2)) and then u = rng.random(n_total). Report the acceptance, the posterior mean and standard deviation of τ\tau and of bb, the Monte Carlo error of the mean of τ\tau from 2020 batch means per chain (combine the four chains in quadrature and divide by four), the effective sample size of τ\tau, and R^\hat R for both parameters. ((d)) Generate the data again with b=0b = 0 (same seed) and fit τ\tau alone, with one chain from τ0=2.0 μs\tau_0 = 2.0\ \mu\mathrm{s}, step 0.19 μs0.19\ \mu\mathrm{s}, 10001000 steps discarded and 50 00050\,000 kept, seed 1111. Compare the posterior standard deviation with the Cramér-Rao value of Chapter 6 for τ=2.197 μs\tau = 2.197\ \mu\mathrm{s} and N=1000N = 1000. Say why the standard deviation of τ\tau in part (c) is larger. ((e)) Chapter 6 found that the mean of a windowed record without background is 2.090 μs2.090\ \mu\mathrm{s}, 4.854.85 per cent below τ\tau. Find the expected sample mean of the times in part (c), and say what it shows about agreeing with τ\tau by the naive mean.

Browse the problem library