#!/usr/bin/env python3 """ Does Mercury retrograde move the Amazon misspelling rate? This is the test called for in the 2026 postscripts to the Mercury retrograde articles -- the one never run in 2015. It cannot be run on the Reddit study's data: that analysis drew from a local SQL table ("May2015") and exported per-post fractions to fractionNMR.csv / fractionMR.csv, none of which survive. Nor is this the daily series behind Figure 1 of the Fourier article -- that is gone too. What survives is fullamazonmisspellingdata.csv: 5,295 consecutive days, 2 Jan 2000 - 1 Jul 2014, one row per day, carrying planetary longitudes, a retrograde flag for each of eight planets, and a misspelling RATE. It spans 49 Mercury retrograde periods. Calendar anchor verified: 10 of 11 known Mercury retrograde onsets across 2000-2014 land within a day of the flags in this file. THE IDEA -------- Do not ask "is the retrograde mean higher?" -- with enough rows something always is. Ask "is the REAL alignment unusual among arbitrary ones?" Rotate the retrograde flag against the rate series by every possible offset. Each rotation preserves exactly: - the autocorrelation of the rate series (lag-1 is +0.41, so this matters), - the number of retrograde days, - the block structure of retrograde periods (~21 days, ~3.2/year), changing only the PHASE -- the one thing under test. With 5,295 days every offset can be enumerated, so this is an exact test rather than a sample. An i.i.d. shuffle is also computed, only to show why it must not be used here. THE CONTROL ----------- The file flags all eight planets' retrogrades, not just Mercury's. If Mercury moves the rate, Mercury should stand out from the other seven. Run identically, it places 5th of 8. The only planet clearing p<0.05 is Uranus -- which is what a false positive looks like, and is also what the Mercury result looked like in 2015. USAGE ----- python3 mercury_permutation_test.py fullamazonmisspellingdata.csv Requires numpy. Runtime a few seconds. """ import csv import sys from datetime import date, timedelta import numpy as np START = date(2000, 1, 2) SEED = 20260716 PLANETS = {"MERRX": "Mercury", "VENRX": "Venus", "MARRX": "Mars", "JUPRX": "Jupiter", "SATRX": "Saturn", "URARX": "Uranus", "NEPRX": "Neptune", "PLURX": "Pluto"} def load(path): rows = list(csv.DictReader(open(path))) rate = np.array([float(r["RATE"]) for r in rows]) flags = {f: np.array([int(r[f]) for r in rows]) for f in PLANETS} dts = [START + timedelta(days=i) for i in range(len(rows))] return rate, flags, dts def mean_diff(series, lab): """Mean(series where flag=1) - Mean(series where flag=0).""" n, k = len(series), int(lab.sum()) s = series[lab == 1].sum() t = series.sum() return s / k - (t - s) / (n - k) def rotation_test(series, lab): """Exact: every circular offset of the flag against the series.""" obs = mean_diff(series, lab) nulls = np.array([mean_diff(series, np.roll(lab, o)) for o in range(1, len(series))]) p = (np.sum(np.abs(nulls) >= abs(obs)) + 1) / (len(nulls) + 1) return obs, nulls, p def iid_test(series, lab, B=20000): """Invalid for an autocorrelated series. Kept to show the size of the mistake.""" rng = np.random.default_rng(SEED) obs = mean_diff(series, lab) x = lab.copy() nulls = np.empty(B) for i in range(B): rng.shuffle(x) nulls[i] = mean_diff(series, x) return obs, (np.sum(np.abs(nulls) >= abs(obs)) + 1) / (B + 1) def predictor(path): """Can the sky predict the misspelling rate? And if so, how? Held out: 2008-2011, training on the years either side. NOT a train-on-the-past, predict-the-future split -- the rate halves across this file as spellcheck and autocorrect spread, so a forward split asks the model to forecast a technology rather than a planet. Holding out the middle avoids that, and a random sample is reported alongside it. Result: all 10 planets + 8 flags predict the RATE at R2 ~ 0.17 -- and predict THE CALENDAR DATE at R2 = 0.9999. Saturn, Neptune and Pluto never complete one circuit in 14.5 years, so their longitude is a date stamp. Keep only the bodies that genuinely recur (Sun..Mars) and the sky predicts nothing: R2 = -0.044. """ rows = list(csv.DictReader(open(path))) n = len(rows) dts = [START + timedelta(days=i) for i in range(n)] def col(nm): v = np.empty(n) for i, r in enumerate(rows): try: v[i] = float(r[nm]) except ValueError: # 15 cells are Mathematica RetrievalFailure expressions v[i] = np.nan ok = ~np.isnan(v) if (~ok).any(): v[~ok] = np.interp(np.flatnonzero(~ok), np.flatnonzero(ok), v[ok]) return v def feats(bodies, flags): F = [] for p in bodies: a = np.deg2rad(col(p)) F += [np.sin(a), np.cos(a)] for f in flags: F.append(col(f)) return np.column_stack(F) y = col("RATE") datenum = np.arange(n, dtype=float) ALL = ["SUN", "MOON", "MERCURY", "VENUS", "MARS", "JUPITER", "SATURN", "URANUS", "NEPTUNE", "PLUTO"] FAST = ["SUN", "MOON", "MERCURY", "VENUS", "MARS"] # these actually go round X_all = feats(ALL, list(PLANETS)) X_fast = feats(FAST, ["MERRX", "VENRX", "MARRX"]) idx = np.arange(n) a = next(i for i, d in enumerate(dts) if d >= date(2008, 1, 1)) b = next(i for i, d in enumerate(dts) if d >= date(2011, 1, 1)) te = idx[a:b] tr = np.concatenate([idx[:a], idx[b:]]) def r2(X, target, lam=1.0): mu, sd = X[tr].mean(0), X[tr].std(0) sd[sd == 0] = 1 A = np.column_stack([np.ones(len(tr)), (X[tr] - mu) / sd]) B = np.column_stack([np.ones(len(te)), (X[te] - mu) / sd]) I = np.eye(A.shape[1]) I[0, 0] = 0 w = np.linalg.solve(A.T @ A + lam * I, A.T @ target[tr]) return 1 - ((target[te] - B @ w) ** 2).sum() / ((target[te] - target[tr].mean()) ** 2).sum() return [("all planets -> misspelling rate", r2(X_all, y)), ("all planets -> THE CALENDAR DATE", r2(X_all, datenum)), ("fast bodies only -> misspelling rate", r2(X_fast, y))] def strawman_baseline(path): """The 2017 deepnet post trained on the first 80% of days and tested on the last 20%, reporting that its error beat methods "based on the mean (average) rate of the training data or an approach assuming random chance". The level of this series falls by more than a third between those two periods (0.293 -> 0.187), most plausibly because spellcheck and autocorrect matured over exactly those years. So a single constant carried forward from the end of training -- no astronomy, no network -- already scores R2 ~ 0.5 against that baseline. Beating "the training mean" here is not evidence; it is the minimum any method clears once it notices the level moved. """ rows = list(csv.DictReader(open(path))) y = np.array([float(r["RATE"]) for r in rows]) n = len(y) cut = int(n * 0.8) tr, te = y[:cut], y[cut:] base = tr.mean() denom = ((te - base) ** 2).sum() out = [("predict the training mean (their baseline)", 1 - ((te - base) ** 2).sum() / denom)] for w in (27, 90, 365): const = y[cut - w:cut].mean() out.append((f"predict a constant: last {w} days of training", 1 - ((te - const) ** 2).sum() / denom)) return tr.mean(), te.mean(), out def adjust(rate, dts): """Remove the weekday effect and the long-run yearly drift.""" grand = rate.mean() dow = {w: rate[[i for i in range(len(rate)) if dts[i].weekday() == w]].mean() for w in range(7)} yrs = {y: rate[[i for i in range(len(rate)) if dts[i].year == y]].mean() for y in {d.year for d in dts}} return np.array([rate[i] - dow[dts[i].weekday()] - yrs[dts[i].year] + 2 * grand for i in range(len(rate))]) def main(): path = sys.argv[1] if len(sys.argv) > 1 else "fullamazonmisspellingdata.csv" rate, flags, dts = load(path) n = len(rate) grand = rate.mean() mer = flags["MERRX"] print(f"{n:,} days: {dts[0]} to {dts[-1]}") blocks = int(np.sum((mer[1:] == 1) & (mer[:-1] == 0))) print(f"Mercury retrograde: {mer.sum()} days ({mer.sum()/n*100:.1f}%), " f"{blocks} blocks, {blocks/(n/365.25):.2f}/yr\n") c = rate - grand ac = lambda k: float(np.sum(c[:-k] * c[k:]) / ((n - k) * rate.var())) print(f"autocorrelation: lag1 {ac(1):+.3f} lag7 {ac(7):+.3f} lag30 {ac(30):+.3f}") print(" -> an i.i.d. shuffle is not a valid null for this series\n") _, p = iid_test(rate, mer) print(f"i.i.d. shuffle (INVALID, for contrast): p = {p:.4f}") for label, series in (("raw RATE", rate), ("weekday+year adjusted", adjust(rate, dts))): obs, nulls, p = rotation_test(series, mer) pct = float(np.mean(nulls < obs) * 100) lo, hi = np.percentile(nulls, [2.5, 97.5]) print(f"\nexact rotation test -- {label}") print(f" observed {obs:+.6f} ({obs/grand*100:+.3f}%) p = {p:.4f} percentile {pct:.1f}") print(f" null 95% band [{lo:+.6f}, {hi:+.6f}] -> observed is " f"{'INSIDE' if lo < obs < hi else 'OUTSIDE'}") dow = {w: rate[[i for i in range(n) if dts[i].weekday() == w]].mean() for w in range(7)} print(f"\nfor scale: weekday spread {max(dow.values())-min(dow.values()):+.5f} " f"vs Mercury {mean_diff(rate, mer):+.5f}") print("\nTHE CONTROL -- every planet in the file, tested identically:") print(f" {'planet':<9} {'observed':>10} {'% of mean':>10} {'p':>7}") res = [] for f, name in PLANETS.items(): obs, _, p = rotation_test(rate, flags[f]) res.append((name, obs, p)) for name, obs, p in sorted(res, key=lambda r: r[2]): mark = " <-- Mercury" if name == "Mercury" else "" print(f" {name:<9} {obs:>+10.5f} {obs/grand*100:>+9.2f}% {p:>7.3f}{mark}") rank = [r[0] for r in sorted(res, key=lambda r: r[2])].index("Mercury") + 1 print(f"\n Mercury places {rank} of 8.") print(f" With 8 tests, at least one clearing p<0.05 by luck alone is expected") print(f" {1-0.95**8:.0%} of the time; a Bonferroni threshold would be {0.05/8:.5f}.") print("\nTHE PREDICTOR -- held out 2008-2011, trained on the years either side:") for label, score in predictor(path): print(f" {label:<40} R2 = {score:>8.4f}") print("\n The sky reconstructs the calendar almost perfectly, and that is the") print(" whole of its apparent skill. Remove the slow planets -- whose longitude") print(" never repeats in 14.5 years and is therefore just a date -- and nothing") print(" is left.") trm, tem, base = strawman_baseline(path) print(f"\nTHE 2017 BASELINE -- first 80% train, last 20% test (their own split):") print(f" training mean rate {trm:.4f} -> test-period mean rate {tem:.4f} (the level moved)") for label, score in base: print(f" {label:<46} R2 = {score:>7.4f}") print("\n A constant carried forward from the end of training beats their baseline") print(" by a mile, with no astronomy in it at all. To count as evidence, the") print(" deepnet had to beat ~0.505 -- not 0.000.") if __name__ == "__main__": main()