#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ drift-versus-selection.py Science Journaling Club, Volume 1 Issue 1, Fall 2024, "Populations and Chance". THE QUESTION ------------ A new mutation that raises its carrier's fitness by a few percent almost always disappears anyway. We want to know how the probability that it survives depends on the selection coefficient s and on the population size N, and how well the two classical approximations to that probability actually hold. WHAT THIS PROGRAM IS -------------------- This is a computation, not an observation. The club has no laboratory, no flies, no bacterial cultures and no sequencing machine. Nothing here was measured in the physical world. Every number below is produced by pseudo-random sampling from a Wright-Fisher Markov chain that we wrote ourselves, seeded once at the top so that anybody re-running the file gets exactly the same figures. The "experiment" is the Monte Carlo run. Where we write "measured" we mean "estimated from our own simulated replicates". THE MODEL --------- A population holds a constant 2N gene copies at one locus (N diploid individuals, or 2N haploids, the arithmetic is the same). Two alleles: wild type, and a mutant whose relative fitness as a gene copy is 1 + s. One generation is: 1. selection, deterministic, acting on gene-copy frequency p: p' = p (1 + s) / (1 + s p) 2. reproduction, the entire stochastic part: i_{t+1} ~ Binomial(2N, p') The chain starts at i = 1, a single new copy, and runs until it hits an absorbing boundary: 0 (lost) or 2N (fixed). We repeat this many times with independent random numbers and count how often it reaches 2N. That count divided by the number of replicates is our estimate of the fixation probability u. The diffusion limit of this chain has infinitesimal mean M(p) = s p (1 - p) and variance V(p) = p (1 - p) / (2N), which is what makes Kimura's formula below the right one to compare against. WHAT WE COMPARE AGAINST ----------------------- Haldane (1927), branching-process limit, large N, small s: u ~= 2s Haldane's underlying branching process itself, without the small-s step: pi solves pi = 1 - exp(-(1+s) pi) (a single copy leaves Poisson(1+s) descendants while it is rare) Kimura (1962), diffusion result at starting frequency p = 1/(2N): u = (1 - exp(-2s)) / (1 - exp(-4 N s)) Neutral case, exact for the Wright-Fisher chain by the martingale/symmetry argument, no approximation involved: u = 1 / (2N) Mean time to fixation, conditional on fixing: neutral, Kimura & Ohta (1969), starting from a single copy: ~ 4N strong selection, deterministic logistic sweep: ~ (2/s) ln(2N) ASSUMPTIONS, STATED PLAINLY --------------------------- * Constant population size. No growth, no crashes, no bottlenecks. * Non-overlapping generations. Everybody reproduces at once and then dies. * Random mating, one panmictic pool. No geography, no population structure. * One locus, two alleles, no recombination, no linked sites. Nothing else in the genome is under selection, so there is no background selection and no genetic draft. * Selection is genic and constant: fitness 1 + s per mutant copy, so in diploid terms the heterozygote is exactly intermediate. No dominance, no overdominance, no epistasis, no frequency dependence. * No recurrent mutation. The mutant arises once and is never created again, and never mutates back. * The environment does not change, so s does not change. * N is the census size and also the effective size. In real populations these differ, usually with Ne well below N. LIMITATIONS ----------- Every one of those assumptions is false of some real population, and several are false of most. A real beneficial mutation in a real species competes with other sweeps, sits on a chromosome that carries deleterious load, experiences a variable environment, and lives in a population whose size fluctuates. The numbers here describe the idealised chain, and the honest claim we can make is about that chain. The comparison with the classical formulas is the part that generalises: the formulas were derived for this same idealisation, so if they fail here they fail for reasons of mathematics rather than biology. RUNTIME ------- About one and a half minutes on a laptop. Requires numpy. """ from __future__ import annotations import math import sys import time import numpy as np # --------------------------------------------------------------------------- # The seed. Stated in the article. Change it and every number below changes # slightly, within the quoted Monte Carlo error. # --------------------------------------------------------------------------- MASTER_SEED = 20240917 N_VALUES = [50, 200, 1000, 5000, 10000] S_VALUES = [0.0, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1] MIN_REPS = 50_000 # floor; the study plan asks for at least 20,000 MAX_REPS = 8_000_000 # ceiling so the whole grid finishes in minutes TARGET_REL_SE = 0.10 # standard error under 10 percent of the estimate SAFETY = 3.0 # replicate multiplier, so a cell that happens to # fix fewer times than predicted still clears the # 10% target rather than landing just outside it # cells whose full replicate-by-replicate record we keep for the convergence plot CONVERGENCE_CELLS = [(1000, 0.01), (1000, 0.0), (10000, 0.002)] # --------------------------------------------------------------------------- # analytic reference values # --------------------------------------------------------------------------- def kimura(N, s): """Kimura's diffusion fixation probability for one new copy out of 2N.""" if s == 0.0: return 1.0 / (2 * N) num = -math.expm1(-2.0 * s) # 1 - exp(-2s), accurate for tiny s den = -math.expm1(-4.0 * N * s) # 1 - exp(-4Ns) return num / den def haldane(s): """Haldane's branching-process approximation, after the small-s step.""" return 2.0 * s def haldane_branching(s): """ The escape probability of Haldane's branching process itself, before the small-s linearisation. While the mutant is rare a single copy leaves a Poisson(1+s) number of copies next generation, and the probability that its lineage never dies out solves pi = 1 - exp(-(1+s) pi). Haldane then expanded this for small s and got 2s. Solved here by bisection. """ if s <= 0.0: return 0.0 lo, hi = 1e-12, 1.0 for _ in range(200): mid = 0.5 * (lo + hi) if 1.0 - math.exp(-(1.0 + s) * mid) - mid > 0.0: lo = mid else: hi = mid return 0.5 * (lo + hi) def wilson_interval(k, n, z=1.959963985): """Wilson score interval. Behaves sensibly when k is small, which Wald does not.""" if n == 0: return (0.0, 0.0) p = k / n d = 1.0 + z * z / n centre = (p + z * z / (2 * n)) / d half = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d return (max(0.0, centre - half), min(1.0, centre + half)) def reps_for(N, s): """Choose the replicate count so the relative standard error clears the target.""" u = kimura(N, s) needed = SAFETY * (1.0 - u) / (u * TARGET_REL_SE ** 2) return int(min(MAX_REPS, max(MIN_REPS, math.ceil(needed)))) # --------------------------------------------------------------------------- # the Wright-Fisher chain # --------------------------------------------------------------------------- def run_cell(N, s, n_reps, rng, keep_record=False): """ Run n_reps independent Wright-Fisher trajectories from a single mutant copy. All replicates advance together: `counts` holds the mutant copy number in every replicate still segregating, and one call to rng.binomial steps the whole surviving set forward one generation. Replicates drop out as they hit 0 or 2N, so the vector shrinks and the late generations are cheap. """ two_n = 2 * N counts = np.ones(n_reps, dtype=np.int32) alive = np.arange(n_reps, dtype=np.int64) fixed = np.zeros(n_reps, dtype=bool) fix_time = np.zeros(n_reps, dtype=np.int32) gen = 0 gen_cap = 60 * N + 2000 # safety valve; we report if it ever bites while alive.size and gen < gen_cap: gen += 1 p = counts[alive] / two_n if s != 0.0: p = p * (1.0 + s) / (1.0 + s * p) new = rng.binomial(two_n, p) counts[alive] = new done = (new == 0) | (new == two_n) if done.any(): just = alive[done] won = just[counts[just] == two_n] if won.size: fixed[won] = True fix_time[won] = gen alive = alive[~done] return { "N": N, "s": s, "reps": n_reps, "n_fixed": int(fixed.sum()), "fix_times": fix_time[fixed], "max_gen": gen, "unresolved": int(alive.size), "record": fixed if keep_record else None, } # --------------------------------------------------------------------------- # reporting helpers # --------------------------------------------------------------------------- def rule(ch="-", width=118): return ch * width def banner(title): print() print(rule("=")) print(title) print(rule("=")) def reldiff(measured, predicted): if predicted == 0.0: return float("nan") return (measured - predicted) / predicted def main(): t_start = time.time() print(rule("=")) print("HOW OFTEN DOES A GOOD GENE LOSE?") print("Wright-Fisher simulation of fixation probability under genic selection") print("Science Journaling Club, Volume 1 Issue 1, Fall 2024") print(rule("=")) print() print("This output is generated by simulation. No organism, population or") print("laboratory measurement is involved anywhere in this file.") print() print("numpy version : " + np.__version__) print("python version : " + sys.version.split()[0]) print("master seed : %d" % MASTER_SEED) print("generator : numpy PCG64, one independent stream per grid cell") print("population sizes N : " + str(N_VALUES)) print("selection coefficients s : " + str(S_VALUES)) print("replicate floor : {:,}".format(MIN_REPS)) print("replicate ceiling : {:,}".format(MAX_REPS)) print("target relative SE : {:.0%}".format(TARGET_REL_SE)) print("replicate safety factor : {:.1f}x the count the target alone implies".format(SAFETY)) print() print("Model, one generation: p' = p(1+s)/(1+sp) then i' ~ Binomial(2N, p').") print("Start i = 1. Absorb at i = 0 (lost) or i = 2N (fixed).") # one independent, reproducible stream per cell seeds = np.random.SeedSequence(MASTER_SEED).spawn(len(N_VALUES) * len(S_VALUES)) results = {} records = {} cell_no = 0 banner("SECTION 1. GRID OF RUNS") print() print("{:>7} {:>8} {:>10} {:>12} {:>11} {:>12} {:>12} {:>8} {:>10} {:>7}".format( "N", "s", "4Ns", "replicates", "fixations", "u_hat", "SE", "rel SE", "gens run", "sec")) print(rule()) for N in N_VALUES: for s in S_VALUES: rng = np.random.default_rng(seeds[cell_no]) cell_no += 1 n_reps = reps_for(N, s) keep = (N, s) in CONVERGENCE_CELLS t0 = time.time() res = run_cell(N, s, n_reps, rng, keep_record=keep) res["seconds"] = time.time() - t0 k, n = res["n_fixed"], res["reps"] u_hat = k / n se = math.sqrt(u_hat * (1 - u_hat) / n) res["u_hat"] = u_hat res["se"] = se res["ci"] = wilson_interval(k, n) results[(N, s)] = res if keep: records[(N, s)] = res["record"] res["record"] = None rel_se = se / u_hat if u_hat > 0 else float("nan") print("{:>7} {:>8.3f} {:>10.1f} {:>12,} {:>11,} {:>12.6f} {:>12.6f} {:>7.1%} {:>10,} {:>7.1f}".format( N, s, 4 * N * s, n, k, u_hat, se, rel_se, res["max_gen"], res["seconds"])) if res["unresolved"]: print(" WARNING: %d replicates hit the generation cap" % res["unresolved"]) sys.stdout.flush() # ----------------------------------------------------------------------- banner("SECTION 2. NEUTRAL CASE AGAINST THE EXACT ANSWER 1/(2N)") print() print("For s = 0 the Wright-Fisher chain gives u = 1/(2N) exactly, because the") print("mutant copy number is a martingale and every one of the 2N copies in the") print("founding generation is equally likely to be the ancestor of the whole") print("future population. This is a check on the code, not on the theory.") print() print("{:>7} {:>12} {:>11} {:>12} {:>12} {:>26} {:>12} {:>10} {:>7}".format( "N", "replicates", "fixations", "u_hat", "SE", "95% CI (Wilson)", "1/(2N)", "rel diff", "z")) print(rule()) neutral_rows = [] for N in N_VALUES: r = results[(N, 0.0)] exact = 1.0 / (2 * N) rd = reldiff(r["u_hat"], exact) se_null = math.sqrt(exact * (1 - exact) / r["reps"]) z = (r["u_hat"] - exact) / se_null lo, hi = r["ci"] inside = "yes" if lo <= exact <= hi else "NO" neutral_rows.append((N, r["reps"], r["n_fixed"], r["u_hat"], r["se"], lo, hi, exact, rd, z, inside)) ci_txt = "[{:.6f}, {:.6f}]".format(lo, hi) print("{:>7} {:>12,} {:>11,} {:>12.6f} {:>12.6f} {:>26} {:>12.6f} {:>+9.2%} {:>+7.2f}".format( N, r["reps"], r["n_fixed"], r["u_hat"], r["se"], ci_txt, exact, rd, z)) print() print("Does the 95% interval contain the exact value?") for row in neutral_rows: print(" N = {:>5}: {}".format(row[0], row[10])) worst = max(abs(r[9]) for r in neutral_rows) print() print("Largest |z| across the five neutral cells: {:.2f}".format(worst)) print("Five independent z scores from a correct simulator should mostly sit") print("inside +/-2. Anything past 3 would mean the code is wrong.") # ----------------------------------------------------------------------- banner("SECTION 3. FULL GRID AGAINST HALDANE AND KIMURA") print() print("u_hat our estimate from the replicates") print("2s Haldane 1927, the large-N branching-process limit") print("Kimura (1 - exp(-2s)) / (1 - exp(-4Ns)), Kimura 1962") print("d_Hal (u_hat - 2s) / 2s") print("d_Kim (u_hat - Kimura) / Kimura") print("z_Kim (u_hat - Kimura) / SE, how many standard errors out we are") print() print("{:>7} {:>7} {:>9} {:>10} {:>8} {:>11} {:>10} {:>10} {:>10} {:>11} {:>9} {:>8}".format( "N", "s", "4Ns", "reps", "fix", "u_hat", "SE", "2s", "d_Hal", "Kimura", "d_Kim", "z_Kim")) print(rule()) grid_rows = [] for N in N_VALUES: for s in S_VALUES: r = results[(N, s)] kim = kimura(N, s) hal = haldane(s) d_kim = reldiff(r["u_hat"], kim) d_hal = reldiff(r["u_hat"], hal) if s > 0 else float("nan") z_kim = (r["u_hat"] - kim) / r["se"] if r["se"] > 0 else float("nan") grid_rows.append((N, s, 4 * N * s, r["reps"], r["n_fixed"], r["u_hat"], r["se"], hal, d_hal, kim, d_kim, z_kim)) hal_txt = "{:>10}".format("--") if s == 0 else "{:>10.6f}".format(hal) dhal_txt = "{:>10}".format("--") if s == 0 else "{:>+9.1%}".format(d_hal) print("{:>7} {:>7.3f} {:>9.1f} {:>10,} {:>8,} {:>11.6f} {:>10.6f} {} {} {:>11.6f} {:>+8.1%} {:>+8.2f}".format( N, s, 4 * N * s, r["reps"], r["n_fixed"], r["u_hat"], r["se"], hal_txt, dhal_txt, kim, d_kim, z_kim)) print(rule(".")) # ----------------------------------------------------------------------- banner("SECTION 4. WHERE THE APPROXIMATIONS BREAK DOWN") print() print("Haldane's 2s assumes the mutant never feels the ceiling: an infinite") print("population in which a lineage either dies out early or escapes for good.") print("Kimura's diffusion keeps the ceiling and so keeps N.") print() print("Haldane's error, sorted by 4Ns, positive s only:") print() print("{:>10} {:>7} {:>7} {:>11} {:>10} {:>10} {:>11} {:>9}".format( "4Ns", "N", "s", "u_hat", "2s", "d_Hal", "Kimura", "d_Kim")) print(rule()) positive = [row for row in grid_rows if row[1] > 0] for row in sorted(positive, key=lambda r: r[2]): print("{:>10.1f} {:>7} {:>7.3f} {:>11.6f} {:>10.6f} {:>+9.1%} {:>11.6f} {:>+8.1%}".format( row[2], row[0], row[1], row[5], row[7], row[8], row[9], row[10])) print() print("Summary by 4Ns band (mean absolute relative error across the cells in the band):") print() bands = [(0.0, 1.0), (1.0, 4.0), (4.0, 20.0), (20.0, 100.0), (100.0, 1e9)] print("{:>18} {:>7} {:>14} {:>14}".format("4Ns band", "cells", "mean |d_Hal|", "mean |d_Kim|")) print(rule()) for lo, hi in bands: sel = [r for r in positive if lo <= r[2] < hi] if not sel: continue mh = sum(abs(r[8]) for r in sel) / len(sel) mk = sum(abs(r[10]) for r in sel) / len(sel) label = "[{:g}, {:g})".format(lo, hi) if hi < 1e8 else "[{:g}, inf)".format(lo) print("{:>18} {:>7} {:>13.1%} {:>13.1%}".format(label, len(sel), mh, mk)) print() print("Worst Haldane failures (largest |d_Hal|):") for row in sorted(positive, key=lambda r: -abs(r[8]))[:6]: print(" N={:>6} s={:.3f} 4Ns={:>7.1f} u_hat={:.6f} 2s={:.6f} d_Hal={:+.1%}".format( row[0], row[1], row[2], row[5], row[7], row[8])) print() print("Worst Kimura disagreements (largest |z_Kim|):") for row in sorted(positive, key=lambda r: -abs(r[11]))[:6]: print(" N={:>6} s={:.3f} 4Ns={:>7.1f} u_hat={:.6f} Kimura={:.6f} d_Kim={:+.1%} z={:+.2f}".format( row[0], row[1], row[2], row[5], row[9], row[10], row[11])) zs = [row[11] for row in grid_rows] print() print("All {} cells, z against Kimura: mean {:+.3f}, max |z| {:.2f}".format( len(zs), sum(zs) / len(zs), max(abs(z) for z in zs))) print("Cells with |z| > 2: {} of {} (about 2 expected by chance at 5%)".format( sum(1 for z in zs if abs(z) > 2), len(zs))) # ----------------------------------------------------------------------- banner("SECTION 4B. IS THE KIMURA GAP AT LARGE s REAL OR IS IT NOISE?") print() print("Each column of the grid shares a value of s across five values of N, and") print("the five runs used five independent random streams. Combining the five z") print("scores by Stouffer's method (sum of z divided by sqrt(5)) turns five weak") print("hints into one sharp test of whether the deviation from Kimura is real.") print() print("{:>8} {:>28} {:>14} {:>12} {:>12}".format( "s", "mean d_Kim across the five N", "sum of z", "pooled z", "verdict")) print(rule()) for sv in S_VALUES: col = [row for row in grid_rows if row[1] == sv] md = sum(r[10] for r in col) / len(col) sz = sum(r[11] for r in col) pz = sz / math.sqrt(len(col)) verdict = "real" if abs(pz) > 3 else ("marginal" if abs(pz) > 2 else "noise") print("{:>8.3f} {:>27.2%} {:>14.2f} {:>12.2f} {:>12}".format(sv, md, sz, pz, verdict)) print() print("Now the same measurements against three different large-N predictions.") print("Only cells with 4Ns >= 100 are used, where finite N is no longer the") print("limiting factor and every prediction should apply.") print() print("2s Haldane's published approximation") print("1 - exp(-2s) Kimura's formula in the large-N limit") print("branching pi Haldane's own branching process, unlinearised:") print(" pi = 1 - exp(-(1+s) pi)") print() print("{:>8} {:>7} {:>12} {:>11} {:>9} {:>13} {:>9} {:>13} {:>9}".format( "s", "cells", "u_hat pooled", "2s", "d", "1-exp(-2s)", "d", "branching pi", "d")) print(rule()) for sv in S_VALUES: if sv == 0.0: continue col = [row for row in grid_rows if row[1] == sv and row[2] >= 100.0] if not col: continue tot_fix = sum(r[4] for r in col) tot_rep = sum(r[3] for r in col) up = tot_fix / tot_rep hal = 2.0 * sv kinf = -math.expm1(-2.0 * sv) pib = haldane_branching(sv) print("{:>8.3f} {:>7} {:>12.6f} {:>11.6f} {:>+8.1%} {:>13.6f} {:>+8.1%} {:>13.6f} {:>+8.1%}".format( sv, len(col), up, hal, reldiff(up, hal), kinf, reldiff(up, kinf), pib, reldiff(up, pib))) print() print("Pooled u_hat is total fixations over total replicates for the cells in") print("that row, which is the right way to combine binomial counts.") banner("SECTION 5. TIME TO FIXATION, CONDITIONAL ON FIXING") print() print("Counted in generations, only over the replicates that actually fixed.") print("Neutral benchmark: Kimura & Ohta 1969 give a mean of about 4N generations") print("for a single neutral copy. Strong-selection benchmark: the deterministic") print("logistic sweep from 1/(2N) to 1 - 1/(2N) takes about (2/s) ln(2N).") print() print("{:>7} {:>7} {:>9} {:>9} {:>10} {:>8} {:>8} {:>8} {:>8} {:>11} {:>8}".format( "N", "s", "4Ns", "n fixed", "mean t", "SE", "median", "p10", "p90", "benchmark", "ratio")) print(rule()) for N in N_VALUES: for s in S_VALUES: r = results[(N, s)] ft = r["fix_times"] if ft.size < 20: print("{:>7} {:>7.3f} {:>9.1f} {:>9,} too few fixations to summarise".format( N, s, 4 * N * s, ft.size)) continue mean_t = float(ft.mean()) se_t = float(ft.std(ddof=1) / math.sqrt(ft.size)) med = float(np.median(ft)) p10 = float(np.percentile(ft, 10)) p90 = float(np.percentile(ft, 90)) if s == 0.0: bench = 4.0 * N blabel = "{:>11.0f}".format(bench) elif 4 * N * s >= 20: bench = (2.0 / s) * math.log(2 * N) blabel = "{:>11.0f}".format(bench) else: bench = float("nan") blabel = "{:>11}".format("--") ratio = mean_t / bench if bench == bench else float("nan") rlabel = "{:>8.2f}".format(ratio) if ratio == ratio else "{:>8}".format("--") print("{:>7} {:>7.3f} {:>9.1f} {:>9,} {:>10.1f} {:>8.1f} {:>8.0f} {:>8.0f} {:>8.0f} {} {}".format( N, s, 4 * N * s, ft.size, mean_t, se_t, med, p10, p90, blabel, rlabel)) print(rule(".")) print() print("Neutral fixation times against the 4N prediction:") for N in N_VALUES: r = results[(N, 0.0)] ft = r["fix_times"] if ft.size < 20: continue mean_t = float(ft.mean()) se_t = float(ft.std(ddof=1) / math.sqrt(ft.size)) print(" N = {:>5}: measured {:>9.1f} +/- {:>6.1f} generations, predicted {:>7.0f}, rel diff {:+.1%}".format( N, mean_t, se_t, 4 * N, reldiff(mean_t, 4.0 * N))) # ----------------------------------------------------------------------- banner("SECTION 6. CONVERGENCE OF THE ESTIMATE AS TRIALS ACCUMULATE") print() print("The running estimate after the first k replicates, with the Monte Carlo") print("standard error at that k. Replicates are independent and identically") print("distributed, so the prefix of the record is itself a valid smaller run.") for (N, s) in CONVERGENCE_CELLS: rec = records.get((N, s)) if rec is None: continue kim = kimura(N, s) csum = np.cumsum(rec) n_tot = rec.size print() print("Cell N = {}, s = {:g}. Kimura value {:.6f}. {:,} replicates.".format( N, s, kim, n_tot)) print() print("{:>12} {:>10} {:>12} {:>12} {:>8} {:>14}".format( "trials k", "fixations", "u_hat(k)", "SE(k)", "rel SE", "u_hat/Kimura")) print(rule()) ks = [] k = 500 while k < n_tot: ks.append(k) nxt = int(k * 1.5) k = nxt if nxt > k else k + 1 ks.append(n_tot) for k in ks: c = int(csum[k - 1]) u = c / k se = math.sqrt(u * (1 - u) / k) if u > 0 else float("nan") rel = se / u if u > 0 else float("nan") if c == 0: print("{:>12,} {:>10,} {:>12.6f} {:>12} {:>8} {:>14}".format( k, c, 0.0, "--", "--", "--")) else: print("{:>12,} {:>10,} {:>12.6f} {:>12.6f} {:>7.1%} {:>14.4f}".format( k, c, u, se, rel, u / kim)) # ----------------------------------------------------------------------- banner("SECTION 7. THE HEADLINE NUMBERS") print() for N in [1000, 10000]: for s in [0.01, 0.05]: r = results[(N, s)] lost = 1 - r["u_hat"] print("N = {:>6}, s = {:.2f}: a single new copy fixes {:.3f}% of the time, so it is lost {:.3f}% of the time (1 in {:.0f} survives).".format( N, s, r["u_hat"] * 100, lost * 100, 1 / r["u_hat"])) print() r = results[(1000, 0.01)] print("A mutation worth a 1% fitness gain in a population of 1000 diploids") print("fixes {:.3f}% of the time in our runs. Haldane says {:.1f}%. Kimura says {:.3f}%.".format( r["u_hat"] * 100, 2 * 0.01 * 100, kimura(1000, 0.01) * 100)) print() r0 = results[(1000, 0.0)] print("The same mutation with no advantage at all fixes {:.4f}% of the time (exact answer {:.4f}%).".format( r0["u_hat"] * 100, 100 / (2 * 1000))) print("So a 1% advantage multiplies the survival odds by {:.1f}, and the mutation still loses {:.2f} times out of 100.".format( r["u_hat"] / r0["u_hat"], (1 - r["u_hat"]) * 100)) banner("SECTION 8. PRECISION AUDIT") print() print("The study plan asks for a standard error under 10% of the estimate in") print("every cell. Realised relative standard errors, worst first:") print() audit = [] for (N, s), r in results.items(): rel = r["se"] / r["u_hat"] if r["u_hat"] > 0 else float("inf") audit.append((rel, N, s, r["reps"], r["n_fixed"])) audit.sort(reverse=True) print("{:>8} {:>7} {:>12} {:>10} {:>9}".format("rel SE", "N", "s", "reps", "fixations")) print(rule()) for rel, N, s, reps, k in audit[:8]: print("{:>7.2%} {:>7} {:>12.3f} {:>10,} {:>9,}".format(rel, N, s, reps, k)) print() worst_rel = audit[0][0] print("Worst relative standard error anywhere in the grid: {:.2%}".format(worst_rel)) print("Target met in all {} cells: {}".format( len(audit), "yes" if worst_rel < 0.10 else "NO")) total_reps = sum(r["reps"] for r in results.values()) total_fix = sum(r["n_fixed"] for r in results.values()) print() print(rule("=")) print("Total replicates simulated : {:,}".format(total_reps)) print("Total fixations observed : {:,}".format(total_fix)) print("Wall clock : {:.1f} s".format(time.time() - t_start)) print(rule("=")) if __name__ == "__main__": main()