Fields

Prelude.Fields collects the classic applied families of stochastic models: spin systems on a lattice, stochastic epidemics, population genetics, bond percolation and the Cramér–Lundberg insurance model. Each family gives you the one-step update to iterate from C# together with the closed-form quantities the subject is known for (the final size of an epidemic, the fixation probability of a mutant, the Lundberg bound on ruin).

The lattice here is a one-dimensional ring of $n$ sites: each site’s neighbours are one step either way, and the last site neighbours the first. Keep in mind that the one-dimensional Ising model has no phase transition at positive temperature; onsager states the two-dimensional critical temperature for reference, and this ring does not reach it.

Randomness is never hidden. Each stochastic step takes its draws as arguments — uniforms in $[0, 1)$ called u and v, or standard normals called z — and you generate them in C#. Spin states are vectors of $\pm 1$; epidemic states are the vector $[S, I, R]$ of counts. Comparisons return 1 or 0, written below with the Iverson bracket $[\,\cdot\,]$.

Notation. Positions in a vector are counted from zero, $0, 1, \dots, n-1$, and on the ring indices are taken modulo $n$, so $s_{-1} = s_{n-1}$ and $s_n = s_0$. Divisors are guarded with $\operatorname{nonzero}(d) = d + 10^{-300}\,[\,\lvert d \rvert < 10^{-300}\,]$ (see Chains), which equals $d$ everywhere except within $10^{-300}$ of zero.

open Prelude.Fields

Lattices

The Ising model on a ring with coupling $J$ and external field $h$ has energy

\[E(s) = -J \sum_i s_i s_{i+1} - h \sum_i s_i .\]

Glauber dynamics picks one site and flips it with probability $1/(1 + e^{\Delta E / T})$, a logistic in the energy change that needs no separate acceptance ratio. At zero field the model is symmetric and the magnetisation wanders around zero.

aroundOf

aroundOf(s) · s spins vector n → neighbour sum at each site vector n

Let $s = (s_0, \dots, s_{n-1}) \in \mathbb{R}^n$ be the spins on a ring of $n \ge 3$ sites, and let $A = \operatorname{ringAhead}(s)$. Then, for $i = 0, \dots, n-1$,

\[\operatorname{aroundOf}(s)_i = \sum_{j=0}^{n-1} \big(A + A^{\top}\big)_{ij}\, s_j = s_{i-1} + s_{i+1} \quad (\text{indices mod } n) .\]

The sum of each site’s two ring neighbours. (On rings of one or two sites the two neighbours coincide and are counted twice.)

Intuition The local field each site feels from its neighbours: $+2$ when both point up, $-2$ when both point down, 0 when they disagree. Everything a nearest-neighbour model needs to know about the surroundings is in this one number per site.

Use it when you compute local energies, flip costs or neighbour majority rules on a ring — flipCost and siteEnergy are built on it — or smooth a periodic 1-D signal by summing neighbours. It builds an n×n neighbour mask, so memory grows as $n^2$; keep rings to thousands of sites, not millions.

neighbour-interactionslattice-simulation

def aroundOf(s) ≔ einsum("ij,j->i", ringAhead(s) + ringBehind(s), s)

flipChance

flipChance(s, coupling, field, temperature) · s spins vector n · coupling J scalar · field h scalar · temperature T scalar → Glauber flip probability of each site vector n

Let $s \in \{-1, 1\}^n$ be the spins, $J \in \mathbb{R}$ (the argument coupling), $h \in \mathbb{R}$ (the argument field) and $T > 0$ (the argument temperature). Let $\Delta E_i = \operatorname{flipCost}(s, J, h)_i = 2\big(J s_i (s_{i-1} + s_{i+1}) + h s_i\big)$. Then, for $i = 0, \dots, n-1$,

\[\operatorname{flipChance}(s, J, h, T)_i = \frac{1}{1 + \exp\!\big(\Delta E_i / \operatorname{nonzero}(T)\big)} .\]

The probability of accepting a flip at each site, where $\Delta E_i$ is the energy change the flip would cause.

Intuition A soft decision to flip: flips that lower the energy are likely, flips that raise it are unlikely, and temperature sets how soft the cut-off is. A flip that changes nothing has probability exactly 1/2; as $T \to 0$ only energy-lowering flips happen, and at high $T$ every flip is a coin toss.

Use it when you want the flip probability of every site at once — to drive flipped, to update all sites of one sublattice in parallel, or to see which sites are “frustrated” (probabilities well above 1/2). Temperature is in units of $J$ (with $k_B = 1$). Unlike a Metropolis rule it is never exactly 1, so dynamics are slightly slower but smooth in $\Delta E$.

lattice-simulationaccept-reject-decisions

def flipChance(s, coupling, field, temperature) ≔ 1 / (1 + exp(flipCost(s, coupling, field) / nonzero(temperature)))

flipped

flipped(s, coupling, field, temperature, u, v) · s spins vector n · coupling J scalar · field h scalar · temperature T scalar · u uniform draw choosing the site scalar · v uniform draw deciding the flip scalar → spins after one Glauber step vector n

Let $s \in \{-1, 1\}^n$, $J, h \in \mathbb{R}$ and $T > 0$ as above, and let $u, v \in [0, 1)$ be independent uniform draws. Let $k = \lfloor u n \rfloor$, $c = \operatorname{chosenSite}(s, u) = e_k$ (the one-hot vector at position $k$) and $q = \operatorname{flipChance}(s, J, h, T)_k$. Then, for $i = 0, \dots, n-1$,

\[\operatorname{flipped}(s, J, h, T, u, v)_i = s_i \,\big(1 - 2\, c_i\, [\,v < q\,]\big) .\]

One Glauber update: site $k$ is chosen, and its spin is flipped if $v$ is below that site’s flip probability. Iterate it with two fresh uniforms per step.

Intuition One tick of a heat bath: pick a random site, let it reconsider its spin given its neighbours and the temperature, and leave everything else alone. Repeated many times, the configurations it visits follow the Boltzmann distribution at that temperature.

Use it when you simulate an Ising ring — equilibrium sampling, relaxation after a quench, domain coarsening, or teaching MCMC on a physical model. Measure time in sweeps of $n$ steps, discard the first several sweeps as burn-in, and draw u and v fresh each step. Remember the ring never orders at $T > 0$ (see ordersAtAll), so long runs at low temperature show large domains that still wander.

lattice-simulationsimulationmagnetic-ordering

def flipped(s, coupling, field, temperature, u, v) ≔ s · (1 − 2 · chosenSite(s, u) · (v < einsum("i,i->", flipChance(s, coupling, field, temperature), chosenSite(s, u))))

magnetisation

magnetisation(s) · s spins vector n → mean spin scalar

Let $s \in \{-1, 1\}^n$ be the spins. Then

\[\operatorname{magnetisation}(s) = \frac{1}{n} \sum_{i=0}^{n-1} s_i .\]

The magnetisation per site.

Intuition How much the spins agree overall: $+1$ all up, $-1$ all down, near 0 for a disordered mix. It is the order parameter — the single number that distinguishes an ordered phase from a disordered one.

Use it when you track a spin simulation over time, measure the response to an external field, or compute susceptibility from its fluctuations. At zero field it averages to 0 by symmetry even in an ordered system that flips as a whole, so average $$

$or$m^2$instead. It is also the mean opinion in voter or consensus models coded as$\pm 1$$.

magnetic-orderinglattice-simulation

def magnetisation(s) ≔ mean(s)

totalEnergy

totalEnergy(s, coupling, field) · s spins vector n · coupling J scalar · field h scalar → energy scalar

Let $s \in \{-1, 1\}^n$ be the spins on a ring of $n \ge 3$ sites, and $J, h \in \mathbb{R}$ (the arguments coupling and field). Then

\[\operatorname{totalEnergy}(s, J, h) = \frac12 \sum_{i=0}^{n-1} \operatorname{siteEnergy}(s, J, h)_i = -J \sum_{i=0}^{n-1} s_i s_{i+1} - \frac{h}{2} \sum_{i=0}^{n-1} s_i .\]

The energy of a configuration, computed as half the sum of the per-site energies so that each bond is counted once. Note that the halving also applies to the field term, so the value equals the energy $E(s)$ above only when $h = 0$.

Intuition Each agreeing pair of neighbours lowers the energy by $J$ and each disagreeing pair raises it by $J$, so the energy counts domain walls: $-nJ$ for a uniform ring, plus $2J$ for every boundary between up and down regions.

Use it when you monitor equilibration (the energy settles before most other observables), compute the heat capacity from energy fluctuations, or compare configurations. Because the field term is halved, use it with $h = 0$, or add $-\tfrac{h}{2} \sum_i s_i$ yourself to get the full energy when a field is present.

energy-evaluationlattice-simulation

def totalEnergy(s, coupling, field) ≔ einsum("i->", siteEnergy(s, coupling, field)) / 2

onsager

onsager → critical temperature scalar

The constant

\[\operatorname{onsager} = 2.269185314213022 \approx \frac{2}{\ln(1 + \sqrt{2})} .\]

Onsager’s critical temperature of the two-dimensional square-lattice Ising model (with $J = k_B = 1$). Given for reference: the ring in this module has no transition.

Intuition On a square grid, below about 2.27 (in units of $J$) spins spontaneously line up into a magnet; above it they stay disordered. It is the most famous exact result in statistical mechanics and the benchmark every Ising simulation is checked against.

Use it when you need a reference temperature: to label a plot, to scale a temperature sweep for a 2-D code written elsewhere, or to test such a code’s estimated critical point. Do not expect the ring in this module to show anything at this temperature; it has no transition at all.

reference-valuesmagnetic-ordering

def onsager       ≔ 2.269185314213022

ordersAtAll

ordersAtAll → 0

The constant

\[\operatorname{ordersAtAll} = 0 .\]

States as a value that the one-dimensional ring does not order at any positive temperature.

Intuition In one dimension a single broken bond costs a fixed $2J$ but can sit anywhere, so the entropy gained always wins at any $T > 0$: domain walls appear and destroy long-range order. The critical temperature of a ring is therefore zero.

Use it when you want the ring’s critical temperature as a value beside onsager — in a table or a test that checks a simulation’s long-run $$ m $$ decays toward 0 on large rings. For a real ordering transition you need two or more dimensions.

reference-valuesmagnetic-ordering

def ordersAtAll   ≔ 0

pottsEnergy

pottsEnergy(s, coupling) · s site states (integers 0 … q−1) vector n · coupling J scalar → energy scalar

Let $s \in \{0, 1, \dots, q-1\}^n$ be the site states on a ring of $n \ge 3$ sites and $J \in \mathbb{R}$ (the argument coupling). Then

\[\operatorname{pottsEnergy}(s, J) = -\frac{J}{2} \sum_{i=0}^{n-1} \operatorname{agreements}(s)_i = -J \sum_{i=0}^{n-1} [\,s_i = s_{i+1}\,] .\]

The energy of a Potts configuration on the ring: each pair of neighbours in the same state lowers the energy by $J$, and any two different states cost nothing. With $q = 2$ it is the Ising model in different clothes. The module does not provide a Potts order parameter.

Intuition Neighbours are rewarded only for matching, whatever the labels are; with more states there are more ways to disagree, so order is harder to establish. It is the natural energy for any labelling where “same as neighbour” is good and all differences are equally bad.

Use it when you score or simulate $q$-state labellings on a ring — cell types, opinions with several options, segmenting a periodic 1-D signal into classes with a smoothness prior. States must be integers $0 \dots q-1$. The module has no Potts update step, so drive a Metropolis step in C# with differences of this energy.

energy-evaluationlattice-simulationneighbour-interactions

def pottsEnergy(s, coupling) ≔ 0 − coupling · einsum("i->", agreements(s)) / 2

Building blocks

ringAhead(s) — Let $s \in \mathbb{R}^n$, with positions $i, j = 0, \dots, n-1$. Then $\operatorname{ringAhead}(s)_{ij} = [\,j - i = 1\,] + [\,j - i = 1 - n\,]$, which is 1 exactly when $j = i + 1 \bmod n$. An n×n 0/1 matrix marking each site’s next neighbour on the ring.
ringBehind(s) — Let $s \in \mathbb{R}^n$ and $A = \operatorname{ringAhead}(s)$. Then $\operatorname{ringBehind}(s) = A^{\top}$, with entry $(i, j)$ equal to 1 exactly when $j = i - 1 \bmod n$. Each site’s previous neighbour.
siteEnergy(s, coupling, field) — Let $s \in \{-1, 1\}^n$ and $J, h \in \mathbb{R}$. Then $\operatorname{siteEnergy}(s, J, h)_i = -J\, s_i\, (s_{i-1} + s_{i+1}) - h\, s_i$. The energy attributed to each site, counting both of its bonds and its field term.
flipCost(s, coupling, field) — Let $s \in \{-1, 1\}^n$ and $J, h \in \mathbb{R}$. Then $\operatorname{flipCost}(s, J, h)_i = \Delta E_i = 2\big(J\, s_i\, (s_{i-1} + s_{i+1}) + h\, s_i\big)$. The energy change of flipping each site; used by flipChance.
chosenSite(s, u) — Let $s \in \mathbb{R}^n$ and $u \in [0, 1)$. Then $\operatorname{chosenSite}(s, u)_i = [\,i = \lfloor u n \rfloor\,]$ for $i = 0, \dots, n-1$, the basis vector $e_{\lfloor u n \rfloor}$. A one-hot vector marking the chosen site.
agreements(s) — Let $s \in \mathbb{R}^n$ and $A = \operatorname{ringAhead}(s)$. Then $\operatorname{agreements}(s)_i = \sum_{j} (A + A^{\top})_{ij}\, [\,s_j = s_i\,] = [\,s_{i-1} = s_i\,] + [\,s_{i+1} = s_i\,]$. For each site, how many of its two neighbours share its state; used by pottsEnergy.

Epidemics

The stochastic SIR model tracks $[S, I, R]$ in a population of size $N$. Infections occur at rate $\beta S I / N$ and recoveries at rate $\gamma I$. The step below is the jump chain of the Gillespie algorithm: it decides which event happens next (you draw the waiting time yourself if you need it). Once $I = 0$ the state is absorbing and the step leaves it unchanged.

The basic reproduction number $R_0 = \beta/\gamma$ decides everything. Below one the outbreak dies out. Above one it either dies out early by chance — with probability $R_0^{-k}$ from $k$ initial infected — or grows to the final size $f$ that solves

\[1 - f = e^{-R_0 f}.\]

infecting

infecting(x, beta, size) · x state [S, I, R] vector 3 · beta transmission rate scalar · size population N scalar → infection rate scalar

Let $x = (S, I, R) \in \mathbb{R}^3$ (the argument x, positions $0, 1, 2$) be the state, $\beta \ge 0$ the transmission rate and $N > 0$ (the argument size) the population. Then

\[\operatorname{infecting}(x, \beta, N) = \frac{\beta\, S\, I}{\operatorname{nonzero}(N)} .\]

The current infection rate.

Intuition Mass action: new infections happen when susceptibles meet infected people, so the rate grows with both counts, and dividing by $N$ turns $I$ into the chance that a given contact is infectious. $\beta$ is the number of infectious contacts one person makes per unit time.

Use it when you build the Gillespie simulation of an SIR model (it is half the total event rate), write the right-hand side of the deterministic SIR equations, or check how fast an outbreak is growing right now. Pair it with recovering; epidemicStep uses both.

epidemic-modellingsimulation

def infecting(x, beta, size)  ≔ beta · elementAt(x, 0) · elementAt(x, 1) / nonzero(size)

recovering

recovering(x, gamma) · x state [S, I, R] vector 3 · gamma recovery rate scalar → recovery rate scalar

Let $x = (S, I, R) \in \mathbb{R}^3$ be the state and $\gamma \ge 0$ the recovery rate. Then

\[\operatorname{recovering}(x, \gamma) = \gamma\, I .\]

The current recovery rate.

Intuition Each infected person recovers independently at rate $\gamma$, so the infectious period is exponential with mean $1/\gamma$ and the total rate is simply proportional to how many are sick.

Use it when you need the recovery half of the total event rate — for the Gillespie waiting time, which is exponential with rate infecting + recovering — or the outflow term of the SIR equations. Set $\gamma$ as one over the mean infectious period in your time unit (a 10-day period gives 0.1 per day).

epidemic-modellingsimulation

def recovering(x, gamma)      ≔ gamma · elementAt(x, 1)

epidemicStep

epidemicStep(x, beta, gamma, size, u) · x state [S, I, R] vector 3 · beta transmission rate scalar · gamma recovery rate scalar · size population N scalar · u uniform draw scalar → next state vector 3

Let $x = (S, I, R) \in \mathbb{R}^3$ be the state, $\beta, \gamma \ge 0$ the rates, $N > 0$ (the argument size) the population, and let $u \in [0, 1)$ be a uniform draw. Let $a = \operatorname{infecting}(x, \beta, N)$, $r = \operatorname{recovering}(x, \gamma)$ and $q = a / \operatorname{nonzero}(a + r)$, and let $e_0, e_1, e_2$ be the basis vectors of $\mathbb{R}^3$. Then

\[\operatorname{epidemicStep}(x, \beta, \gamma, N, u) = x + [\,a + r > 0\,] \Big( [\,u < q\,]\,(e_1 - e_0) + [\,u \ge q\,]\,(e_2 - e_1) \Big) .\]

One event of the epidemic. If $u$ is below the share of the total rate due to infection, one susceptible becomes infected; otherwise one infected recovers. With no events possible (in particular with no infected) the state is returned unchanged, so counts never go negative.

Intuition Two competing clocks, infection and recovery; whichever rings first decides the next event, and its chance of ringing first is its share of the total rate. Following one event at a time gives the exact randomness of a small outbreak — including early chance extinction — that the smooth SIR equations miss.

Use it when you simulate outbreaks in small or moderate populations, estimate the probability of a major outbreak, or study the spread of final sizes. Call it with a fresh uniform per event; an epidemic ends within at most $2N$ events, so that many calls always suffice. If you need calendar time, add an exponential waiting time with rate infecting + recovering yourself. For the large-population answer use finalSize.

epidemic-modellingsimulationoutbreak-risk

def epidemicStep(x, beta, gamma, size, u) ≔ x + spreading(x, beta, gamma, size) · ((u < infectionShare(x, beta, gamma, size)) · ((positions(x) = 1) − (positions(x) = 0)) + (u ≥ infectionShare(x, beta, gamma, size)) · ((positions(x) = 2) − (positions(x) = 1)))

reproduction

reproduction(beta, gamma) · beta transmission rate scalar · gamma recovery rate scalar → R₀ scalar

Let $\beta \ge 0$ be the transmission rate and $\gamma > 0$ the recovery rate. Then

\[\operatorname{reproduction}(\beta, \gamma) = R_0 = \frac{\beta}{\operatorname{nonzero}(\gamma)} .\]

The basic reproduction number.

Intuition How many people one case infects in a fully susceptible population: transmissions per unit time ($\beta$) times the time spent infectious ($1/\gamma$). Above 1 each generation of cases is larger than the last; below 1 it shrinks.

Use it when you decide whether an outbreak can take off, compare pathogens or interventions (an intervention that cuts contacts by a fraction $c$ multiplies $R_0$ by $1 - c$), or estimate the immune share needed to stop spread, $1 - 1/R_0$. It is the input to diesOut and finalSize.

epidemic-modellingoutbreak-riskintervention-planning

def reproduction(beta, gamma)  ≔ beta / nonzero(gamma)

diesOut

diesOut(beta, gamma, start) · beta transmission rate scalar · gamma recovery rate scalar · start initial number infected k scalar → probability of early extinction scalar

Let $\beta, \gamma > 0$ be the rates, $R_0 = \operatorname{reproduction}(\beta, \gamma)$, and $k \in \mathbb{N}$ (the argument start) the initial number infected. Then

\[\operatorname{diesOut}(\beta, \gamma, k) = \min\big(R_0^{-k},\ 1\big) .\]

The probability that the outbreak dies out by chance — the branching-process approximation, which is 1 whenever $R_0 \le 1$.

Intuition Early on, each case either infects someone before recovering (odds $\beta : \gamma$) or does not, so a lone case’s chain of transmission is lost with probability $1/R_0$. With $k$ independent starting cases every chain must fail, hence $R_0^{-k}$. Even a very transmissible disease often fizzles when introduced once.

Use it when you estimate the risk that an imported case starts an outbreak ($1 - \text{diesOut}$), plan how many introductions are needed before an outbreak is likely, or check a stochastic simulation’s share of minor outbreaks. It assumes a large population and exponential infectious periods; it ignores depletion of susceptibles.

outbreak-riskepidemic-modelling

def diesOut(beta, gamma, start) ≔ min(reproduction(beta, gamma) ^ (0 − start), 1)

finalSize

finalSize(beta, gamma) · beta transmission rate scalar · gamma recovery rate scalar → final fraction infected scalar

Let $\beta, \gamma > 0$ be the rates and $R_0 = \operatorname{reproduction}(\beta, \gamma)$. Let $f_0 = 0.5$ and $f_{k+1} = 1 - e^{-R_0 f_k}$ for $k = 0, 1, \dots, 199$. Then

\[\operatorname{finalSize}(\beta, \gamma) = [\,R_0 > 1\,]\; f_{200} .\]

The fraction of the population eventually infected in a major outbreak, the positive root of $1 - f = e^{-R_0 f}$, found by 200 fixed-point iterations. It is 0 when $R_0 \le 1$.

Intuition Someone escapes infection only if they escape every case, and the chance of that falls exponentially with the total infectious pressure $R_0 f$. The final size is the self-consistent point where the share never infected, $1 - f$, equals that escape probability. It does not depend on how fast the epidemic runs, only on $R_0$.

Use it when you estimate the attack rate of an unmitigated outbreak, the burden avoided by lowering $R_0$, or check the typical final size of stochastic runs that did not fizzle. It is a large-population, homogeneous-mixing answer given that a major outbreak occurs. Just above $R_0 = 1$ the fixed-point iteration converges slowly, so 200 iterations may undershoot there.

epidemic-modellingattack-rate

def finalSize(beta, gamma) ≔ (reproduction(beta, gamma) > 1) · sizeLoop(200, reproduction(beta, gamma), 0.5)

Building blocks

epidemicPace(x, beta, gamma, size) — Let $x = (S, I, R) \in \mathbb{R}^3$, $\beta, \gamma \ge 0$ and $N > 0$. Then $\operatorname{epidemicPace}(x, \beta, \gamma, N) = \beta S I / \operatorname{nonzero}(N) + \gamma I$. The total event rate, infection plus recovery.
infectionShare(x, beta, gamma, size) — With $x, \beta, \gamma, N$ as above, let $a = \operatorname{infecting}(x, \beta, N)$ and $\rho = \operatorname{epidemicPace}(x, \beta, \gamma, N)$. Then $\operatorname{infectionShare}(x, \beta, \gamma, N) = a / \operatorname{nonzero}(\rho)$. The fraction of the total rate due to infection.
spreading(x, beta, gamma, size) — With $x, \beta, \gamma, N$ as above, $\operatorname{spreading}(x, \beta, \gamma, N) = [\,\operatorname{epidemicPace}(x, \beta, \gamma, N) > 0\,]$. 1 while any event can happen, 0 once the epidemic is over.
sizeLoop(k, r, f) — Let $k \in \mathbb{N}$ and $r, f \in \mathbb{R}$. Then $\operatorname{sizeLoop}(0, r, f) = f$ and $\operatorname{sizeLoop}(k, r, f) = \operatorname{sizeLoop}(k - 1, r, 1 - e^{-r f})$ for $k \ge 1$. The fixed-point iteration $f \leftarrow 1 - e^{-r f}$ repeated $k$ times; used by finalSize.

Genetics

Wright–Fisher resamples a population of $N$ gene copies each generation: an allele at frequency $p$ is first changed by selection and mutation, then drifts. With no selection the frequency is a martingale, yet fixation is certain, because the only absorbing states are 0 and 1.

The Moran model replaces one individual at a time, which makes its fixation probability exact. With selection coefficient $s$ and $r = 1/(1+s)$, a mutant present in $i$ of $N$ copies fixes with probability

\[\frac{1 - r^i}{1 - r^N},\]

the gambler’s-ruin formula again, and $i/N$ in the neutral case $s = 0$.

selected

selected(p, selection) · p allele frequency scalar · selection selection coefficient s scalar → frequency after selection scalar

Let $p \in [0, 1]$ be the allele frequency and $s > -1$ (the argument selection) the selection coefficient. Then

\[\operatorname{selected}(p, s) = \frac{p\,(1 + s)}{\operatorname{nonzero}\big(p\,(1 + s) + 1 - p\big)} .\]

The frequency after one round of selection with relative fitness $1+s$.

Intuition Carriers of the allele leave $1+s$ offspring for every 1 left by the rest, so its share is reweighted by fitness and renormalised. The change is largest at intermediate frequencies and vanishes when the allele is absent or fixed.

Use it when you compute the deterministic effect of selection before drift, predict how fast a beneficial allele spreads in a large population, or build your own generation step (it is the first stage of binomialCount and nextFrequency). $s$ is per generation and must exceed $-1$. It models haploid selection, with no dominance.

allele-frequency-dynamicspopulation-genetics

def selected(p, selection)      ≔ p · (1 + selection) / nonzero(p · (1 + selection) + 1 − p)

mutated

mutated(p, mutation) · p allele frequency scalar · mutation symmetric mutation rate μ scalar → frequency after mutation scalar

Let $p \in [0, 1]$ be the allele frequency and $\mu \in [0, 1]$ (the argument mutation) the mutation rate. Then

\[\operatorname{mutated}(p, \mu) = p\,(1 - \mu) + (1 - p)\,\mu .\]

The frequency after symmetric mutation at rate $\mu$ in both directions.

Intuition Some copies of the allele mutate away and some copies of the other allele mutate into it. With equal rates both ways the frequency is pulled toward 1/2, the balance point where the two flows cancel.

Use it when you add recurrent mutation to a generation step, study mutation–drift balance (mutation reintroduces lost variation), or model any two-state quantity that flips with the same small probability each generation. With $\mu > 0$ the frequencies 0 and 1 are no longer absorbing, so fixation is only temporary. Realistic per-locus rates are tiny; large $\mu$ swamps selection.

allele-frequency-dynamicspopulation-genetics

def mutated(p, mutation)        ≔ p · (1 − mutation) + (1 − p) · mutation

binomialCount

binomialCount(p, u, selection, mutation) · p allele frequency scalar · u uniform draws, one per gene copy vector N · selection s scalar · mutation μ scalar → new allele frequency scalar

Let $p \in [0, 1]$, $s > -1$ (the argument selection) and $\mu \in [0, 1]$ (the argument mutation), let $u \in [0, 1)^N$ be independent uniform draws, and let $p'' = \operatorname{mutated}(\operatorname{selected}(p, s), \mu)$. Then

\[\operatorname{binomialCount}(p, u, s, \mu) = \frac{1}{N} \sum_{i=1}^{N} [\,u_i < p''\,] .\]

One exact Wright–Fisher generation: a binomial resample divided by $N$. It needs $N$ fresh uniforms every generation; reusing the same draws repeats the same generation.

Intuition The next generation’s $N$ gene copies each pick a parent at random from a pool shaped by selection and mutation; counting how many picked the allele is a binomial draw. Chance alone moves the frequency, more strongly in small populations.

Use it when you want exact Wright–Fisher drift — small populations, frequencies near 0 or 1, or checking the diffusion approximation. It costs $N$ uniforms per generation, so for large $N$ or long runs use nextFrequency, which needs one normal.

allele-frequency-dynamicspopulation-geneticssimulation

def binomialCount(p, u, selection, mutation) ≔ einsum("i->", u < mutated(selected(p, selection), mutation)) / count(u)

nextFrequency

nextFrequency(p, z, selection, mutation, population) · p allele frequency scalar · z a standard normal draw scalar · selection s scalar · mutation μ scalar · population N scalar → new allele frequency scalar

Let $p \in [0, 1]$, $s > -1$ (the argument selection), $\mu \in [0, 1]$ (the argument mutation) and $N > 0$ (the argument population), let $z \in \mathbb{R}$ be a standard normal draw, and let $p'' = \operatorname{mutated}(\operatorname{selected}(p, s), \mu)$. Then

\[\operatorname{nextFrequency}(p, z, s, \mu, N) = \min\!\Big(\max\!\Big(p'' + \sqrt{\max\big(p(1 - p)/\operatorname{nonzero}(N),\ 0\big)}\; z,\ 0\Big),\ 1\Big) .\]

One Wright–Fisher generation in the diffusion approximation. The drift scale uses the frequency $p$ before selection and mutation. It needs one normal per generation instead of $N$ uniforms, matches Wright–Fisher’s variance and, with no mutation, keeps 0 and 1 absorbing, so fixation still happens.

Intuition For large $N$ a binomial count is close to normal, so one generation is the deterministic change plus Gaussian noise with the binomial’s variance $p(1-p)/N$. The noise vanishes at 0 and 1, and the result is clipped to $[0, 1]$, so without mutation an allele that is lost or fixed stays that way.

Use it when you simulate drift over many generations or large populations cheaply: fixation times, loss of variation, many replicate populations in parallel. It is least accurate when $Np$ or $N(1-p)$ is small, where binomialCount is exact. With mutation $\mu > 0$ the boundaries are not absorbing, since mutated moves 0 and 1 inward.

allele-frequency-dynamicspopulation-geneticssimulation

def nextFrequency(p, z, selection, mutation, population) ≔ min(max(mutated(selected(p, selection), mutation) + driftScale(p, population) · z, 0), 1)

fixationChance

fixationChance(p) · p allele frequency scalar → fixation probability scalar

Let $p \in [0, 1]$ be the frequency of a neutral allele. Then

\[\operatorname{fixationChance}(p) = p .\]

The probability that a neutral allele eventually fixes is its frequency.

Intuition Eventually every gene copy descends from a single ancestor, and with no selection each current copy is equally likely to be that ancestor. The allele wins if the ancestor carries it, which happens with probability equal to its current share. A new mutant, one copy in $N$, fixes with probability $1/N$.

Use it when you need the neutral baseline — to test whether an allele’s success needs selection, to validate a drift simulation (the share of runs ending at 1 should match), or as the $s = 0$ case of moranFixation. Under selection use moranFixation instead.

fixation-probabilitypopulation-genetics

def fixationChance(p)           ≔ p

heterozygosity

heterozygosity(p) · p allele frequency scalar → heterozygosity scalar

Let $p \in [0, 1]$ be the allele frequency. Then

\[\operatorname{heterozygosity}(p) = 2\,p\,(1 - p) .\]

The expected heterozygosity under random mating.

Intuition The chance that two gene copies drawn at random carry different alleles: one of each, in either order. It peaks at 0.5 when both alleles are equally common and is zero once one is lost, so it measures how much variation a population still has.

Use it when you summarise genetic diversity at a locus, track its decay under drift across simulated generations, or compare observed and expected heterozygosity (a deficit points to inbreeding or population structure). For a two-state variable in general it is twice the variance of a Bernoulli indicator.

genetic-diversitypopulation-genetics

def heterozygosity(p)           ≔ 2 · p · (1 − p)

driftTime

driftTime(population) · population N scalar → generations scalar

Let $N > 0$ (the argument population) be the population size. Then

\[\operatorname{driftTime}(N) = 4N .\]

The time scale of neutral genetic drift, in generations.

Intuition Drift is slow in big populations: the time for chance alone to replace variation grows in proportion to population size. $4N$ is the classic figure for a diploid population of $N$ individuals — for example the mean time for a new neutral mutation that does fix to get there.

Use it when you choose how many generations a drift simulation must run, or judge whether drift can matter on a given timescale compared with selection (roughly $1/s$) or mutation ($1/\mu$). Mind the convention: binomialCount and nextFrequency count $N$ gene copies, for which the corresponding scale is about $2N$.

population-geneticstimescale-estimation

def driftTime(population)       ≔ 4 · population

moranFixation

moranFixation(count, size, selection) · count initial mutant copies i scalar · size population N scalar · selection s scalar → fixation probability scalar

Let $i \in \{0, \dots, N\}$ (the argument count) be the initial number of mutant copies, $N \ge 1$ (the argument size) the population and $s > -1$ (the argument selection) the selection coefficient. Let $\nu = [\,\lvert s \rvert < 10^{-12}\,]$, $s' = s + 0.01\,\nu$ and $r = 1/\operatorname{nonzero}(1 + s')$. Then

\[\operatorname{moranFixation}(i, N, s) = \nu\, \frac{i}{\operatorname{nonzero}(N)} + (1 - \nu)\, \frac{1 - r^{i}}{\operatorname{nonzero}(1 - r^{N})} .\]

The exact Moran fixation probability. When $\lvert s \rvert < 10^{-12}$ it returns the neutral value $i/N$ instead of evaluating the formula at its $0/0$ singularity; the shift to $s' = s + 0.01$ only keeps the unused branch finite.

Intuition Mutant count moves up or down one at a time, with up-steps $1+s$ times as likely as down-steps: a gambler with a biased coin who wins when reaching $N$ and loses at 0. Even a clear advantage is often lost while the mutant is rare; a single copy with small $s$ fixes with probability only about $s/(1+s)$ in a large population.

Use it when you ask whether a beneficial or harmful mutation will take over a finite population — evolution of resistance, spread of an innovation in a fixed-size group, invasion analysis in evolutionary game theory. It is exact for the Moran model; use fixationChance for the neutral shorthand. Time in Moran is one birth–death event, about $N$ events per generation.

fixation-probabilitypopulation-geneticsinvasion-analysis

def moranFixation(count, size, selection) ≔ neutralSelection(selection) · neutralFixation(count, size) + (1 − neutralSelection(selection)) · biasedFixation(count, size, clearOf(selection))

Building blocks

driftScale(p, population) — Let $p \in [0, 1]$ and $N > 0$. Then $\operatorname{driftScale}(p, N) = \sqrt{\max\big(p(1 - p)/\operatorname{nonzero}(N),\ 0\big)}$. The standard deviation of one generation of drift.
moranOdds(selection) — Let $s > -1$. Then $\operatorname{moranOdds}(s) = r = 1/\operatorname{nonzero}(1 + s)$. The ratio used by the Moran formula.
neutralFixation(count, size) — Let $i \ge 0$ and $N \ge 1$. Then $\operatorname{neutralFixation}(i, N) = i / \operatorname{nonzero}(N)$. The neutral fixation probability.
neutralSelection(selection) — Let $s \in \mathbb{R}$. Then $\operatorname{neutralSelection}(s) = [\,\lvert s \rvert < 10^{-12}\,]$. 1 when the selection coefficient is zero to within $10^{-12}$.
clearOf(selection) — Let $s \in \mathbb{R}$. Then $\operatorname{clearOf}(s) = s + 0.01\,[\,\lvert s \rvert < 10^{-12}\,]$. Moves a zero selection coefficient to 0.01 so the unused biased branch never divides 0 by 0.
biasedFixation(count, size, selection) — Let $i \ge 0$, $N \ge 1$, $s > -1$ and $r = \operatorname{moranOdds}(s)$. Then $\operatorname{biasedFixation}(i, N, s) = (1 - r^{i}) / \operatorname{nonzero}(1 - r^{N})$. The gambler’s-ruin fixation formula.

Percolation

Bond percolation on the ring: each bond is open when its uniform draw is below $p$. On a ring a cluster is a run of consecutive open bonds.

runLengths

runLengths(u, p) · u one uniform draw per bond vector n · p open probability scalar → running count vector n

Let $u \in [0, 1)^n$ be independent uniform draws, one per bond at positions $0, \dots, n-1$, let $p \in [0, 1]$, and let $o_j = [\,u_j < p\,]$. Then, for $i = 0, \dots, n-1$,

\[\operatorname{runLengths}(u, p)_i = \sum_{j=0}^{n-1} [\,i \ge j\,]\, o_j\, o_j = \sum_{j=0}^{i} o_j .\]

At each position $i$, the number of open bonds among positions $0, \dots, i$. The count does not reset at a closed bond, so it is a cumulative count of open bonds rather than the length of the current run.

Intuition Walk along the ring and keep a tally of open bonds seen so far. On a ring, clusters are runs of consecutive open bonds, and a counter that restarted at each closed bond would give each run’s length; this one does not restart, so it only tells you how many open bonds precede each position.

Use it when you need the running count of open bonds (for example the open share of each prefix, dividing by $i + 1$). Do not use it for cluster sizes: compute run lengths in C#, as the page’s note on largestRun recommends.

connectivity-analysisrunning-totals

def runLengths(u, p)    ≔ einsum("ij,j->i", (held(positions(u)) ≥ laid(positions(u))) · openAt(u, p), openAt(u, p))

largestRun

largestRun(u, p) · u one uniform draw per bond vector n · p open probability scalar → count scalar

Let $u \in [0, 1)^n$ be independent uniform draws, $p \in [0, 1]$, and $o_j = [\,u_j < p\,]$. Then

\[\operatorname{largestRun}(u, p) = \max_{0 \le i \le n-1}\; o_i \sum_{j=0}^{i} o_j = \sum_{j=0}^{n-1} o_j .\]

The largest value of runLengths at an open bond. Because runLengths does not reset at closed bonds, this is the total number of open bonds, not the length of the longest run: for open pattern 1, 1, 0, 1, 1, 1 it returns 5 rather than 3. Compute the longest run in C# if you need it.

Intuition Percolation asks whether open connections form a path across a system; the largest cluster is how close you came. On a ring that cluster would be the longest run of open bonds, which reaches all the way round only when every bond is open — why the 1-D threshold is $p = 1$.

Use it when — as implemented it returns the total number of open bonds, not the longest run. Use it only as that count (it equals $n \cdot$ openShare), and compute the longest run in C# for cluster-size questions such as network robustness or the probability that a chain of links holds.

connectivity-analysisnetwork-robustness

def largestRun(u, p)    ≔ greatest(runLengths(u, p) · openAt(u, p))

openShare

openShare(u, p) · u one uniform draw per bond vector n · p open probability scalar → fraction open scalar

Let $u \in [0, 1)^n$ be independent uniform draws and $p \in [0, 1]$. Then

\[\operatorname{openShare}(u, p) = \frac{1}{n} \sum_{i=0}^{n-1} [\,u_i < p\,] .\]

The fraction of bonds that are open.

Intuition The realised density of open links in one random configuration. It scatters around $p$ with standard deviation $\sqrt{p(1-p)/n}$, so it tells you how typical this particular draw is.

Use it when you check a percolation sample against its nominal $p$, plot cluster statistics against the density actually realised rather than the target one, or estimate a failure rate from simulated links. It says nothing about connectivity: the same share can be one long run or many short ones.

connectivity-analysismonte-carlo-estimation

def openShare(u, p)     ≔ mean(openAt(u, p))

threshold

threshold → 0.5

The constant

\[\operatorname{threshold} = 0.5 .\]

The bond-percolation threshold of the two-dimensional square lattice, exactly $1/2$ (Kesten), given for reference.

Intuition On an infinite square grid where each bond is open with probability $p$, a connected cluster spanning the whole grid appears exactly when $p$ passes 1/2. By self-duality, open bonds on the grid and closed bonds on the dual grid play symmetric roles, which pins the threshold at the midpoint.

Use it when you need the 2-D reference value: scaling a sweep of $p$ for a grid simulation written elsewhere, testing such a simulation’s estimated threshold, or reasoning about the robustness of a mesh-like network. It does not apply to the ring here, whose threshold is 1.

reference-valuesconnectivity-analysis

def threshold           ≔ 0.5

Building blocks

openAt(u, p) — Let $u \in [0, 1)^n$ and $p \in [0, 1]$. Then $\operatorname{openAt}(u, p)_i = [\,u_i < p\,]$. 1 for each bond whose draw is below p, 0 otherwise.

Insurance

The Cramér–Lundberg model: an insurer starts with surplus $u$, earns premium at rate $c$, and pays claims arriving as a Poisson process at rate $\lambda$ with mean size $m$. Ruin is the surplus going negative. A simulation over a finite horizon can only report that ruin has not happened yet; the Lundberg inequality bounds the infinite-horizon ruin probability,

\[\psi(u) \le e^{-R u},\]

where $R$ is the adjustment coefficient. The coefficient used here, $R = 1/m - \lambda/c$, is the exact one for exponentially distributed claims.

safetyLoading

safetyLoading(premium, rate, claimMean) · premium premium rate c scalar · rate claim arrival rate λ scalar · claimMean mean claim size m scalar → loading θ scalar

Let $c > 0$ (the argument premium) be the premium rate, $\lambda > 0$ (the argument rate) the claim arrival rate and $m > 0$ (the argument claimMean) the mean claim size. Then

\[\operatorname{safetyLoading}(c, \lambda, m) = \theta = \frac{c}{\operatorname{nonzero}(\lambda m)} - 1 .\]

How much the premium exceeds the expected claim outflow. Ruin is certain unless $\theta > 0$.

Intuition The insurer’s margin: premiums bring in $c$ per unit time while claims take out $\lambda m$ on average. With $\theta \le 0$ the surplus has no upward drift, so its random swings eventually take it below zero; any positive margin gives it a chance to escape for good.

Use it when you price a policy portfolio (set $c = (1 + \theta)\lambda m$ for a chosen loading), check that a model insurer is viable before simulating it, or compare portfolios on how much cushion their premiums carry. It uses only means, so it says nothing about how risky large claims are; for that use adjustment and lundberg.

premium-pricingruin-probability

def safetyLoading(premium, rate, claimMean) ≔ premium / nonzero(rate · claimMean) − 1

adjustment

adjustment(premium, rate, claimMean) · premium c scalar · rate λ scalar · claimMean m scalar → adjustment coefficient R scalar

Let $c > 0$ (the argument premium), $\lambda > 0$ (the argument rate) and $m > 0$ (the argument claimMean). Then

\[\operatorname{adjustment}(c, \lambda, m) = R = \max\!\left(\frac{1}{\operatorname{nonzero}(m)} - \frac{\lambda}{\operatorname{nonzero}(c)},\ 0\right) .\]

The adjustment coefficient for exponential claims.

Intuition The rate at which ruin risk decays as starting capital grows: each extra unit of surplus multiplies the risk by about $e^{-R}$. It grows with the premium and with smaller claims, and it is clamped to 0 when the loading is not positive, which makes the Lundberg bound the useless value 1.

Use it when you compare how quickly extra capital buys safety across portfolios, or need $R$ for lundberg. The closed form is exact only for exponentially distributed claims; for other claim-size distributions solve $\lambda(M_X(R) - 1) = cR$ numerically, and heavy-tailed claims (no moment-generating function) have no adjustment coefficient at all.

ruin-probabilityrisk-measurement

def adjustment(premium, rate, claimMean) ≔ max(1 / nonzero(claimMean) − rate / nonzero(premium), 0)

lundberg

lundberg(surplus, premium, rate, claimMean) · surplus initial surplus u scalar · premium c scalar · rate λ scalar · claimMean m scalar → upper bound on ruin probability scalar

Let $u \ge 0$ (the argument surplus) be the initial surplus, $c, \lambda, m > 0$ as above, and $R = \operatorname{adjustment}(c, \lambda, m)$. Then

\[\operatorname{lundberg}(u, c, \lambda, m) = e^{-R u} .\]

The Lundberg bound on the probability of ever being ruined.

Intuition Ruin risk falls exponentially with starting capital. To go broke from surplus $u$ the claims must overcome a cushion that keeps growing with premiums, and the probability of such an excursion shrinks by the factor $e^{-R}$ per unit of cushion.

Use it when you size capital reserves for a target ruin probability — solve $u = -\ln(\text{target})/R$ — or check a finite-horizon simulation, whose ruin rate must stay below this bound. It is an upper bound: for exponential claims the exact value is $\frac{1}{1+\theta} e^{-Ru}$, so it overstates the risk by the factor $1 + \theta$. To see ruin paths, simulate with surplusAfter.

ruin-probabilityrisk-measurementcapital-reserving

def lundberg(surplus, premium, rate, claimMean) ≔ exp(0 − adjustment(premium, rate, claimMean) · surplus)

surplusAfter

surplusAfter(surplus, premium, gap, claim) · surplus current surplus scalar · premium premium rate scalar · gap time until the next claim scalar · claim size of that claim scalar → surplus just after the claim scalar

Let $u \in \mathbb{R}$ (the argument surplus) be the current surplus, $c > 0$ (the argument premium) the premium rate, $\Delta t \ge 0$ (the argument gap) the time until the next claim and $X \ge 0$ (the argument claim) its size. Then

\[\operatorname{surplusAfter}(u, c, \Delta t, X) = u + c\,\Delta t - X .\]

One step of a surplus simulation: premium earned over the gap, less the claim.

Intuition Between claims the surplus only rises, so the only moments ruin can happen are just after a claim. Jumping from claim to claim therefore loses nothing: check the sign after each step.

Use it when you simulate surplus paths — estimating finite-horizon ruin probabilities, testing reinsurance or dividend rules, or checking lundberg. Draw the gap as exponential with rate $\lambda$ ($-\ln(1-u)/\lambda$ from a uniform) and the claim from any size distribution, including heavy-tailed ones the closed forms cannot handle. Stop a path at the first negative value, and remember a run that survives the horizon has only not been ruined yet.

ruin-probabilitysimulationrisk-measurement

def surplusAfter(surplus, premium, gap, claim) ≔ surplus + premium · gap − claim

Examples

The examples below live in one file, Lattice.hb, so the class is Lattice.

Energy, magnetisation and a Glauber step on an Ising ring

Four aligned spins have four satisfied bonds, energy −4. Flipping the second spin breaks two bonds, energy 0. A Glauber step with $u = 0.3$ picks site 1; with $v = 0$ the flip is accepted (its probability is $1/(1+e^4) \approx 0.018$), with $v = 0.5$ it is not.

open Prelude.Fields

fn ringEnergy(s, coupling, field) ↦ totalEnergy(s, coupling, field)
fn ringMagnet(s) ↦ magnetisation(s)
fn glauberStep(s, coupling, field, temperature, u, v) ↦ flipped(s, coupling, field, temperature, u, v)

using Examples;

var up    = new[] { 1d, 1d, 1d, 1d };
var mixed = new[] { 1d, -1d, 1d, 1d };

var e1 = Lattice.RingEnergy(up, 1.0, 0.0);
var e2 = Lattice.RingEnergy(mixed, 1.0, 0.0);
var m  = Lattice.RingMagnet(mixed);
var s1 = Lattice.GlauberStep(up, 1.0, 0.0, 1.0, 0.3, 0.0);
var s2 = Lattice.GlauberStep(up, 1.0, 0.0, 1.0, 0.3, 0.5);

Output (verified):

energy = -4
energy2 = 0
magnet = 0.5
flip hot = [1, -1, 1, 1]
flip no = [1, 1, 1, 1]

How big an epidemic gets, and how likely it is to fizzle

With $\beta = 0.3$ and $\gamma = 0.1$, $R_0 = 3$: about 94% of the population is eventually infected in a major outbreak, and a single case dies out by chance a third of the time. At $\beta = 0.05$ ($R_0 = 0.5$) there is no major outbreak.

open Prelude.Fields

fn basicReproduction(beta, gamma) ↦ reproduction(beta, gamma)
fn attackRate(beta, gamma) ↦ finalSize(beta, gamma)
fn extinctionChance(beta, gamma, start) ↦ diesOut(beta, gamma, start)

using Examples;

var r0      = Lattice.BasicReproduction(0.3, 0.1);
var final   = Lattice.AttackRate(0.3, 0.1);
var minor   = Lattice.AttackRate(0.05, 0.1);
var fizzles = Lattice.ExtinctionChance(0.3, 0.1, 1.0);

Output (verified, rounded to six digits):

R0 = 3
final = 0.94048
final sub = 0
diesOut = 0.333333

Stepping a stochastic SIR epidemic

From $[99, 1, 0]$ in a population of 100, the infection rate is 0.297 and the recovery rate 0.1, so infection has share about 0.748. A draw of 0.5 gives an infection, a draw of 0.9 a recovery; once nobody is infected the state no longer changes.

open Prelude.Fields

fn sirStep(x, beta, gamma, size, u) ↦ epidemicStep(x, beta, gamma, size, u)

using Examples;

var infection = Lattice.SirStep(new[] { 99d, 1d, 0d }, 0.3, 0.1, 100.0, 0.5);
var recovery  = Lattice.SirStep(new[] { 99d, 1d, 0d }, 0.3, 0.1, 100.0, 0.9);
var over      = Lattice.SirStep(new[] { 99d, 0d, 1d }, 0.3, 0.1, 100.0, 0.9);

Output (verified):

sir = [98, 2, 0]
sir rec = [99, 0, 1]
sir over = [99, 0, 1]

Will a mutant take over?

One mutant copy in a Moran population of 10: neutral, it fixes with probability 1/10; with a 10% fitness advantage, about 0.148.

open Prelude.Fields

fn moranOf(count, size, selection) ↦ moranFixation(count, size, selection)

using Examples;

var neutral    = Lattice.MoranOf(1.0, 10.0, 0.0);
var advantaged = Lattice.MoranOf(1.0, 10.0, 0.1);

Output (verified, rounded to six digits):

moran neutral = 0.1
moran s = 0.14795

One generation of genetic drift

An allele at frequency 0.5 in a population of 100 with a normal draw of $z = 1$ moves by $\sqrt{0.25/100} = 0.05$. The exact resample with four uniform draws counts how many fall below 0.5.

open Prelude.Fields

fn wrightFisher(p, z, selection, mutation, population) ↦ nextFrequency(p, z, selection, mutation, population)
fn resampled(p, u, selection, mutation) ↦ binomialCount(p, u, selection, mutation)

using Examples;

var diffusion = Lattice.WrightFisher(0.5, 1.0, 0.0, 0.0, 100.0);
var exact     = Lattice.Resampled(0.5, new[] { 0.1, 0.4, 0.6, 0.9 }, 0.0, 0.0);

Output (verified):

wf = 0.55
resampled = 0.5

Bounding an insurer’s chance of ruin

Claims of mean size 1 arrive at rate 1 and the premium rate is 1.2, a 20% safety loading. With a starting surplus of 10, the adjustment coefficient is $1 - 1/1.2 = 1/6$ and the probability of ever being ruined is at most $e^{-10/6}$.

open Prelude.Fields

fn ruinBound(surplus, premium, rate, claimMean) ↦ lundberg(surplus, premium, rate, claimMean)
fn loadingOf(premium, rate, claimMean) ↦ safetyLoading(premium, rate, claimMean)

using Examples;

var bound   = Lattice.RuinBound(10.0, 1.2, 1.0, 1.0);
var loading = Lattice.LoadingOf(1.2, 1.0, 1.0);

Output (verified, rounded to six digits):

lundberg = 0.188876
loading = 0.2