# ---- chunk: packages ----
import sys, subprocess, importlib
needed = ["numpy", "scipy", "pandas", "plotnine", "matplotlib"]
to_install = [m for m in needed if importlib.util.find_spec(m) is None]
if to_install:
    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", *to_install])

import math, warnings
import numpy as np
import pandas as pd
from scipy import stats
from scipy.special import gammaln, logsumexp, comb
from scipy.optimize import brentq
from IPython.display import display, HTML, Markdown
from plotnine import *
warnings.filterwarnings("ignore")
pd.set_option("display.max_columns", 50)

# ---- chunk: helpers ----
# ---- colour roles (fixed per entity, never cycled) ---------------------------
col = {
    "Studentska lista":  "#2a78d6",   # blue
    "Ujedinjena Srbija": "#eb6834",   # orange
    "Manjinske liste":   "#1baf7a",   # aqua
    "Other":             "#898781",   # muted
    "RIK (official)":    "#2a78d6",
    "Ipsos (poll)":      "#eb6834",
}

# ---- a quiet ggplot theme ----------------------------------------------------
def theme_report(base_size=13, figure_size=(9, 5.5)):
    return (theme_minimal(base_size=base_size) +
            theme(plot_background=element_rect(fill="#fcfcfb", color="none"),
                  panel_grid_minor=element_blank(),
                  panel_grid_major=element_line(color="#e1e0d9", size=0.4),
                  axis_line_x=element_line(color="#c3c2b7", size=0.4),
                  axis_ticks=element_blank(),
                  text=element_text(color="#52514e"),
                  plot_title=element_text(color="#0b0b0b", weight="bold", size=base_size + 2, ha="left"),
                  plot_subtitle=element_text(color="#52514e", size=base_size - 1, ha="left"),
                  plot_caption=element_text(color="#898781", size=base_size - 3, ha="left"),
                  legend_position="top", legend_title=element_blank(),
                  strip_text=element_text(color="#0b0b0b", weight="bold"),
                  figure_size=figure_size, dpi=110))

# ---- formatting helpers -------------------------------------------------------
def fmt_n(x):   return f"{int(round(float(x))):,}"
def fmt_pct(x, d=2): return f"{float(x):.{d}f}%"
# turn a log10 probability into readable text, e.g. "2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros)"
# (both halves describe the same number, each rounded to one significant digit:
#  p ≈ m × 10^f  and  1/p ≈ lead × 10^e; note that m is NOT the same digit as lead unless m ≈ √10)
def fmt_log10p(l10):
    l10 = float(l10)
    if np.isnan(l10): return "not computable"
    if not np.isfinite(l10): return "below 10<sup>−300</sup> (underflow)"
    if l10 > -3: return f"{10 ** l10:.4f}"
    f = math.floor(l10);  m = round(10 ** (l10 - f))                     # p   = m × 10^f
    if m == 10: m, f = 1, f + 1
    e = math.floor(-l10); lead = round(10 ** (-l10 - e))                 # 1/p = lead × 10^e
    if lead == 10: lead, e = 1, e + 1
    return f"{m} × 10<sup>−{-f}</sup> (about 1 in a {lead} followed by {e} zeros)"
# a power of ten for prose, e.g. fmt_pow(-35.49) -> "10<sup>−35</sup>"
def fmt_pow(l10, d=0):
    s = f"{round(float(l10), d):.{d}f}".replace("-", "−")
    return f"10<sup>{s}</sup>"
# integer counts for a poll of size n from published percentages (sum to n)
def pct_to_counts(pct, n, remainder_to=None):
    pct = np.asarray(pct, dtype=float)
    k = np.rint(pct / 100 * n).astype(int)
    i = len(pct) - 1 if remainder_to is None else remainder_to
    k[i] += n - k.sum()
    return k
# log10 of a multinomial density, vectorised over the columns of a matrix (categories × polls)
def log10_dmultinom_mat(X, p):
    X = np.asarray(X, dtype=float); p = np.asarray(p, dtype=float)
    n = X.sum(axis=0)
    return (gammaln(n + 1) - gammaln(X + 1).sum(axis=0) + (X * np.log(p)[:, None]).sum(axis=0)) / np.log(10)
def log10_dmultinom(k, p):
    return float(log10_dmultinom_mat(np.asarray(k, dtype=float)[:, None], p)[0])
# the multinomial cumulative probability P(lower < N <= upper) for the rectangular region used in Method 3:
# exact, by Levin's representation written as a conditional-binomial sum (the role of pmultinom in R)
def pmultinom_rect(k1_max, k2_min, n, p):
    p1, p2 = p[0], p[1]
    n1 = np.arange(0, k1_max + 1)
    la = stats.binom.logpmf(n1, n, p1)
    lb = stats.binom.logsf(k2_min - 1, n - n1, p2 / (1 - p1))
    return float(np.exp(logsumexp(la + lb)))
def md(text): display(Markdown(text))
def kable_report(df, caption, align=None):
    # text columns left-aligned, numeric-looking columns right-aligned; HTML (superscripts) passed through
    import re
    def looks_num(s):
        s = s.astype(str)
        return bool((s.str.match(r"^[-−+]?[0-9.,%× ]*$") | s.str.contains("<sup>") | (s == "")).all() and s.str.contains("[0-9]").any())
    if align is None:
        align = ["r" if (pd.api.types.is_numeric_dtype(df[c]) or looks_num(df[c])) else "l" for c in df.columns]
    al = {"l": "left", "r": "right", "c": "center"}
    head = "".join(f'<th style="text-align:{al[a]}">{c}</th>' for c, a in zip(df.columns, align))
    rows = []
    for _, r in df.iterrows():
        rows.append("<tr>" + "".join(f'<td style="text-align:{al[a]}">{"" if (v is None or (isinstance(v, float) and np.isnan(v))) else v}</td>'
                                     for v, a in zip(r.values, align)) + "</tr>")
    display(HTML(f'<table class="table"><caption>{caption}</caption><thead><tr>{head}</tr></thead><tbody>{"".join(rows)}</tbody></table>'))
def fig(g, caption):
    display(g); display(HTML(f'<p class="caption">{caption}</p>'))
pct_lab = lambda l: [f"{v:g}%" for v in l]
comma   = lambda l: [f"{v:,.0f}" for v in l]
pow10   = lambda l: [f"$10^{{{int(v)}}}$" for v in l]

# ---- chunk: ipsos-data ----
ipsos = pd.DataFrame([
    ( 1, "Aleksandar Vučić – Ujedinjena Srbija",   "Ujedinjena Srbija",        "Ujedinjena Srbija", 47.4, 45.3, 49.5),
    ( 2, "Ivica Dačić – Faktor stabilnosti!",      "Faktor stabilnosti (SPS)", "Other",              3.2,  2.5,  3.9),
    ( 3, "Studentska lista – Studenti pobeđuju",   "Studentska lista",         "Studentska lista",  40.0, 37.9, 42.1),
    ( 4, "Rasim Ljajić – Ljudski (SDP – CDP)",     "Ljajić – Ljudski",         "Manjinske liste",    0.5,  0.2,  0.8),
    ( 5, "Savez vojvođanskih Mađara – dr Balint Pastor", "SVM",                "Manjinske liste",    1.2,  0.7,  1.7),
    ( 6, "Usame Zukorlić (SPP) Pokret za ujedinjenje", "Zukorlić (SPP)",       "Other",              0.3,  0.1,  0.5),
    ( 7, "Ruska stranka – Srbija u BRIKSU",        "Ruska stranka",            "Other",              0.2,  0.0,  0.4),
    ( 8, "Autentična desnica – dr Miloš Jovanović","Autentična desnica",       "Other",              1.7,  1.1,  2.3),
    ( 9, "Evropska Srbija – Marinika Tepić – Zdravko Ponoš", "Evropska Srbija", "Other",             0.8,  0.4,  1.2),
    (10, "SDA Sandžaka – Čuvari Sandžaka",         "SDA Sandžaka",             "Manjinske liste",    0.3,  0.1,  0.5),
    (11, "Bogoljub Karić – Pokret Snaga Srbije – BK", "Karić – PSS",           "Other",              0.6,  0.3,  0.9),
    (12, "Liberalno zelena stranka",               "Liberalno zelena",         "Other",              0.1,  0.0,  0.2),
    (13, "Biram pravdu – Biram budućnost – Koalicija Snaga", "Biram pravdu",   "Other",              0.1,  0.0,  0.2),
    (14, "Stranka rusko-srpskog jedinstva",        "Rusko-srpsko jedinstvo",   "Other",              0.2,  0.0,  0.4),
    (15, "Ćale ovo je za tebe – Petar Đurić",      "Ćale ovo je za tebe",      "Other",              0.2,  0.0,  0.4),
    (16, "Vladimir Gajić – Alternativa za Srbiju", "Gajić – Alternativa",      "Other",              0.3,  0.1,  0.5),
    (17, "Albanska koalicija – Koalicioni Shqiptar","Albanska koalicija",      "Manjinske liste",    0.3,  0.1,  0.5),
    (18, "Majka Srbija – Gruda",                   "Majka Srbija",             "Other",              0.1,  0.0,  0.2),
    (None, "(assumed invalid ballots)",            "Invalid ballots",          "Other",              2.5, None, None),
], columns=["no", "list", "short", "group", "ipsos_pct", "ci_lo", "ci_hi"])
ipsos["no"] = ipsos["no"].astype("Int64")
ipsos_n       = 3624
ipsos_turnout = {"est": 64.7, "lo": 62.9, "hi": 66.5}
assert abs(ipsos["ipsos_pct"].sum() - 100) < 1e-9

kable_report(ipsos.rename(columns={"no": "No.", "list": "Electoral list", "group": "Group", "ipsos_pct": "Ipsos %",
                                   "ci_lo": "95% interval, low", "ci_hi": "95% interval, high"})
             [["No.", "Electoral list", "Group", "Ipsos %", "95% interval, low", "95% interval, high"]],
             caption="Table 1. Ipsos final pre-election estimate (base: all ballots cast, incl. an assumed 2.5% invalid). The 95% intervals are Ipsos's own; what they mean is explained in Section 6. Some list names are abbreviated.")

# ---- chunk: rik-data ----
rik_meta = dict(
    snapshot       = "31 October 2026, 11:00 CET",
    stations_total = 8392, stations_done = 8339,
    registered     = 6451292,   # all 8,392 polling stations
    turned_out     = 4433332,   # voters who were handed a ballot, at the 8,339 processed stations
    turnout_pct    = 68.98,     # RIK's own turnout figure (relative to the registered voters at processed stations)
    ballots        = 4430659,   # ballots found in boxes (valid + invalid)
    valid          = 4374187,
    invalid        = 56472,
)

# national votes, in ballot order (no. 1..18)
rik_votes = np.array([1726554, 94571, 2349777, 20415, 35004, 12132, 6492, 39485, 12562,
                      12381, 20031, 2683, 1186, 4284, 2018, 3704, 27681, 3227])
assert rik_votes.sum() == rik_meta["valid"] and rik_meta["valid"] + rik_meta["invalid"] == rik_meta["ballots"]

# regional sub-totals Ipsos did NOT sample: Kosovo i Metohija and abroad
kim_votes    = np.array([25660, 1170, 7353, 13, 2, 10, 34, 360, 14, 12, 101, 6, 5, 23, 8, 10, 3, 7])
kim_invalid  = 863           # 35,654 ballots at 2.42% invalid
abroad_votes = np.array([19763, 386, 57756, 35, 17, 387, 29, 419, 278, 29, 78, 14, 3, 24, 11, 33, 8, 13])
abroad_invalid = 304         # RIK 73rd session, 29 Oct 2026: valid 79,283, invalid 304
assert kim_votes.sum() == 34791 and abroad_votes.sum() == 79283

rik = ipsos[ipsos["no"].notna()][["no", "list", "short", "group"]].copy()
rik["votes"], rik["kim"], rik["abroad"] = rik_votes, kim_votes, abroad_votes
rik = pd.concat([rik, pd.DataFrame([{"no": pd.NA, "list": "(invalid ballots)", "short": "Invalid ballots", "group": "Other",
                                      "votes": rik_meta["invalid"], "kim": kim_invalid, "abroad": abroad_invalid}])],
                ignore_index=True)
rik["no"] = rik["no"].astype("Int64")
rik["rik_pct"] = 100 * rik["votes"] / rik_meta["ballots"]

seats = {"Studentska lista": 141, "Ujedinjena Srbija": 103, "SVM": 2,
         "Albanska koalicija": 2, "Ljajić – Ljudski": 1, "SDA Sandžaka": 1}

kable_report(pd.DataFrame({"No.": rik["no"], "Electoral list": rik["list"], "Group": rik["group"],
                           "Votes": rik["votes"].map(fmt_n), "% of ballots": rik["rik_pct"].map(fmt_pct),
                           "Seats": rik["short"].map(seats).fillna(0).astype(int)}),
             caption=f"Table 2. RIK preliminary results, {rik_meta['snapshot']} ({fmt_n(rik_meta['stations_done'])} of {fmt_n(rik_meta['stations_total'])} polling stations). Seats are the preliminary allocation.")

# ---- chunk: rik-summary ----
kable_report(pd.DataFrame({
    "Quantity": ["Registered voters", "Voters who turned out", "Turnout",
                 "Ballots in boxes (valid + invalid)", "Valid ballots", "Invalid ballots", "Polling stations processed"],
    "Value": [fmt_n(rik_meta["registered"]), fmt_n(rik_meta["turned_out"]),
              f"{fmt_pct(rik_meta['turnout_pct'])} (RIK figure, relative to registered voters at processed stations)",
              f"{fmt_n(rik_meta['ballots'])} (the {fmt_n(rik_meta['turned_out'] - rik_meta['ballots'])} voters who were handed a ballot but did not deposit one are the difference from the row above)",
              fmt_n(rik_meta["valid"]),
              f"{fmt_n(rik_meta['invalid'])} ({fmt_pct(100 * rik_meta['invalid'] / rik_meta['ballots'])})",
              f"{fmt_n(rik_meta['stations_done'])} of {fmt_n(rik_meta['stations_total'])}"]}),
    caption="Table 3. Turnout and ballot totals (RIK, same snapshot).", align=["l", "l"])

# ---- chunk: fig-dumbbell ----
cmp = ipsos[["short", "group", "ipsos_pct"]].merge(rik[["short", "rik_pct"]], on="short")
cmp["panel"] = np.where(cmp["short"].isin(["Studentska lista", "Ujedinjena Srbija"]), "Two large\nlists", "All other lists\nand invalid\nballots")
cmp["panel"] = pd.Categorical(cmp["panel"], categories=["Two large\nlists", "All other lists\nand invalid\nballots"])
cmp["short"] = pd.Categorical(cmp["short"], categories=cmp.sort_values("rik_pct")["short"])   # reorder(short, rik_pct)

cmp_long = cmp.melt(id_vars=["short", "group", "panel"], value_vars=["ipsos_pct", "rik_pct"], var_name="source", value_name="pct")
cmp_long["source"] = pd.Categorical(cmp_long["source"].map({"ipsos_pct": "Ipsos (poll)", "rik_pct": "RIK (official)"}),
                                    categories=["RIK (official)", "Ipsos (poll)"])

g = (ggplot() +
     geom_segment(cmp, aes(x="ipsos_pct", xend="rik_pct", y="short", yend="short"), color="#c3c2b7", size=1.2) +
     geom_point(cmp_long, aes(x="pct", y="short", color="source"), size=3.2) +
     facet_grid("panel", scales="free", space="free") +
     scale_color_manual(values=col) +
     scale_x_continuous(labels=pct_lab) +
     labs(x="", y="", title="Where the poll put each list, and where the voters did",
          subtitle="Share of all ballots cast. Ipsos, 20–22 Oct 2026 (n = 3,624) vs RIK preliminary count, 31 Oct 2026") +
     theme_report(figure_size=(9, 7)) +
     theme(strip_text_y=element_text(rotation=0, ha="left")))
fig(g, "Figure 1. Ipsos estimate (orange) versus official result (blue) for every electoral list and for invalid ballots. The panels use different horizontal scales.")

# ---- chunk: scheme-a ----
grp_levels = ["Studentska lista", "Ujedinjena Srbija", "Manjinske liste", "Other"]

scheme_a = (ipsos.groupby("group", as_index=False)["ipsos_pct"].sum()
            .merge(rik.groupby("group", as_index=False)[["votes", "kim", "abroad"]].sum(), on="group", how="outer"))
scheme_a["group"] = pd.Categorical(scheme_a["group"], categories=grp_levels)
scheme_a = scheme_a.sort_values("group").reset_index(drop=True)
scheme_a["rik_pct"] = 100 * scheme_a["votes"] / rik_meta["ballots"]
scheme_a["rik_p"]   = scheme_a["votes"] / rik_meta["ballots"]

p_rik_A   = pd.Series(scheme_a["rik_p"].values, index=grp_levels)
pct_ips_A = pd.Series(scheme_a["ipsos_pct"].values, index=grp_levels)

kable_report(pd.DataFrame({
    "Category": grp_levels,
    "Lists included": ["Studentska lista – Studenti pobeđuju", "Aleksandar Vučić – Ujedinjena Srbija",
                       "SVM, Albanska koalicija, Ljajić – Ljudski, SDA Sandžaka (the four lists that won seats)",
                       "The remaining 12 lists plus invalid ballots"],
    "RIK votes": scheme_a["votes"].map(fmt_n), "RIK %": scheme_a["rik_pct"].map(fmt_pct),
    "Ipsos %": scheme_a["ipsos_pct"].map(lambda x: fmt_pct(x, 1)),
    "Gap (Ipsos − RIK)": (scheme_a["ipsos_pct"] - scheme_a["rik_pct"]).map(lambda x: f"{x:+.1f}".replace("-", "−"))}),
    caption="Table 4. Scheme A: the four-category comparison that drives the main analysis.")

# ---- chunk: fig-scheme-a ----
sa_long = scheme_a.melt(id_vars="group", value_vars=["ipsos_pct", "rik_pct"], var_name="source", value_name="pct")
sa_long["source"] = pd.Categorical(sa_long["source"].map({"ipsos_pct": "Ipsos (poll)", "rik_pct": "RIK (official)"}),
                                   categories=["RIK (official)", "Ipsos (poll)"])
sa_long["lbl"] = sa_long["pct"].map(lambda x: f"{x:.1f}")
g = (ggplot(sa_long, aes(x="group", y="pct", fill="source")) +
     geom_col(position=position_dodge(width=0.7), width=0.62) +
     geom_text(aes(label="lbl"), position=position_dodge(width=0.7), nudge_y=1.2, size=10, color="#52514e") +
     scale_fill_manual(values=col) +
     scale_y_continuous(labels=pct_lab, expand=(0, 0, 0.1, 0)) +
     labs(x="", y="", title="Scheme A: four categories", subtitle="Share of all ballots cast") +
     theme_report(figure_size=(9, 4.5)))
fig(g, "Figure 2. The four-category comparison. Ipsos (orange) against RIK (blue).")

# ---- chunk: scheme-b ----
scheme_b = rik[["short", "group", "votes"]].merge(ipsos[["short", "ipsos_pct"]], on="short")
scheme_b["rik_p"] = scheme_b["votes"] / rik_meta["ballots"]
p_rik_B   = pd.Series(scheme_b["rik_p"].values, index=scheme_b["short"])
pct_ips_B = pd.Series(scheme_b["ipsos_pct"].values, index=scheme_b["short"])
assert abs(p_rik_B.sum() - 1) < 1e-12 and abs(pct_ips_B.sum() - 100) < 1e-9

# ---- chunk: fig-sim-intro ----
rng_demo = np.random.default_rng(2026)
n_demo = 2000
sim_demo = rng_demo.multinomial(n_demo, p_rik_A.values, size=100_000).T      # categories × polls, as rmultinom
sim_df = pd.DataFrame({"Studentska lista":  100 * sim_demo[0] / n_demo,
                       "Ujedinjena Srbija": 100 * sim_demo[1] / n_demo}).melt(var_name="group", value_name="pct")

ref_lines = pd.DataFrame({"group": ["Studentska lista", "Ujedinjena Srbija"] * 2,
                          "what":  ["RIK (official)"] * 2 + ["Ipsos (poll)"] * 2,
                          "x":     list(100 * p_rik_A.values[:2]) + list(pct_ips_A.values[:2])})

g = (ggplot(sim_df, aes(x="pct")) +
     geom_histogram(aes(fill="group"), binwidth=0.25, color="#fcfcfb", size=0.2) +
     geom_vline(ref_lines, aes(xintercept="x", linetype="what", color="what"), size=0.9) +
     facet_wrap("group", scales="free_x") +
     scale_fill_manual(values=col, guide=None) +
     # the histograms carry the list colours, so the reference lines are kept neutral
     scale_color_manual(values={"RIK (official)": "#52514e", "Ipsos (poll)": "#0b0b0b"}) +
     scale_linetype_manual(values={"RIK (official)": "dashed", "Ipsos (poll)": "solid"}) +
     scale_x_continuous(labels=pct_lab, expand=(0, 1)) +
     scale_y_continuous(labels=comma) +
     labs(x="Share found by a simulated random poll", y="Number of simulated polls",
          title="What chance alone produces",
          subtitle="100,000 random polls of 2,000 decided voters, drawn from the actual RIK electorate") +
     theme_report(figure_size=(9, 5)))
fig(g, "Figure 3. Shares obtained by 100,000 simulated random polls of 2,000 decided voters from the actual electorate (RIK proportions). Histograms in each list's colour; the dashed line marks the true value, the solid black line what Ipsos published.")

# ---- chunk: sim-intro-range ----
rng = sim_df.groupby("group")["pct"].agg(lo="min", hi="max", sd="std").reset_index()

# ---- chunk: neff ----
neff_tbl = ipsos[ipsos["ci_lo"].notna() & (ipsos["ipsos_pct"] >= 0.5)].copy()     # tiny shares are rounded too coarsely to be informative
neff_tbl["p"] = neff_tbl["ipsos_pct"] / 100
neff_tbl["half_width"] = (neff_tbl["ci_hi"] - neff_tbl["ci_lo"]) / 2 / 100
neff_tbl["n_eff"] = neff_tbl["p"] * (1 - neff_tbl["p"]) / (neff_tbl["half_width"] / 1.96) ** 2
neff_tbl = neff_tbl.sort_values("ipsos_pct", ascending=False).reset_index(drop=True)

_t = {k: v / 100 for k, v in ipsos_turnout.items()}
n_turnout_eff = _t["est"] * (1 - _t["est"]) / (((_t["hi"] - _t["lo"]) / 2) / 1.96) ** 2
deff_full = ipsos_n / n_turnout_eff

n_primary = int(round(neff_tbl["n_eff"].median(), -1))

kable_report(pd.DataFrame({"Electoral list": neff_tbl["short"], "Ipsos %": neff_tbl["ipsos_pct"].map(lambda x: fmt_pct(x, 1)),
                           "95% interval": [f"{lo:.1f} – {hi:.1f}" for lo, hi in zip(neff_tbl["ci_lo"], neff_tbl["ci_hi"])],
                           "Half-width (pts)": (neff_tbl["half_width"] * 100).map(lambda x: f"{x:.1f}"),
                           "Implied n<sub>eff</sub>": neff_tbl["n_eff"].map(fmt_n)}),
             caption="Table 5. Effective sample size implied by each of Ipsos's published confidence intervals (lists at 0.5% or more).")

# ---- chunk: fig-neff ----
nt = neff_tbl.copy()
nt["short"] = pd.Categorical(nt["short"], categories=nt.sort_values("n_eff")["short"])   # reorder(short, n_eff)
g = (ggplot(nt, aes(x="short", y="n_eff")) +
     geom_segment(aes(xend="short", y=0, yend="n_eff"), color="#c3c2b7", size=1) +
     geom_point(size=3.4, color="#2a78d6") +
     geom_hline(yintercept=n_primary, color="#0b0b0b", size=0.6) +
     geom_hline(yintercept=ipsos_n, color="#898781", linetype="dashed", size=0.6) +
     annotate("text", x=len(nt) + 0.45, y=n_primary * 1.01, label=f"adopted n = {fmt_n(n_primary)}", ha="left", va="top", size=10, color="#0b0b0b") +
     annotate("text", x=0.55, y=ipsos_n * 0.99, label=f"interviews reported = {fmt_n(ipsos_n)}", ha="right", va="bottom", size=10, color="#52514e") +
     coord_flip() +
     scale_y_continuous(labels=comma, limits=(0, ipsos_n * 1.05)) +
     labs(x="", y="Implied effective sample size", title="How many respondents do Ipsos's confidence intervals imply?") +
     theme_report(figure_size=(9, 4.5)))
fig(g, "Figure 4. Effective sample size implied by each published interval. The solid line is the median, which we adopt as the primary n; the dashed line is the number of interviews Ipsos reports.")

# ---- chunk: ipsos-counts ----
k_ips_A = pd.Series(pct_to_counts(pct_ips_A.values, n_primary, remainder_to=3), index=grp_levels)   # any rounding remainder goes to "Other"
expected_A = n_primary * p_rik_A
sd_A = np.sqrt(n_primary * p_rik_A * (1 - p_rik_A))

kable_report(pd.DataFrame({"Category": grp_levels,
                           "Ipsos count (out of n)": k_ips_A.map(fmt_n).values,
                           "Expected under H₀": expected_A.map(lambda x: f"{x:,.1f}").values,
                           "Difference": (k_ips_A - expected_A).map(lambda x: f"{x:+.0f}").values,
                           "In standard deviations": ((k_ips_A - expected_A) / sd_A).map(lambda x: f"{x:+.1f}").values}),
             caption=f"Table 6. Ipsos's result as counts out of n = {fmt_n(n_primary)}, against what the RIK electorate would produce on average.")

# ---- chunk: m1-point ----
log10_point_ips = log10_dmultinom(k_ips_A.values, p_rik_A.values)
# (approximately) the most likely single outcome: the rounded expected counts
k_mode = pct_to_counts(100 * p_rik_A.values, n_primary, remainder_to=3)
log10_point_mode = log10_dmultinom(k_mode, p_rik_A.values)

# ---- chunk: m2-binomial ----
def binom_tail(k, n, p, direction="low"):
    if direction == "low":
        return stats.binom.logcdf(k, n, p) / np.log(10)
    return stats.binom.logsf(k - 1, n, p) / np.log(10)

m2 = pd.DataFrame({"Category": ["Studentska lista", "Ujedinjena Srbija"], "direction": ["low", "high"],
                   "k": k_ips_A.values[:2], "p": p_rik_A.values[:2]})
m2["expected"] = n_primary * m2["p"]
m2["sd"] = np.sqrt(n_primary * m2["p"] * (1 - m2["p"]))
m2["z"] = (m2["k"] - m2["expected"]) / m2["sd"]
m2["log10_p"] = [binom_tail(k, n_primary, p, d) for k, p, d in zip(m2["k"], m2["p"], m2["direction"])]
m2["question"] = np.where(m2["direction"] == "low", "Share this low or lower?", "Share this high or higher?")
kable_report(pd.DataFrame({"Category": m2["Category"], "Ipsos count": m2["k"].map(fmt_n), "Expected": m2["expected"].map(lambda x: f"{x:,.1f}"),
                           "z-score": m2["z"].map(lambda x: f"{x:+.1f}".replace("-", "−")),
                           "Question": m2["question"], "Probability": m2["log10_p"].map(fmt_log10p)}),
             caption="Table 7. Binomial tail probabilities for the two large lists, taken one at a time.")

# ---- chunk: m3-pmultinom ----
# pmultinom: P(lower < N <= upper); lower is exclusive, upper inclusive
# (Python: the same rectangle, P(N1 <= k1 and N2 >= k2), computed exactly by the conditional-binomial form of Levin's representation)
p_norm = p_rik_A.values / p_rik_A.values.sum()
joint_exact = pmultinom_rect(k1_max=int(k_ips_A.iloc[0]), k2_min=int(k_ips_A.iloc[1]), n=n_primary, p=p_norm)

# an independent check by brute-force simulation
rng_A = np.random.default_rng(2026)
n_sim = 2_000_000
sim_A = rng_A.multinomial(n_primary, p_rik_A.values, size=n_sim).T          # categories × polls
joint_mc = np.mean((sim_A[0] <= k_ips_A.iloc[0]) & (sim_A[1] >= k_ips_A.iloc[1]))
# correlation of the two counts under the multinomial
rho = -np.sqrt(p_rik_A.iloc[0] * p_rik_A.iloc[1] / ((1 - p_rik_A.iloc[0]) * (1 - p_rik_A.iloc[1])))
# given Studentska's count, where would Ujedinjena's count be expected to land?
us_cond_p   = p_rik_A.iloc[1] / (1 - p_rik_A.iloc[0])
us_cond_exp = (n_primary - k_ips_A.iloc[0]) * us_cond_p
us_cond_sd  = np.sqrt((n_primary - k_ips_A.iloc[0]) * us_cond_p * (1 - us_cond_p))

# ---- chunk: m3-table ----
kable_report(pd.DataFrame({
    "Approach": ["exact algorithm (Levin's representation; pmultinom in R)", f"Monte Carlo, {fmt_n(n_sim)} simulated polls"],
    "P(Studentska ≤ 40.0% and Ujedinjena ≥ 47.4%)": [
        fmt_log10p(np.log10(joint_exact)),
        f"0 of {fmt_n(n_sim)} (so below {fmt_log10p(np.log10(1 / n_sim))})" if joint_mc == 0 else f"{joint_mc:.2e}"]}),
    caption="Table 8. The joint tail probability computed two ways.")

# ---- chunk: m4-exact-enum ----
# The exact test by exhaustive enumeration: visit every count vector (n1, n2, n3, n4) with
# n1 + n2 + n3 + n4 = n, and add up P(vector | H0) over those no more probable than Ipsos's.
# One loop over n1; for each n1 the (n2, n3) grid is handled as a matrix, n4 being what is left.
def exact_tail_enum(k_obs, p, n):
    k_obs = np.asarray(k_obs, dtype=float); p = np.asarray(p, dtype=float)
    thr   = gammaln(n + 1) - gammaln(k_obs + 1).sum() + (k_obs * np.log(p)).sum()   # log P(Ipsos's vector)
    x     = np.arange(n + 1)
    a     = [x * np.log(p[i]) - gammaln(x + 1) for i in range(4)]                   # per-category pieces of log P
    const = gammaln(n + 1)
    E     = n * p
    tot = corner = 0.0
    for n1 in range(n + 1):
        m   = n - n1
        idx = np.arange(m + 1)
        lp  = const + a[0][n1] + a[1][idx][:, None] + a[2][idx][None, :]
        n4  = m - idx[:, None] - idx[None, :]
        ok  = n4 >= 0                                        # n4 must be a count
        lp  = np.where(ok, lp + a[3][np.clip(n4, 0, n)], -np.inf)
        sel = ok & (lp <= thr)                               # at least as improbable as Ipsos's
        tot += np.exp(lp[sel]).sum()
        # the part of that mass in the "Ipsos-like" corner: Studentska below and Ujedinjena above expectation
        if n1 < E[0]:
            corner += np.exp(lp[sel & (idx[:, None] >= E[1])]).sum()
    return {"p": tot, "corner_share": corner / tot, "n_outcomes": comb(n + 3, 3, exact=True)}

enum = exact_tail_enum(k_ips_A.values, p_rik_A.values, n_primary)
p_exact       = enum["p"]
log10_p_exact = np.log10(p_exact)

# ---- chunk: m4-exact-mc ----
log10_sim = log10_dmultinom_mat(sim_A, p_rik_A.values)
n_as_extreme = int(np.sum(log10_sim <= log10_point_ips))

# importance sampling: draw from an electorate centred on Ipsos's counts, re-weight to H0.
# The weights are heavy-tailed, so instead of trusting a single standard error we run
# several independent replicates and report the spread of their exponents.
rng_is = np.random.default_rng(2026)
n_rep, n_per = 10, 1_000_000; n_is = n_rep * n_per
q_prop = k_ips_A.values / n_primary
p_is_rep = np.empty(n_rep)
for i in range(n_rep):
    sim_is = rng_is.multinomial(n_primary, q_prop, size=n_per).T
    l10_h0 = log10_dmultinom_mat(sim_is, p_rik_A.values)   # log10 P(poll | real electorate)
    l10_q  = log10_dmultinom_mat(sim_is, q_prop)           # log10 P(poll | made-up electorate)
    p_is_rep[i] = np.mean(10 ** (l10_h0 - l10_q) * (l10_h0 <= log10_point_ips))
p_is       = p_is_rep.mean()
log10_p_is = np.log10(p_is)
is_spread  = (np.log10(p_is_rep).min(), np.log10(p_is_rep).max())   # exponents of the individual replicates

# ---- chunk: fig-m4 ----
hist_df = pd.DataFrame({"l10": log10_sim})
ymax = np.histogram(log10_sim, bins=np.arange(np.floor(log10_sim.min()), np.ceil(log10_sim.max()) + 0.25, 0.25))[0].max()
xmid = (log10_point_ips + log10_sim.min()) / 2
g = (ggplot(hist_df, aes(x="l10")) +
     geom_histogram(binwidth=0.25, fill="#2a78d6", color="#fcfcfb", size=0.2) +
     geom_vline(xintercept=log10_point_ips, color="#eb6834", size=1) +
     annotate("label", x=log10_point_ips + 0.3, y=ymax * 0.96, label=f"Ipsos's result: $10^{{{log10_point_ips:.0f}}}$",
              ha="left", va="top", color="#eb6834", fill="#fcfcfb", label_size=0, size=10.5) +
     annotate("text", x=xmid, y=ymax * 0.75, label=f"about $10^{{{log10_sim.min() - log10_point_ips:.0f}}}$ times less likely than",
              color="#52514e", size=9.5) +
     annotate("text", x=xmid, y=ymax * 0.68, label="the least likely of two million simulated polls", color="#52514e", size=9.5) +
     scale_x_continuous(labels=pow10) +
     scale_y_continuous(labels=comma) +
     labs(x="Probability of the simulated poll's result (log scale)", y="Number of simulated polls",
          title="How unusual is Ipsos's result, measured by its own probability?",
          subtitle="Two million random polls from the RIK electorate (Scheme A, n as adopted)") +
     theme_report(figure_size=(9, 4.5)))
fig(g, "Figure 5. Log-probability (base 10) of two million simulated random polls under H₀. Ipsos's result lies far to the left of every one of them.")

# ---- chunk: m5-gtest ----
def g_test(obs, p):
    obs = np.asarray(obs, dtype=float); p = np.asarray(p, dtype=float)
    n = obs.sum(); E = n * p
    G  = 2 * np.sum(np.where(obs > 0, obs * np.log(np.where(obs > 0, obs, 1) / E), 0))
    X2 = np.sum((obs - E) ** 2 / E)
    df = len(obs) - 1
    return {"G": G, "X2": X2, "df": df,
            "log10_p_G":  stats.chi2.logsf(G,  df) / np.log(10),
            "log10_p_X2": stats.chi2.logsf(X2, df) / np.log(10)}
gA = g_test(k_ips_A.values, p_rik_A.values)

# Scheme B: all 19 categories, merging those whose expected count is below 5 (the usual rule of thumb)
k_ips_B = pd.Series(pct_to_counts(pct_ips_B.values, n_primary, remainder_to=list(pct_ips_B.index).index("Invalid ballots")), index=pct_ips_B.index)
small = n_primary * p_rik_B < 5
k_ips_B2 = pd.concat([k_ips_B[~small], pd.Series({"Small lists (merged)": k_ips_B[small].sum()})])
p_rik_B2 = pd.concat([p_rik_B[~small], pd.Series({"Small lists (merged)": p_rik_B[small].sum()})])
gB = g_test(k_ips_B2.values, p_rik_B2.values)

kable_report(pd.DataFrame({
    "Scheme": ["A: four categories", f"B: all lists ({len(k_ips_B2)} categories after merging {int(small.sum())} tiny ones)"],
    "G statistic": [f"{gA['G']:.0f}", f"{gB['G']:.0f}"],
    "Degrees of freedom": [gA["df"], gB["df"]],
    "p-value (G)": [fmt_log10p(gA["log10_p_G"]), fmt_log10p(gB["log10_p_G"])],
    "p-value (Pearson χ², less reliable this far out)": [fmt_log10p(gA["log10_p_X2"]), fmt_log10p(gB["log10_p_X2"])]}),
    caption="Table 9. Likelihood-ratio and Pearson chi-square tests of Ipsos's result against the RIK electorate (asymptotic approximations; orders of magnitude).")

# ---- chunk: methods-summary ----
summary_tbl = pd.DataFrame({
    "Method": ["Method 1 – Point probability of Ipsos's exact vector (for reference; not a tail)",
               "Method 2a – Binomial tail: Studentska lista this low or lower (exact)",
               "Method 2b – Binomial tail: Ujedinjena Srbija this high or higher (exact)",
               "Method 3 – Joint tail: both at once (exact; plain Monte Carlo)",
               "Method 4a – Exact multinomial test, exhaustive enumeration (exact; the headline)",
               "Method 4b – Exact multinomial test, plain Monte Carlo",
               "Method 4c – Exact multinomial test, importance sampling (biased low here; see text)",
               "Method 5a – Likelihood-ratio G test, Scheme A (chi-square approximation)",
               "Method 5b – Likelihood-ratio G test, Scheme B, all lists (chi-square approximation)"],
    "Probability under random sampling": [
        fmt_log10p(log10_point_ips), fmt_log10p(m2.log10_p[0]), fmt_log10p(m2.log10_p[1]),
        f"exact: {fmt_log10p(np.log10(joint_exact))}; Monte Carlo: 0 of {fmt_n(n_sim)}",
        f"<b>{fmt_log10p(log10_p_exact)}</b>",
        f"below {fmt_log10p(np.log10(1 / n_sim))} (0 of {fmt_n(n_sim)})" if n_as_extreme == 0 else f"{n_as_extreme / n_sim:.2e}",
        f"{fmt_log10p(log10_p_is)} ({n_rep} runs: exponents {is_spread[0]:.1f} to {is_spread[1]:.1f})",
        fmt_log10p(gA["log10_p_G"]), fmt_log10p(gB["log10_p_G"])]})
kable_report(summary_tbl, caption=f"Table 10. Summary of all methods, Scheme A unless stated, n = {fmt_n(n_primary)}.", align=["l", "l"])
# the range of tail probabilities across the methods that can reach the tail (the importance-sampling
# estimate is excluded: Method 4 shows it to be biased low)
tail_range = (min(m2.log10_p[0], np.log10(joint_exact), log10_p_exact, gA["log10_p_G"], gB["log10_p_G"]),
              max(m2.log10_p[0], np.log10(joint_exact), log10_p_exact, gA["log10_p_G"], gB["log10_p_G"]))

# ---- chunk: own-precision ----
own = pd.DataFrame({"Category": ["Studentska lista", "Ujedinjena Srbija"],
                    "ipsos": [40.0, 47.4], "lo": [37.9, 45.3], "hi": [42.1, 49.5],
                    "actual": 100 * p_rik_A.values[:2]})
own["se_own"] = (own["hi"] - own["lo"]) / 2 / 1.96
own["outside_by"] = np.where(own["actual"] > own["hi"], own["actual"] - own["hi"], own["lo"] - own["actual"])
own["z_own"] = (own["actual"] - own["ipsos"]) / own["se_own"]
own["log10_p_own"] = stats.norm.logcdf(-np.abs(own["z_own"])) / np.log(10)
kable_report(pd.DataFrame({"Category": own["Category"], "Ipsos 95% interval": [f"{lo:.1f} – {hi:.1f}" for lo, hi in zip(own["lo"], own["hi"])],
                           "Actual result": own["actual"].map(fmt_pct), "Outside the interval by (pts)": own["outside_by"].map(lambda x: f"{x:.1f}"),
                           "Ipsos's own SE (pts)": own["se_own"].map(lambda x: f"{x:.2f}"), "Distance in SEs": own["z_own"].abs().map(lambda x: f"{x:.1f}"),
                           "One-sided probability": own["log10_p_own"].map(fmt_log10p)}),
             caption="Table 11. The actual result measured against Ipsos's own published margins of error.")

# ---- chunk: fig-cloud ----
cloud = pd.DataFrame({"sl": 100 * sim_A[0, :20000] / n_primary, "us": 100 * sim_A[1, :20000] / n_primary})
points = pd.DataFrame({
    "who": ["RIK (official)", "Ipsos (poll)"],
    "sl":  [100 * p_rik_A.iloc[0], 40.0],
    "us":  [100 * p_rik_A.iloc[1], 47.4],
    "note": ["Official result\n(RIK, preliminary)",
             f"Ipsos, fieldwork 20–22 Oct\n{fmt_n(ipsos_n)} interviews, effective n ≈ {fmt_n(n_primary)}"],
    # label positions: RIK below-left of the cloud, Ipsos to the right of its box
    "hjust": ["right", "left"], "lx": [52.3, 42.6], "ly": [36.2, 47.4]})
points["lbl"] = [f"{n}\n({s:.1f}%, {u:.1f}%)" for n, s, u in zip(points["note"], points["sl"], points["us"])]
points["who"] = pd.Categorical(points["who"], categories=["RIK (official)", "Ipsos (poll)"])

g = (ggplot() +
     geom_point(cloud, aes(x="sl", y="us"), color="#898781", alpha=0.15, size=0.9) +
     annotate("rect", xmin=37.9, xmax=42.1, ymin=45.3, ymax=49.5, fill="none", color="#eb6834", linetype="dashed", size=0.6) +
     annotate("text", x=42.1, y=49.75, label="Ipsos's own 95% intervals", ha="right", va="bottom", size=9.5, color="#eb6834") +
     geom_point(points, aes(x="sl", y="us", color="who", shape="who"), size=4, stroke=1.6) +
     geom_text(points, aes(x="lx", y="ly", label="lbl", ha="hjust"), size=9.5, color="#52514e", lineheight=0.95) +
     scale_color_manual(values=col) +
     scale_shape_manual(values={"RIK (official)": "x", "Ipsos (poll)": "o"}) +
     scale_x_continuous(labels=pct_lab) + scale_y_continuous(labels=pct_lab) +
     labs(x="Studentska lista share", y="Ujedinjena Srbija share",
          title="Everything chance can produce, and where Ipsos landed",
          subtitle=f"Grey: 20,000 random polls of n = {fmt_n(n_primary)} from the RIK electorate. Each point is labelled (Studentska, Ujedinjena).") +
     theme_report(figure_size=(9, 6.5)))
fig(g, "Figure 6. Where ideal random polls land (grey cloud, 20,000 simulations from the RIK electorate), Ipsos's published result (orange), Ipsos's own 95% box (orange dashed rectangle), and the official result (blue cross).")

# ---- chunk: sens-n ----
n_grid = [1500, n_primary, int(round(0.8 * ipsos_n, -2)), ipsos_n]
rows = []
for n in n_grid:
    k = pct_to_counts(pct_ips_A.values, n, 3)
    g_ = g_test(k, p_rik_A.values)
    rows.append({"n": n, "G test (Scheme A)": g_["log10_p_G"],
                 "Studentska this low":  binom_tail(k[0], n, p_rik_A.iloc[0], "low"),
                 "Ujedinjena this high": binom_tail(k[1], n, p_rik_A.iloc[1], "high")})
sens_n = pd.DataFrame(rows)
kable_report(pd.DataFrame({"Assumed n": sens_n["n"].map(fmt_n),
                           **{c: sens_n[c].map(fmt_log10p) for c in sens_n.columns[1:]}}),
             caption="Table 12. Probability of Ipsos's result for different assumed sample sizes.")

# ---- chunk: fig-sens-n ----
sl = sens_n.melt(id_vars="n", var_name="method", value_name="log10_p")
sl["method"] = pd.Categorical(sl["method"], categories=["G test (Scheme A)", "Studentska this low", "Ujedinjena this high"])
g = (ggplot(sl, aes(x="n", y="log10_p", color="method")) +
     geom_hline(yintercept=np.log10(0.05), linetype="dotted", color="#898781") +
     geom_hline(yintercept=-6, linetype="dotted", color="#898781") +
     annotate("text", x=1500, y=np.log10(0.05) + 1.2, label="p = 0.05", ha="left", va="bottom", size=10, color="#52514e") +
     annotate("text", x=1500, y=-6 + 1.2, label="p = one in a million", ha="left", va="bottom", size=10, color="#52514e") +
     geom_line(size=1) + geom_point(size=3) +
     scale_color_manual(values={"G test (Scheme A)": "#0b0b0b", "Studentska this low": "#2a78d6", "Ujedinjena this high": "#eb6834"}) +
     scale_x_continuous(breaks=n_grid, labels=comma) +
     scale_y_continuous(labels=pow10, expand=(0.02, 0, 0.1, 0)) +
     labs(x="Assumed effective sample size", y="Probability under random sampling",
          title="The conclusion does not depend on the sample size we assume") +
     theme_report(figure_size=(9, 4.5)))
fig(g, "Figure 7. Log-probability of Ipsos's result as a function of the assumed effective sample size. Every line is far below the conventional 5% and 1-in-a-million thresholds at every n.")

# ---- chunk: sens-pop ----
scheme_a_adj = scheme_a.copy()
scheme_a_adj["votes_adj"] = scheme_a_adj["votes"] - scheme_a_adj["kim"] - scheme_a_adj["abroad"]
scheme_a_adj["rik_p_adj"] = scheme_a_adj["votes_adj"] / scheme_a_adj["votes_adj"].sum()
scheme_a_adj["rik_pct_adj"] = 100 * scheme_a_adj["rik_p_adj"]
p_rik_A_adj = pd.Series(scheme_a_adj["rik_p_adj"].values, index=grp_levels)
g_adj = g_test(k_ips_A.values, p_rik_A_adj.values)
bin_adj = [binom_tail(k_ips_A.iloc[0], n_primary, p_rik_A_adj.iloc[0], "low"),
           binom_tail(k_ips_A.iloc[1], n_primary, p_rik_A_adj.iloc[1], "high")]

kable_report(pd.DataFrame({"Category": grp_levels, "RIK % (all voters)": scheme_a_adj["rik_pct"].map(fmt_pct),
                           "RIK % (Serbia proper, without Kosovo i Metohija and abroad)": scheme_a_adj["rik_pct_adj"].map(fmt_pct),
                           "Ipsos %": scheme_a_adj["ipsos_pct"].map(lambda x: fmt_pct(x, 1))}),
             caption=f"Table 13. The population redefined to match Ipsos's sampling frame (removes {fmt_n((scheme_a['kim'] + scheme_a['abroad']).sum())} ballots).")

# ---- chunk: sens-valid ----
# Scheme A on the valid-vote base: drop the invalid row from both sources
ips_valid = (ipsos[ipsos["short"] != "Invalid ballots"].groupby("group")["ipsos_pct"].sum() / 0.975).reindex(grp_levels)
rik_valid = rik[rik["short"] != "Invalid ballots"].groupby("group")["votes"].sum().reindex(grp_levels)
p_rik_A_valid   = rik_valid / rik_valid.sum()
pct_ips_A_valid = ips_valid
k_ips_A_valid   = pct_to_counts(pct_ips_A_valid.values, n_primary, 3)
g_valid   = g_test(k_ips_A_valid, p_rik_A_valid.values)
bin_valid = binom_tail(k_ips_A_valid[0], n_primary, p_rik_A_valid.iloc[0], "low")

# Scheme B on the valid-vote base
sb_valid = scheme_b[scheme_b["short"] != "Invalid ballots"].copy()
sb_valid["p"] = sb_valid["votes"] / sb_valid["votes"].sum()
sb_valid["pct"] = sb_valid["ipsos_pct"] / 0.975
k_B_valid = pd.Series(pct_to_counts(sb_valid["pct"].values, n_primary, int(np.argmax(sb_valid["pct"].values))), index=sb_valid["short"])
p_B_valid = pd.Series(sb_valid["p"].values, index=sb_valid["short"])
small_v = n_primary * p_B_valid < 5
k_B_valid2 = pd.concat([k_B_valid[~small_v], pd.Series({"Small lists (merged)": k_B_valid[small_v].sum()})])
p_B_valid2 = pd.concat([p_B_valid[~small_v], pd.Series({"Small lists (merged)": p_B_valid[small_v].sum()})])
gB_valid = g_test(k_B_valid2.values, p_B_valid2.values)

kable_report(pd.DataFrame({
    "Test": ["Studentska lista this low or lower (binomial)", "G test, Scheme A", "G test, Scheme B"],
    "All ballots (main analysis)": [fmt_log10p(x) for x in (m2.log10_p[0], gA["log10_p_G"], gB["log10_p_G"])],
    "Valid votes only": [fmt_log10p(x) for x in (bin_valid, g_valid["log10_p_G"], gB_valid["log10_p_G"])]}),
    caption="Table 14. The key tests with invalid ballots removed from both sources.")
# the overall range of tail probabilities across methods AND bases, used in the conclusions
tail_range_all = (min(tail_range[0], g_valid["log10_p_G"], gB_valid["log10_p_G"]), max(tail_range[1], g_valid["log10_p_G"], gB_valid["log10_p_G"]))

# ---- chunk: swing ----
p_ips_A = pct_ips_A / 100
# population a fraction t of the way from the actual result towards Ipsos's picture
def swing_log10p(t, k=None, p_from=None, p_to=None):
    k = k_ips_A.values if k is None else k
    p_from = p_rik_A.values if p_from is None else np.asarray(p_from)
    p_to = p_ips_A.values if p_to is None else np.asarray(p_to)
    p2 = (1 - t) * p_from + t * p_to
    return g_test(k, p2)["log10_p_G"]
t_05 = brentq(lambda t: swing_log10p(t) - np.log10(0.05), 0, 1)
t_50 = brentq(lambda t: swing_log10p(t) - np.log10(0.50), 0, 1)
to_sl   = lambda t: t * (p_rik_A.iloc[0] - p_ips_A.iloc[0]) * rik_meta["ballots"]     # voters who must move TO Studentska
from_us = lambda t: t * (p_ips_A.iloc[1] - p_rik_A.iloc[1]) * rik_meta["ballots"]     # ... of whom from Ujedinjena Srbija
voters_05, voters_50 = to_sl(t_05), to_sl(t_50)

swing_df = pd.DataFrame({"t": np.arange(0, 1.0001, 0.01)})
swing_df["voters"] = to_sl(swing_df["t"])
swing_df["log10_p"] = [swing_log10p(t) for t in swing_df["t"]]

# how much does the answer depend on the assumed sample size?
voters_05_by_n = np.array([to_sl(brentq(lambda t, k=pct_to_counts(pct_ips_A.values, n, 3): swing_log10p(t, k=k) - np.log10(0.05), 0, 1))
                           for n in n_grid])

# ---- chunk: fig-swing ----
g = (ggplot(swing_df, aes(x="voters", y="log10_p")) +
     geom_hline(yintercept=np.log10(0.05), linetype="dotted", color="#898781") +
     annotate("text", x=0, y=np.log10(0.05) + 0.8, label="p = 0.05 (“could be chance”)", ha="left", va="bottom", size=10, color="#52514e") +
     geom_line(color="#2a78d6", size=1.1) +
     geom_vline(xintercept=voters_05, color="#eb6834", size=0.7) +
     annotate("label", x=voters_05 * 0.985, y=swing_df["log10_p"].min() * 0.5,
              label=f"about {fmt_n(round(voters_05, -3))} voters\nwould have had to move\nto Studentska lista",
              ha="right", size=10, color="#eb6834", fill="#fcfcfb", label_size=0) +
     scale_x_continuous(labels=comma, expand=(0.02, 0, 0.08, 0)) +
     scale_y_continuous(labels=pow10, expand=(0.02, 0, 0.1, 0)) +
     labs(x="Voters assumed to have moved to Studentska lista between 22 and 25 October",
          y="Probability of Ipsos's result", title="Scenario 1: how big a late swing would rescue the poll?") +
     theme_report(figure_size=(9, 4.5)))
fig(g, "Figure 8. How the probability of Ipsos's result rises as we assume ever more voters moved to Studentska lista after the poll (the electorate drifting from Ipsos's picture towards the actual result).")

# ---- chunk: turnout ----
registered_processed = rik_meta["turned_out"] / (rik_meta["turnout_pct"] / 100)        # registered voters at the processed stations
extra_voters = (rik_meta["turnout_pct"] - ipsos_turnout["est"]) / 100 * registered_processed
votes_t = scheme_a["votes"].values.astype(float).copy()
votes_t[0] = votes_t[0] - extra_voters                       # all extra voters assumed to be Studentska voters
p_turnout = votes_t / votes_t.sum()
g_turn = g_test(k_ips_A.values, p_turnout)

kable_report(pd.DataFrame({"Category": grp_levels, "RIK % (actual)": [fmt_pct(100 * x) for x in p_rik_A],
                           "% if all extra turnout were Studentska voters": [fmt_pct(100 * x) for x in p_turnout],
                           "Ipsos %": [fmt_pct(x, 1) for x in pct_ips_A]}),
             caption=f"Table 15. Scenario 2: the electorate Ipsos expected, under the extreme assumption that all {fmt_n(round(extra_voters, -3))} voters it did not expect voted Studentska lista.")

# ---- chunk: combined ----
t_05_comb = brentq(lambda t: swing_log10p(t, p_from=p_turnout) - np.log10(0.05), 0, 1)
voters_05_comb = t_05_comb * (p_turnout[0] - p_ips_A.iloc[0]) * votes_t.sum()

# ---- chunk: undecided ----
rows = []
for u in [0.20, 0.25, 0.30, 0.40, 0.50]:
    q = (p_rik_A.values - (1 - u) * p_ips_A.values) / u          # how the undecided must have voted
    rows.append({"u": u, "Studentska lista": 100 * q[0], "Ujedinjena Srbija": 100 * q[1], "Manjinske liste": 100 * q[2], "Other": 100 * q[3]})
undec = pd.DataFrame(rows)
kable_report(pd.DataFrame({"Undecided share": undec["u"].map(lambda u: fmt_pct(100 * u, 0)),
                           **{c: undec[c].map(lambda x: f"{x:.0f}%".replace("-", "−")) for c in undec.columns[1:]}}),
             caption="Table 16. Scenario 4: how the undecided voters would have had to split for Ipsos's decided picture to have been exactly right. Negative values are impossible.")

# ---- chunk: words ----
# the headline exponent in words (short scale, as used in English; long scale, as used in Serbian)
def pow_words(l10):
    # the number 1/p written as  lead × 10^e, then named: e.g. 3 × 10^36 = "three undecillion"
    e = math.floor(-l10); lead = round(10 ** (-l10 - e))
    if lead == 10: lead, e = 1, e + 1
    grp, mult = e // 3, lead * 10 ** (e % 3)                        # mult in 1..900
    short_names = {10: "nonillion", 11: "decillion", 12: "undecillion", 13: "duodecillion"}
    # Serbian forms: singular (1), paucal (2–4), plural (5+)
    long_names = {10: ("kvintilion", "kvintiliona", "kvintiliona"), 11: ("kvintilijarda", "kvintilijarde", "kvintilijardi"),
                  12: ("sekstilion", "sekstiliona", "sekstiliona"), 13: ("sekstilijarda", "sekstilijarde", "sekstilijardi")}
    if grp not in short_names: return f"one in a {lead} followed by {e} zeros"
    en_small = ["one", "two", "three", "four", "five", "six", "seven", "eight", "nine"]
    en_tens  = ["ten", "twenty", "thirty", "forty", "fifty", "sixty", "seventy", "eighty", "ninety"]
    en_mult  = en_small[mult - 1] if mult < 10 else (en_tens[mult // 10 - 1] if mult < 100 else f"{en_small[mult // 100 - 1]} hundred")
    sr_small = ["", "dva", "tri", "četiri", "pet", "šest", "sedam", "osam", "devet"]
    sr_tens  = ["deset", "dvadeset", "trideset", "četrdeset", "pedeset", "šezdeset", "sedamdeset", "osamdeset", "devedeset"]
    sr_hund  = ["sto", "dvesta", "trista", "četiristo", "petsto", "šeststo", "sedamsto", "osamsto", "devetsto"]
    sr_mult  = sr_small[mult - 1] if mult < 10 else (sr_tens[mult // 10 - 1] if mult < 100 else sr_hund[mult // 100 - 1])
    sr_form  = 0 if mult == 1 else (1 if (mult % 10 in (2, 3, 4) and mult % 100 not in (12, 13, 14)) else 2)
    sr_name  = long_names[grp][sr_form]
    return (f"**one in {en_mult} {short_names[grp]}** (a {lead} followed by {e} zeros; in the long-scale naming used in Serbian "
            f"and most of continental Europe, *{(sr_mult + ' ' + sr_name).strip()}*)")

# ---- chunk: session-info ----
import platform, scipy, matplotlib, plotnine
print(f"Python {platform.python_version()} on {platform.platform()}")
for m in (np, scipy, pd, matplotlib, plotnine):
    print(f"{m.__name__:<12}{m.__version__}")

Python 3.14.6 on macOS-26.6.2-arm64-arm-64bit-Mach-O
numpy       2.5.4
scipy       1.18.1
pandas      3.1.0
matplotlib  3.11.2
plotnine    0.15.8

Category	Ipsos count (out of n)	Expected under H₀	Difference	In standard deviations
Studentska lista	844	1,119.0	-275	-12.0
Ujedinjena Srbija	1,000	822.2	+178	+7.9
Manjinske liste	49	45.5	+4	+0.5
Other	217	123.3	+94	+8.7

Category	Ipsos 95% interval	Actual result	Outside the interval by (pts)	Ipsos's own SE (pts)	Distance in SEs	One-sided probability
Studentska lista	37.9 – 42.1	53.03%	10.9	1.07	12.2	2 × 10⁻³⁴ (about 1 in a 4 followed by 33 zeros)
Ujedinjena Srbija	45.3 – 49.5	38.97%	6.3	1.07	7.9	2 × 10⁻¹⁵ (about 1 in a 6 followed by 14 zeros)

Assumed n	G test (Scheme A)	Studentska this low	Ujedinjena this high
1,500	7 × 10⁻²⁶ (about 1 in a 1 followed by 25 zeros)	3 × 10⁻²⁴ (about 1 in a 3 followed by 23 zeros)	2 × 10⁻¹¹ (about 1 in a 5 followed by 10 zeros)
2,110	3 × 10⁻³⁶ (about 1 in a 3 followed by 35 zeros)	2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros)	2 × 10⁻¹⁵ (about 1 in a 4 followed by 14 zeros)
2,900	9 × 10⁻⁵⁰ (about 1 in a 1 followed by 49 zeros)	4 × 10⁻⁴⁵ (about 1 in a 3 followed by 44 zeros)	2 × 10⁻²⁰ (about 1 in a 6 followed by 19 zeros)
3,624	3 × 10⁻⁶² (about 1 in a 4 followed by 61 zeros)	8 × 10⁻⁵⁶ (about 1 in a 1 followed by 55 zeros)	4 × 10⁻²⁵ (about 1 in a 3 followed by 24 zeros)

Undecided share	Studentska lista	Ujedinjena Srbija	Manjinske liste	Other
20%	105%	5%	2%	−12%
25%	92%	14%	2%	−8%
30%	83%	19%	2%	−5%
40%	73%	26%	2%	−1%
50%	66%	31%	2%	1%

No.	Electoral list	Group	Ipsos %	95% interval, low	95% interval, high
1	Aleksandar Vučić – Ujedinjena Srbija	Ujedinjena Srbija	47.4	45.3	49.5
2	Ivica Dačić – Faktor stabilnosti!	Other	3.2	2.5	3.9
3	Studentska lista – Studenti pobeđuju	Studentska lista	40.0	37.9	42.1
4	Rasim Ljajić – Ljudski (SDP – CDP)	Manjinske liste	0.5	0.2	0.8
5	Savez vojvođanskih Mađara – dr Balint Pastor	Manjinske liste	1.2	0.7	1.7
6	Usame Zukorlić (SPP) Pokret za ujedinjenje	Other	0.3	0.1	0.5
7	Ruska stranka – Srbija u BRIKSU	Other	0.2	0.0	0.4
8	Autentična desnica – dr Miloš Jovanović	Other	1.7	1.1	2.3
9	Evropska Srbija – Marinika Tepić – Zdravko Ponoš	Other	0.8	0.4	1.2
10	SDA Sandžaka – Čuvari Sandžaka	Manjinske liste	0.3	0.1	0.5
11	Bogoljub Karić – Pokret Snaga Srbije – BK	Other	0.6	0.3	0.9
12	Liberalno zelena stranka	Other	0.1	0.0	0.2
13	Biram pravdu – Biram budućnost – Koalicija Snaga	Other	0.1	0.0	0.2
14	Stranka rusko-srpskog jedinstva	Other	0.2	0.0	0.4
15	Ćale ovo je za tebe – Petar Đurić	Other	0.2	0.0	0.4
16	Vladimir Gajić – Alternativa za Srbiju	Other	0.3	0.1	0.5
17	Albanska koalicija – Koalicioni Shqiptar	Manjinske liste	0.3	0.1	0.5
18	Majka Srbija – Gruda	Other	0.1	0.0	0.2
	(assumed invalid ballots)	Other	2.5

No.	Electoral list	Group	Votes	% of ballots	Seats
1	Aleksandar Vučić – Ujedinjena Srbija	Ujedinjena Srbija	1,726,554	38.97%	103
2	Ivica Dačić – Faktor stabilnosti!	Other	94,571	2.13%	0
3	Studentska lista – Studenti pobeđuju	Studentska lista	2,349,777	53.03%	141
4	Rasim Ljajić – Ljudski (SDP – CDP)	Manjinske liste	20,415	0.46%	1
5	Savez vojvođanskih Mađara – dr Balint Pastor	Manjinske liste	35,004	0.79%	2
6	Usame Zukorlić (SPP) Pokret za ujedinjenje	Other	12,132	0.27%	0
7	Ruska stranka – Srbija u BRIKSU	Other	6,492	0.15%	0
8	Autentična desnica – dr Miloš Jovanović	Other	39,485	0.89%	0
9	Evropska Srbija – Marinika Tepić – Zdravko Ponoš	Other	12,562	0.28%	0
10	SDA Sandžaka – Čuvari Sandžaka	Manjinske liste	12,381	0.28%	1
11	Bogoljub Karić – Pokret Snaga Srbije – BK	Other	20,031	0.45%	0
12	Liberalno zelena stranka	Other	2,683	0.06%	0
13	Biram pravdu – Biram budućnost – Koalicija Snaga	Other	1,186	0.03%	0
14	Stranka rusko-srpskog jedinstva	Other	4,284	0.10%	0
15	Ćale ovo je za tebe – Petar Đurić	Other	2,018	0.05%	0
16	Vladimir Gajić – Alternativa za Srbiju	Other	3,704	0.08%	0
17	Albanska koalicija – Koalicioni Shqiptar	Manjinske liste	27,681	0.62%	2
18	Majka Srbija – Gruda	Other	3,227	0.07%	0
	(invalid ballots)	Other	56,472	1.27%	0

Quantity	Value
Registered voters	6,451,292
Voters who turned out	4,433,332
Turnout	68.98% (RIK figure, relative to registered voters at processed stations)
Ballots in boxes (valid + invalid)	4,430,659 (the 2,673 voters who were handed a ballot but did not deposit one are the difference from the row above)
Valid ballots	4,374,187
Invalid ballots	56,472 (1.27%)
Polling stations processed	8,339 of 8,392

Electoral list	Ipsos %	95% interval	Half-width (pts)	Implied n_eff
Ujedinjena Srbija	47.4%	45.3 – 49.5	2.1	2,172
Studentska lista	40.0%	37.9 – 42.1	2.1	2,091
Faktor stabilnosti (SPS)	3.2%	2.5 – 3.9	0.7	2,429
Autentična desnica	1.7%	1.1 – 2.3	0.6	1,783
SVM	1.2%	0.7 – 1.7	0.5	1,822
Evropska Srbija	0.8%	0.4 – 1.2	0.4	1,905
Karić – PSS	0.6%	0.3 – 0.9	0.3	2,546
Ljajić – Ljudski	0.5%	0.2 – 0.8	0.3	2,124

Category	Ipsos count	Expected	z-score	Question	Probability
Studentska lista	844	1,119.0	−12.0	Share this low or lower?	2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros)
Ujedinjena Srbija	1,000	822.2	+7.9	Share this high or higher?	2 × 10⁻¹⁵ (about 1 in a 4 followed by 14 zeros)

Approach	P(Studentska ≤ 40.0% and Ujedinjena ≥ 47.4%)
exact algorithm (Levin's representation; pmultinom in R)	2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros)
Monte Carlo, 2,000,000 simulated polls	0 of 2,000,000 (so below 5 × 10⁻⁷ (about 1 in a 2 followed by 6 zeros))

Scheme	G statistic	Degrees of freedom	p-value (G)	p-value (Pearson χ², less reliable this far out)
A: four categories	168	3	3 × 10⁻³⁶ (about 1 in a 3 followed by 35 zeros)	3 × 10⁻³⁸ (about 1 in a 3 followed by 37 zeros)
B: all lists (13 categories after merging 7 tiny ones)	188	12	1 × 10⁻³³ (about 1 in a 9 followed by 32 zeros)	3 × 10⁻³⁷ (about 1 in a 3 followed by 36 zeros)

Method	Probability under random sampling
Method 1 – Point probability of Ipsos's exact vector (for reference; not a tail)	1 × 10⁻⁴¹ (about 1 in a 1 followed by 41 zeros)
Method 2a – Binomial tail: Studentska lista this low or lower (exact)	2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros)
Method 2b – Binomial tail: Ujedinjena Srbija this high or higher (exact)	2 × 10⁻¹⁵ (about 1 in a 4 followed by 14 zeros)
Method 3 – Joint tail: both at once (exact; plain Monte Carlo)	exact: 2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros); Monte Carlo: 0 of 2,000,000
Method 4a – Exact multinomial test, exhaustive enumeration (exact; the headline)	2 × 10⁻³⁶ (about 1 in a 5 followed by 35 zeros)
Method 4b – Exact multinomial test, plain Monte Carlo	below 5 × 10⁻⁷ (about 1 in a 2 followed by 6 zeros) (0 of 2,000,000)
Method 4c – Exact multinomial test, importance sampling (biased low here; see text)	2 × 10⁻³⁷ (about 1 in a 7 followed by 36 zeros) (10 runs: exponents -37.0 to -36.6)
Method 5a – Likelihood-ratio G test, Scheme A (chi-square approximation)	3 × 10⁻³⁶ (about 1 in a 3 followed by 35 zeros)
Method 5b – Likelihood-ratio G test, Scheme B, all lists (chi-square approximation)	1 × 10⁻³³ (about 1 in a 9 followed by 32 zeros)

Category	RIK % (all voters)	RIK % (Serbia proper, without Kosovo i Metohija and abroad)	Ipsos %
Studentska lista	53.03%	52.94%	40.0%
Ujedinjena Srbija	38.97%	38.96%	47.4%
Manjinske liste	2.16%	2.21%	2.3%
Other	5.84%	5.89%	10.3%

Test	All ballots (main analysis)	Valid votes only
Studentska lista this low or lower (binomial)	2 × 10⁻³³ (about 1 in a 4 followed by 32 zeros)	1 × 10⁻³¹ (about 1 in a 8 followed by 30 zeros)
G test, Scheme A	3 × 10⁻³⁶ (about 1 in a 3 followed by 35 zeros)	8 × 10⁻³³ (about 1 in a 1 followed by 32 zeros)
G test, Scheme B	1 × 10⁻³³ (about 1 in a 9 followed by 32 zeros)	1 × 10⁻³⁰ (about 1 in a 7 followed by 29 zeros)

Category	RIK % (actual)	% if all extra turnout were Studentska voters	Ipsos %
Studentska lista	53.03%	49.93%	40.0%
Ujedinjena Srbija	38.97%	41.55%	47.4%
Manjinske liste	2.16%	2.30%	2.3%
Other	5.84%	6.23%	10.3%

1 Introduction¶

1.1 Background¶

1.2 What this document does¶

1.3 What this analysis can and cannot tell you¶

2 Setup¶

3 The data¶

3.1 Why we need both sources on the record¶

3.2 The Ipsos poll (published 22 October 2026)¶

3.3 The official results (RIK, snapshot of 31 October 2026, 11:00)¶

3.4 Poll versus reality, list by list¶

4 Making the two sources comparable¶

4.1 Why this step matters¶

4.2 Scheme A: four categories¶

4.3 Scheme B: all nineteen categories¶

5 The multinomial distribution¶

5.1 Why we need a probability model¶

Concept: the multinomial distribution

5.2 Seeing it: what random polls from the real electorate look like¶

6 How big was Ipsos's sample, really?¶

6.1 Why this is the crux¶

Concept: confidence intervals, standard error and effective sample size

6.2 Backing out the effective sample size from Ipsos's own intervals¶

7 The probability of Ipsos's result¶

7.1 Setting up the hypothesis¶

7.2 Why "the probability of exactly this result" is the wrong question¶

Concept: point probability versus tail probability

7.3 Method 1: the point probability, for reference only¶

7.4 Method 2: one list at a time (binomial tails)¶

Concept: the binomial tail and the z-score

7.5 Method 3: both lists at once (the multinomial cumulative distribution)¶

7.6 Method 4: the exact multinomial test¶

Concept: "at least as unusual" measured by likelihood

7.7 Method 5: the likelihood-ratio (G) test, which can reach the extreme tail¶

Concept: the G statistic and the chi-square distribution

7.8 All methods side by side¶

7.9 The same conclusion from Ipsos's own stated precision¶

7.10 The sampling cloud¶

8 Sensitivity: does the conclusion depend on our choices?¶

8.1 Why check¶

8.2 Varying the sample size¶

8.3 Varying the population: excluding Kosovo i Metohija and the diaspora¶

8.4 Varying the base: valid votes only¶

9 What would it take to rescue the poll?¶

9.1 Why ask¶

9.2 Scenario 1: a late swing between 22 and 25 October¶

9.3 Scenario 2: differential turnout (the net surplus)¶

9.4 Scenario 3: both together¶

9.5 Scenario 4: undecided voters¶

10 Summary and conclusions¶

10.1 What we did¶

10.2 What we found¶

10.3 What it means, and what it does not¶

10.4 Caveats¶

11 Appendix¶

11.1 Sources¶

11.2 References¶

11.3 Reproducibility¶