Probability & Bayes

The exact value is 2.28%, so the 68-95-99.7 rule gives a good approximation. This is used everywhere: 2σ events are "unlikely" (p ≈ 0.025 per tail), 3σ events are "rare" (p ≈ 0.0015).

The mean is np = 6, so getting 7 heads is only 1 above average — quite plausible at 21.5%. The binomial peaks near the mean and decays on both sides.

About 10% chance. The Poisson distribution models the count of rare, independent events in a fixed interval. It's the go-to for traffic rates, mutation counts, and photon arrivals.

Technically Binomial(500, p) works, but when n is large and p is small, the Poisson with λ = np is the standard approximation. The Poisson is the limit of the Binomial as n → ∞ and p → 0 with np = λ held constant. Both (B) and (C) are defensible, but Poisson is the canonical choice for "count of rare events in many trials."

So we expect 6 ± 1.55 heads. Using the Gaussian approximation: ~68% of the time we'd get between 4.45 and 7.55 heads, i.e., 5-7 heads. The Gaussian approximation works well when np and n(1-p) are both ≥ 5.

javascript
function gaussianPDF(x, mu, sigma) {
  const coeff = 1 / Math.sqrt(2 * Math.PI * sigma * sigma);
  const exponent = -(x - mu) * (x - mu) / (2 * sigma * sigma);
  return coeff * Math.exp(exponent);
}

The MLE for the Gaussian mean is simply the sample mean. This follows from taking the derivative of the log-likelihood with respect to μ, setting it to zero, and solving.

Note: the unbiased estimate uses 1/(n-1) = 1/4, giving σ² = 10 and σ = 3.16. MLE is slightly biased — it underestimates the true variance. For large n, the difference is negligible.

The Bernoulli MLE is simply the fraction of successes: p̂ = k/n. Intuitively obvious, but the derivation proves it's optimal — no other estimate makes the observed data more probable.

The Poisson MLE is the sample mean — just like the Gaussian MLE for μ. This makes sense: λ IS the mean of the Poisson distribution, so the best estimate of the mean is... the sample mean.

javascript
function mlEstimate(data) {
  const n = data.length;
  const mu = data.reduce((s, x) => s + x, 0) / n;
  const variance = data.reduce((s, x) => s + (x - mu) ** 2, 0) / n;
  return { mu, sigma: Math.sqrt(variance) };
}

MLE divides by n, unbiased by (n-1). The ratio is (n-1)/n. At n=2 the MLE is half the truth. At n=100 the MLE is 99% of the truth. At n=1000 it's 99.9%. In ML with millions of data points, the difference is utterly negligible — MLE is the standard choice.

The Beta(2,2) prior pulls the estimate from 0.70 toward 0.50 (the prior mean), landing at 0.667. The prior acts like 2 "phantom" observations of each outcome. Think of α-1 as prior heads and β-1 as prior tails.

The strong prior (equivalent to 38 phantom observations) overwhelms the 10 real observations, pulling the estimate close to 0.5. With Beta(20,20), you'd need ~40+ real observations to significantly move away from the prior belief of fairness.

MAP maximizes P(data|θ) × P(θ). If P(θ) is constant, the prior drops out and MAP reduces to maximizing P(data|θ) alone — which is exactly MLE. For the Beta prior: Beta(1,1) is uniform on [0,1], so MAP with Beta(1,1) gives p̂ = (k+0)/(n+0) = k/n = MLE.

Also: as n → ∞, the data overwhelms any finite prior, so MAP converges to MLE asymptotically. But they are only exactly equal when the prior is flat.

Small λ (weak regularization) = large σ (wide prior = vague belief about weight magnitudes). Large λ (strong regularization) = small σ (tight prior = strong belief weights should be near zero). This is the Bayesian interpretation of every regularizer you've ever used.

With n=100 data points, you need about 70-71 heads to get a MAP of 0.70. Compare to MLE: you'd need exactly 70 heads (70/100). The prior barely shifts the answer — with enough data, MAP and MLE converge. The 2 phantom observations are negligible against 100 real ones.

Position advances by velocity × dt: 10 + 2×1 = 12. Velocity stays at 2 (constant-velocity model). The predict step is just applying the physics model forward.

Position variance grew from 4 to 5.25. The extra 1 came from velocity uncertainty propagating into position uncertainty (the off-diagonal cross term), and 0.25 from process noise. Uncertainty always grows in the predict step — that's the whole reason we need measurements.

K = 0.724: we trust the measurement about 72% and our prediction about 28%. This makes sense — our prediction variance (5.25) is larger than our measurement noise (2), so we lean more toward the measurement.

K is literally a trust slider. R = 0 means perfect measurements, so K = 1 (just use the measurement). R = ∞ means useless measurements, so K = 0 (ignore them, trust the model). This is the single most important intuition for the Kalman filter.

The innovation (z - x̂^-) = -0.5 says the measurement is below the prediction. K = 0.724 says we trust the measurement more, so we move 72.4% of the way toward the measurement: 12 → 11.638 (closer to 11.5 than to 12).

The correct order is: Predict state → Predict covariance → Kalman gain → Update state → Update covariance.

Predict must happen before update (we need x̂^- and P^- before computing K). The Kalman gain must come before state/covariance update (K is needed for both). State and covariance updates can happen in either order since they're independent, but conventionally state comes first.

K = 0.8 means we put 80% weight on the measurement and 20% on the prediction. Our prediction is 4x noisier than our measurement (P^-/R = 4), so we lean heavily toward the measurement.

Variance dropped from 4 to 0.8 — an 80% reduction! Verify with information form: 1/P = 1/P^- + 1/R = 1/4 + 1/1 = 0.25 + 1 = 1.25, so P = 1/1.25 = 0.8. Fusing two sources always produces less uncertainty than either alone (0.8 < min(4, 1)).

The fused variance (3.45) is less than the better sensor alone (4). The GPS barely helps here (it's so noisy), but it still shaves off 0.55 from the LIDAR's variance. In information form, every additional sensor adds information (1/P), so fusion always helps.

With P^- > 0 and R > 0: the numerator P^- is strictly positive, the denominator P^- + R is strictly larger than P^-, so K = P^-/(P^-+R) is in (0, 1). It approaches 0 as R → ∞ and approaches 1 as R → 0, but never exactly reaches either (since both P^- and R are strictly positive, not zero or infinity).

javascript
function kalmanUpdate1D(xPred, pPred, z, R) {
  const K = pPred / (pPred + R);    // Kalman gain
  const x = xPred + K * (z - xPred); // updated state
  const p = (1 - K) * pPred;         // updated variance
  return { x, p };
}

Line 4 is the bug. It says (1 + K) instead of (1 - K). With (1+K), the variance increases at every update instead of shrinking. Since 0 < K < 1, multiplying by (1+K) gives a value between 1x and 2x of P^- — guaranteed divergence. The correct formula P = (1-K) × P^- always shrinks variance because (1-K) ∈ (0,1).

Seeing a "High" observation makes the Loaded state more likely (0.28 vs 0.12), even though the prior favored Fair (0.6 vs 0.4). The emission probability flipped the balance because Loaded strongly favors High rolls.

For completeness: α₂(Fair) = (0.12 × 0.7 + 0.28 × 0.4) × 0.2 = (0.084 + 0.112) × 0.2 = 0.196 × 0.2 = 0.0392. Two consecutive High rolls: α(Loaded) = 0.1428 dominates α(Fair) = 0.0392.

There is an 18.2% probability of seeing two consecutive High observations under this HMM. This value is used in model training (Baum-Welch) and model comparison (which HMM best explains the data?).

At each timestep t, for each of the N states j, we sum over all N previous states i: ∑_i α_t(i) × A[i][j]. That's N multiplications for each of N states = N² per timestep. Over T timesteps: O(N²T) = 25 × 200 = 5,000. Compare to brute-force enumeration of 5²⁰⁰ paths — a number with 140 digits.

javascript
function hmmForward(pi, A, B, obs) {
  const N = pi.length;
  // Initialize
  let alpha = pi.map((p, i) => p * B[i][obs[0]]);

  // Recurse
  for (let t = 1; t < obs.length; t++) {
    const newAlpha = [];
    for (let j = 0; j < N; j++) {
      let sum = 0;
      for (let i = 0; i < N; i++) sum += alpha[i] * A[i][j];
      newAlpha.push(sum * B[j][obs[t]]);
    }
    alpha = newAlpha;
  }

  // Terminate
  return alpha.reduce((s, a) => s + a, 0);
}

Viterbi initialization is identical to forward initialization. The divergence happens at the recursion step (max vs sum). Loaded is the most likely starting state given a High observation.

The best path to Loaded at t=2 came from Loaded at t=1 (backpointer = Loaded). So the Viterbi path so far is [Loaded, Loaded]. Two High observations in a row — the loaded die is the most plausible explanation.

Best final state: argmax(δ₃) = Fair (0.0376 > 0.0212). Backtrace: ψ₃(Fair) = Loaded (0.04704 > 0.02744). ψ₂(Loaded) = Loaded. So the Viterbi path is [Loaded, Loaded, Fair]. The casino used the loaded die for the two High rolls, then switched to the fair die for the Low roll. Intuitive!

Line 8 is the bug. It uses val += delta[i] * A[i][j] (summing — that's the forward algorithm!). It should use val = Math.max(val, delta[i] * A[i][j]) to find the maximum. The crucial difference between Forward and Viterbi is sum vs max at the recursion step.

Log-domain Viterbi: log δ_t+1(j) = max_i[log δ_t(i) + log A[i][j]] + log B[j][o_t+1]. Products become sums, max of products becomes max of sums. Values stay in a numerically stable range (e.g., -500 to 0 instead of 10^-217 to 10⁰). This is why log-probabilities are the default in all production HMM implementations.

A fair coin has maximum entropy for a binary variable: 1 bit. This means you need exactly 1 yes/no question to determine the outcome. It's the most uncertain a 2-outcome distribution can be.

Less than half a bit — much less uncertain than a fair coin. When outcomes are predictable (p near 0 or 1), entropy drops toward 0. Maximum entropy (1 bit) occurs at p = 0.5.

KL divergence measures how much information is lost when using q to approximate p. It's always ≥ 0 and equals 0 only when p = q.

KL(p‖q) = 0.368 but KL(q‖p) = 0.511 — they're different! KL divergence is not symmetric. This is why it's called a "divergence" not a "distance." The direction matters: KL(p‖q) penalizes places where p is high but q is low.

Since q = [0.5, 0.5], the cross-entropy simplifies to ln(2) regardless of p. When your model predicts uniform (maximum ignorance), the cross-entropy is always ln(number of classes). Training reduces this toward H(p), the irreducible entropy of the true distribution.

javascript
function entropy(p) {
  return -p.reduce((s, pi) =>
    s + (pi > 0 ? pi * Math.log(pi) : 0), 0);
}

function klDivergence(p, q) {
  return p.reduce((s, pi, i) =>
    s + (pi > 0 ? pi * Math.log(pi / q[i]) : 0), 0);
}

The first GPS reading is powerful: K = 0.8 because our initial uncertainty (100) dwarfs the GPS noise (25). After the first update, P drops from 100 to 20 — an 80% reduction in uncertainty from a single measurement.

K dropped from 0.80 to 0.44 — we're now trusting the prediction more relative to GPS. Variance continued to drop: 100 → 20 → 11.1. Each GPS fix reduces uncertainty further.

Steady-state variance = 0.495, much less than either sensor alone (GPS: 25, IMU accumulated drift). The steady-state Kalman gain is K_ss = (0.495 + 0.01) / (0.495 + 0.01 + 25) ≈ 0.505/25.505 ≈ 0.0198. At steady state, each GPS update only nudges by ~2% — the filter is highly confident.

Without GPS updates, each predict step adds Q to the variance: P_t = P_t-1 + Q. After 100 steps: P = 0.5 + 100 × 0.01 = 1.5. Uncertainty grows linearly with time — this is IMU drift. After 10,000 steps it would be 100.5. Without occasional GPS corrections, the estimate degrades. This is why real navigation systems have "GPS denied" alerts: the filter knows its own uncertainty is growing.

The correct order is: Initialize → IMU predict → Compute Kalman gain → GPS update → Output fused position.

In practice, the IMU runs at ~100 Hz (predict 100x per second), while GPS arrives at ~1 Hz. Between GPS updates, you run 100 predict steps. When a GPS measurement arrives, you compute K and update. The output is the fused estimate after each step (predict-only between GPS, full update when GPS is available).

javascript
function runKalman1D(x0, P0, Q, R, imuDisplacements, gpsReadings) {
  let x = x0, p = P0;
  for (let i = 0; i < imuDisplacements.length; i++) {
    // Predict
    x = x + imuDisplacements[i];
    p = p + Q;
    // Update
    const K = p / (p + R);
    x = x + K * (gpsReadings[i] - x);
    p = (1 - K) * p;
  }
  return { x, p };
}