marimo-team/marimo-learn

1# /// script2# requires-python = ">=3.13"3# dependencies = [4#     "altair",5#     "asimpy",6#     "marimo",7#     "polars==1.24.0",8# ]9# ///10 11import marimo12 13__generated_with = "0.20.4"14app = marimo.App(width="medium")15 16 17@app.cell(hide_code=True)18def _():19    import marimo as mo20    import random21    import statistics22 23    import altair as alt24    import polars as pl25 26    from asimpy import Environment, Process27 28    return Environment, Process, alt, mo, pl, random, statistics29 30 31@app.cell(hide_code=True)32def _(mo):33    mo.md(r"""34    # The Inspector's Paradox35 36    ## *Why the Bus Is Always Late*37 38    Buses arrive at a stop with some average headway (gap between buses) of $\mu$ minutes. A passenger arrives at a uniformly random time and waits for the next bus. How long do they wait? The naive answer is $\mu / 2$: on average you land in the middle of a gap. The correct answer is almost always longer—sometimes much longer.39 40    The expected wait is not $\mu/2$ but:41 42    $$E[\text{wait}] = \frac{\mu}{2} + \frac{\sigma^2}{2\mu}$$43 44    where $\sigma^2 = \text{Var}[\text{headway}]$. The second term is always non-negative, so higher variance always means longer expected waits, even when the mean headway is unchanged.45 46    ### Three Bus Schedules with Mean Headway $\mu = 10$47 48    | Schedule    | $\sigma^2$ | Predicted wait | Naive wait |49    |-------------|-----------|----------------|-----------|50    | Regular     | 0         | 5.0            | 5.0       |51    | Exponential | 100       | 10.0           | 5.0       |52    | Clustered   | 64        | 8.2            | 5.0       |53 54    For exponentially distributed headways, $\sigma^2 = \mu^2$, so:55 56    $$E[\text{wait}] = \frac{\mu}{2} + \frac{\mu^2}{2\mu} = \mu$$57 58    A passenger waits on average for an *entire* mean headway — twice the naive expectation.59 60    ## Why This Happens: Length-Biased Sampling61 62    A passenger arriving at a random time is more likely to land inside a *long* gap than a short one, because long gaps occupy more time on the clock. This is called *length-biased sampling*. The interval containing your arrival is not a random headway: it is drawn from the length-biased distribution with density:63 64    $$f^*(h) = \frac{h \cdot f(h)}{\mu}$$65 66    The mean of this biased distribution is $\mu + \sigma^2/\mu$, and you arrive uniformly within it, giving expected wait $(\mu + \sigma^2/\mu)/2$.67 68    The same phenomenon explains why the average class size experienced by a student exceeds the average class size reported by the university (large classes have more students to report them).69 70    ## Why "Inspector's Paradox"?71 72    The name comes from quality control, where an inspector arrives at a random time to sample a production process and systematically encounters longer-than-average intervals. The paradox is that a random observer is more likely to land inside a long gap than a short one, so their experienced mean interval exceeds the true mean interval. It feels paradoxical because you'd expect a random arrival to see the average gap, but length-biased sampling guarantees they see worse-than-average gaps whenever there's any variance at all.73    """)74    return75 76 77@app.cell(hide_code=True)78def _(mo):79    mo.md(r"""80    ## Implementation81 82    A `BusService` process generates buses under three headway distributions (regular, exponential, clustered bimodal) and records their arrival times. After the simulation, passenger wait times are estimated by sampling $N$ uniformly random arrival times and finding the next bus for each, without needing explicit `Passenger` processes.83    """)84    return85 86 87@app.cell(hide_code=True)88def _(mo):89    sim_time_slider = mo.ui.slider(90        start=0,91        stop=100_000,92        step=1_000,93        value=20_000,94        label="Simulation time",95    )96 97    mean_headway_slider = mo.ui.slider(98        start=5.0,99        stop=30.0,100        step=1.0,101        value=10.0,102        label="Mean headway",103    )104 105    seed_input = mo.ui.number(106        value=192,107        step=1,108        label="Random seed",109    )110 111    run_button = mo.ui.run_button(label="Run simulation")112 113    mo.vstack([114        sim_time_slider,115        mean_headway_slider,116        seed_input,117        run_button,118    ])119    return mean_headway_slider, seed_input, sim_time_slider120 121 122@app.cell123def _(mean_headway_slider, seed_input, sim_time_slider):124    SIM_TIME = int(sim_time_slider.value)125    MEAN_HEADWAY = float(mean_headway_slider.value)126    SEED = int(seed_input.value)127    N_PASSENGERS = 20_000128    return MEAN_HEADWAY, N_PASSENGERS, SEED, SIM_TIME129 130 131@app.cell132def _(MEAN_HEADWAY, Process, random):133    class BusService(Process):134        def init(self, mode, bus_arrivals):135            self.mode = mode136            self.bus_arrivals = bus_arrivals137 138        async def run(self):139            while True:140                if self.mode == "regular":141                    headway = MEAN_HEADWAY142                elif self.mode == "exponential":143                    headway = random.expovariate(1.0 / MEAN_HEADWAY)144                elif self.mode == "clustered":145                    headway = MEAN_HEADWAY * 0.2 if random.random() < 0.5 else MEAN_HEADWAY * 1.8146                else:147                    raise ValueError(f"Unknown mode: {self.mode}")148                await self.timeout(headway)149                self.bus_arrivals.append(self.now)150 151    return (BusService,)152 153 154@app.cell155def _(BusService, Environment, SIM_TIME):156    def collect_buses(mode):157        bus_arrivals = []158        env = Environment()159        BusService(env, mode, bus_arrivals)160        env.run(until=SIM_TIME)161        return bus_arrivals162 163    return (collect_buses,)164 165 166@app.cell167def _(N_PASSENGERS, random, statistics):168    def expected_wait(bus_arrivals, n=N_PASSENGERS):169        max_t = bus_arrivals[-1]170        waits = []171        for _ in range(n):172            t = random.uniform(0.0, max_t * 0.95)173            for b in bus_arrivals:174                if b > t:175                    waits.append(b - t)176                    break177        return statistics.mean(waits) if waits else 0.0178 179    return (expected_wait,)180 181 182@app.cell183def _(statistics):184    def headway_variance(bus_arrivals):185        headways = [b - a for a, b in zip(bus_arrivals, bus_arrivals[1:])]186        return statistics.variance(headways) if len(headways) > 1 else 0.0187 188    return (headway_variance,)189 190 191@app.cell(hide_code=True)192def _(MEAN_HEADWAY, mo):193    mu = MEAN_HEADWAY194    naive = MEAN_HEADWAY / 2.0195    var_exp = mu ** 2196    var_clustered = 0.5 * (mu * 0.2 - mu) ** 2 + 0.5 * (mu * 1.8 - mu) ** 2197    mo.md(f"""198    ## Results199 200    Mean headway: {MEAN_HEADWAY} → naive expected wait = {naive:.1f}201 202    - **Exponential** (Var ≈ {var_exp:.1f}): predicted = {mu / 2 + var_exp / (2 * mu):.1f} (= full mean headway!)203    - **Clustered** (Var ≈ {var_clustered:.1f}): predicted = {mu / 2 + var_clustered / (2 * mu):.1f}204    """)205    return (naive,)206 207 208@app.cell209def _(MEAN_HEADWAY, collect_buses, expected_wait, headway_variance, naive, pl):210    def run_models():211        rows = []212        for mode in ["regular", "exponential", "clustered"]:213            buses = collect_buses(mode)214            var_h = headway_variance(buses)215            mean_w = expected_wait(buses)216            rows.append({217                "mode": mode,218                "var_headway": round(var_h, 4),219                "mean_wait": round(mean_w, 4),220                "predicted": round(MEAN_HEADWAY / 2.0 + var_h / (2.0 * MEAN_HEADWAY), 4),221                "ratio": round(mean_w / naive, 4),222            })223        return pl.DataFrame(rows)224 225    return (run_models,)226 227 228@app.cell229def _(SEED, random, run_models):230    random.seed(SEED)231    df = run_models()232    df233    return (df,)234 235 236@app.cell237def _(alt, df, naive, pl):238    chart = (239        alt.Chart(df)240        .mark_bar()241        .encode(242            x=alt.X("mode:N", title="Bus schedule type"),243            y=alt.Y("mean_wait:Q", title="Mean passenger wait"),244            color=alt.Color("mode:N", legend=None),245            tooltip=["mode:N", "mean_wait:Q", "ratio:Q"],246        )247        .properties(title="Inspector's Paradox: Mean Wait by Schedule Type")248    )249    naive_line = (250        alt.Chart(pl.DataFrame({"naive": [naive]}))251        .mark_rule(strokeDash=[4, 4], color="gray")252        .encode(y="naive:Q")253    )254    (chart + naive_line)255    return256 257 258@app.cell(hide_code=True)259def _(mo):260    mo.md(r"""261    ## Understanding the Math262 263    ### Length-biased sampling264 265    Suppose buses run on an irregular schedule where gaps between buses are either 2 minutes or 18 minutes, each with probability 1/2. The mean gap is $\mu = (2 + 18)/2 = 10$ minutes. Now ask: if you arrive at a completely random moment, which gap are you most likely to land inside?266 267    A 2-minute gap occupies only 2 minutes on the clock, but an 18-minute gap occupies 18. Out of every 20 minutes of clock time on average, 2 minutes belong to a short gap and 18 to a long one. So a random arrival lands in a short gap with probability $2/(2+18) = 1/10$ and in a long gap with probability $18/20 = 9/10$. The expected gap length you experience is:268 269    $$E[\text{gap experienced}] = \frac{1}{10} \cdot 2 + \frac{9}{10} \cdot 18 = 0.2 + 16.2 = 16.4 \text{ minutes}$$270 271    That is far above the mean gap of 10 minutes. You are disproportionately likely to land inside a long gap simply because it takes up more time.272 273    ### The wait formula274 275    Once you are inside a gap, you arrive uniformly within it, so on average you land in the middle. Your expected wait is half the gap length you experience. The full formula is:276 277    $$E[\text{wait}] = \frac{\mu}{2} + \frac{\sigma^2}{2\mu}$$278 279    Here $\mu$ is the mean gap and $\sigma^2 = \text{Var}[\text{gap}]$ is the variance of gap lengths. The first term, $\mu/2$, is what you would get if every gap were exactly $\mu$ (deterministic buses — arrive in the middle every time). The second term, $\sigma^2/(2\mu)$, is the extra waiting from length-biased sampling. It is always non-negative, so irregular buses always make you wait longer than regular buses with the same mean headway.280 281    ### Why variance matters282 283    The variance $\sigma^2$ measures how spread out the gap sizes are. A perfectly regular bus schedule has $\sigma^2 = 0$ and gives the naive answer $\mu/2$. An exponentially distributed schedule has $\sigma^2 = \mu^2$, which doubles the expected wait to $\mu$. More irregular buses, higher penalty.284 285    ### Connecting to expected values286 287    The formula arises from a standard result: the expected length of the gap containing a random arrival is $\mu + \sigma^2/\mu$. You can think of this as the mean gap plus a correction term proportional to the variance divided by the mean. Dividing by 2 (uniform arrival within the gap) gives the wait formula above.288    """)289    return290 291 292if __name__ == "__main__":293    app.run()294