# Deriving the Distributions from a Coin and a Clock

An online shop keeps asking questions that sound unrelated. Of the next 200 visitors, how many will buy? How long until the next order? How many support tickets will arrive between noon and one? Is a 12% conversion in a small test real, or luck? What will today's revenue be, and why is the average customer worth more than the typical one? Each question has a named distribution, and the names are usually learned as a list. They are better learned as one process looked at in different ways.

<figure class="mfig ds-figure ds-wide" id="ds-process" aria-label="One stream of random events, with the number of events per window and the gaps between events">
<div class="ds-plot"></div>
<figcaption><b>One process, two distributions.</b> The dots are events. Count the dots in each window and the counts follow a Poisson law; measure the gaps and they follow an exponential one. Drag the handle to change the rate.</figcaption>
</figure>

Three ingredients are enough: a coin flip, a clock tick and a draw that is equally likely anywhere in $[0,1]$. Three moves build everything else from them. **Composing** adds things up, multiplies them, waits for them, or lets the ticks get small. **Marginalising** averages over something that was random but unseen. **Relaxing an assumption** changes one rule of the ticks and sees what follows; it is used once, for the Weibull. The named distributions are what those moves produce, and the last figure is a map of which comes from which. One distribution does most of the connecting: the gamma, which is a wait for several events, returns in the beta, the negative binomial, Student's $t$ and the Tweedie. Each section opens with a question the shop might ask, and answers it once the distribution is built. One colour code runs through the text and the figures: <span class="k-obs">blue</span> is what we observe, <span class="k-hid">copper</span> is a hidden quantity such as a rate or a probability, and <span class="k-der">violet</span> is the distribution we derive.

## Assumptions

1. **Ticks.** Time is cut into small steps of length $\Delta$. In each step an event happens or it does not.
2. **Independence.** What happens in one step does not change the chance in another.
3. **A constant chance.** In every step the chance of an event is $\lambda\Delta$: it is proportional to the length of the step, with the same $\lambda$ everywhere. We call $\lambda$ the rate.
4. **Hidden quantities have their own laws.** When a rate or a probability is unknown, we treat it as a draw from a distribution and average over it.
5. **Sizes are drawn independently of counts.** When each event carries a random size, the sizes do not depend on how many events there were.

Assumptions 2 and 3 are where real data most often disagree, and the last section comes back to them; assumptions 4 and 5 are used where the hidden rates and the sizes appear.

## Adding coin flips gives the binomial

**Teaser.** Of 200 visitors who see a new checkout, 31 buy. If the true conversion rate is 12%, how surprising is 31?

A **Bernoulli** variable is one yes-or-no outcome with chance $p$ of yes. It has mean $p$ and variance $p(1-p)$. It is the smallest thing that can be random, and everything below is built from copies of it.

Add $n$ independent copies and look at the count of yeses $k$. There are $\binom nk$ ways to choose which flips were yes, and each has probability $p^k(1-p)^{n-k}$, so

$$
\htmlClass{eq-der}{\Pr(K=k)}=\binom nk p^k(1-p)^{n-k}.
$$

This is the **binomial** law. Nothing was fitted: it is the sum of $n$ Bernoulli variables, counted by hand.

<figure class="mfig ds-figure" id="ds-binom" aria-label="The binomial distribution built from the number of coin flips">
<div class="ds-plot"></div>
<figcaption><b>A sum of coin flips.</b> Each bar is the chance of that many yeses in $n$ flips, each with chance 0.35. Drag the handle to add flips: the bars move right and spread, and the mean $np$ and the variance $np(1-p)$ grow in step.</figcaption>
</figure>


<details class="answer"><summary>Worked example: is 31 buyers out of 200 surprising?</summary>

The mean is $200\times0.12=24$ buyers and the standard deviation is $\sqrt{200\times0.12\times0.88}=4.6$, so 31 is $(31-24)/4.6=1.5$ standard deviations above the mean. The exact binomial gives $\Pr(K\ge31)=0.082$, about one chance in twelve. That is unusual but not decisive: a 12% rate is not ruled out by 31 buyers, and 31 of 200 is an observed rate of 15.5%.

</details>

## A clock is a coin flipped faster: the exponential

**Teaser.** A support desk gets a ticket every three minutes on average. The last one came nine minutes ago. Is the next one overdue?

Keep flipping, and ask how many flips it takes to get the first yes. It takes $k$ flips if the first $k-1$ are no and the next is yes, so $\Pr(K=k)=(1-p)^{k-1}p$. This is the **geometric** law.

Now make the flips fast. Cut one unit of time into $m$ ticks and let the chance per tick be $p=\lambda/m$, which keeps the rate at $\lambda$ per unit. The wait is longer than time $t$ if the first $mt$ ticks all fail:

$$
\Pr(T>t)=\Big(1-\frac{\lambda}{m}\Big)^{mt}\ \xrightarrow{\ m\to\infty\ }\ e^{-\lambda t}.
$$

The limit is a continuous wait with density $\lambda e^{-\lambda t}$, the **exponential**. It carries one property that all the later constructions use: it has no memory. Given that you have already waited $s$, the chance of waiting $t$ more is the same as it was at the start, because each tick is independent of the ones before it.

<figure class="mfig ds-figure" id="ds-geometric" aria-label="The geometric waiting time converging to the exponential as ticks get shorter">
<div class="ds-plot"></div>
<figcaption><b>Shrink the tick and the geometric becomes the exponential.</b> Bars are the geometric waiting time, rescaled to time units. The curve is the exponential. Drag the handle to cut each unit of time into more ticks.</figcaption>
</figure>

<details class="more"><summary>Derivation: the limit and the lack of memory</summary>

Take logarithms: $mt\log(1-\lambda/m)$. For large $m$, $\log(1-\lambda/m)=-\lambda/m+O(m^{-2})$, so the product tends to $-\lambda t$. For the lack of memory, $\Pr(T>s+t\mid T>s)=e^{-\lambda(s+t)}/e^{-\lambda s}=e^{-\lambda t}$, which does not depend on $s$.

</details>


<details class="answer"><summary>Worked example: is the next ticket overdue?</summary>

A mean gap of three minutes is a rate $\lambda=1/3$ per minute. A nine-minute wait was unlikely from the start: $\Pr(T>9)=e^{-9/3}=0.050$. But the desk has already waited, and the exponential forgets: the chance of waiting at least three more minutes is $\Pr(T>9+3\mid T>9)=\Pr(T>3)=e^{-1}=0.368$, exactly what it was nine minutes ago, and the expected extra wait is still three minutes. A quiet spell is not a reason to expect a ticket sooner.

</details>

## Counting events in a window gives the Poisson

**Teaser.** The desk gets 20 tickets an hour on average, and each agent can handle 5 an hour. Four agents cover the average. How often will an hour bring more than they can handle, and what do a fifth and a sixth agent buy?

Use the same ticks and count the events in a window of length $t$. That is $n=t/\Delta$ coin flips, each yes with chance $\lambda\Delta$, so the count is binomial with $n=t/\Delta$ and mean $\mu=\lambda t$. Let the ticks shrink with the mean held fixed. The binomial becomes

$$
\htmlClass{eq-der}{\Pr(N=k)}=\frac{\mu^k e^{-\mu}}{k!},
$$

the **Poisson** law. Its mean and its variance are both $\mu$, which is the fingerprint to check against when a count is said to be Poisson.

<figure class="mfig ds-figure" id="ds-poisson" aria-label="The binomial with many small chances converging to the Poisson with the same mean">
<div class="ds-plot"></div>
<figcaption><b>Many rare chances give a Poisson.</b> Bars are the binomial with $n$ trials and chance $3/n$ each; the dots are the Poisson with the same mean, 3. Drag the handle to increase $n$.</figcaption>
</figure>

With $n=10$ trials the total variation distance between the two is $0.086$; with $n=100$ it is $0.0076$, and with $n=1000$ it is $0.0008$.

<details class="more"><summary>Derivation: the binomial becomes the Poisson</summary>

With $p=\mu/n$, $\binom nk p^k(1-p)^{n-k}=\frac{n(n-1)\cdots(n-k+1)}{n^k}\cdot\frac{\mu^k}{k!}\cdot(1-\mu/n)^{n}\cdot(1-\mu/n)^{-k}$. As $n\to\infty$ the first factor tends to 1, the third to $e^{-\mu}$, and the fourth to 1, leaving $\mu^ke^{-\mu}/k!$.

</details>

The exponential and the Poisson are not two models. They are two ways of reading the same ticks, one measuring the gaps and one counting the events: the chance of no event in a window of length $t$ is $e^{-\mu}$ with $\mu=\lambda t$, and it is also the chance that the first wait exceeds $t$.


<details class="answer"><summary>Worked example: how many agents?</summary>

The hour is a Poisson count with mean 20. Four agents have capacity 20, and 21 or more tickets arrive in $44\%$ of hours, so staffing for the mean fails nearly half the time. A fifth agent (capacity 25) is overrun by 26 or more tickets in $11.2\%$ of hours, and a sixth (capacity 30) by 31 or more in $1.3\%$. The Poisson says how much headroom buys how much safety.

</details>

## Waiting for $k$ events gives the gamma

**Teaser.** A flash sale is about to open. Orders arrive at 4 per minute on average. How long until the fifth order, and how much does that vary?

Three facts about the gamma are enough for everything that follows, and each has a plain story.

**Fact 1: waiting for $k$ events is $k$ exponential waits added up.** After each event the clock starts afresh, because the exponential has no memory, so the wait for the $k$-th event is a sum of $k$ independent exponential waits. For two, the density of the total at $t$ adds up every way of splitting it: the first wait ends at $s$, the second ends $t-s$ later.

$$
f_2(t)=\int_0^t\lambda e^{-\lambda s}\,\lambda e^{-\lambda(t-s)}\,ds=\lambda^2e^{-\lambda t}\int_0^t ds=\lambda^2\,t\,e^{-\lambda t}.
$$

The two exponentials combine into $e^{-\lambda t}$ whatever $s$ is, so only the length of the interval is left. Each further wait integrates one more power of $t$ and divides by it, and for $k$ waits the density is

$$
\htmlClass{eq-der}{f_k(t)}=\frac{\lambda^k t^{k-1}e^{-\lambda t}}{(k-1)!}.
$$

<figure class="mfig ds-figure" id="ds-gamma" aria-label="Simulated waiting times for the k-th event against the gamma density">
<div class="ds-plot"></div>
<figcaption><b>The wait for the $k$-th event.</b> The bars are 400 simulated waits, each the sum of $k$ exponential gaps; the curve is the gamma density with shape $k$. Drag the handle to wait for more events.</figcaption>
</figure>

<details class="more"><summary>Derivation: adding one more wait</summary>

If $f_k$ is the density of the first $k$ waits, the next wait is independent and exponential, so $f_{k+1}(t)=\int_0^t f_k(s)\,\lambda e^{-\lambda(t-s)}ds=\lambda^{k+1}e^{-\lambda t}\int_0^t\frac{s^{k-1}}{(k-1)!}ds=\frac{\lambda^{k+1}t^{k}e^{-\lambda t}}{k!}$. The case $k=1$ is the exponential, and induction gives every $k$.

</details>

<details class="more"><summary>Derivation: the same density from the Poisson count</summary>

Waiting more than $t$ for the $k$-th event is the same as seeing fewer than $k$ events in the window, so $\Pr(T_k>t)=S_k(t)=\sum_{j<k}(\lambda t)^je^{-\lambda t}/j!$. Then $\frac{d}{dt}\big[(\lambda t)^je^{-\lambda t}/j!\big]=\lambda\big[(\lambda t)^{j-1}/(j-1)!-(\lambda t)^j/j!\big]e^{-\lambda t}$, and in the sum over $j<k$ every term cancels except $-\lambda(\lambda t)^{k-1}e^{-\lambda t}/(k-1)!$. The density is minus the derivative of $S_k$, which is the formula above. The two routes agree, as they must.

</details>

**Fact 2: the shape counts events and the rate sets the speed of the clock.** A faster clock only rescales time: the wait for $k$ events at rate $\lambda$ is the wait at rate 1 divided by $\lambda$. The mean is $k/\lambda$ and the variance is $k/\lambda^2$, since each wait has mean $1/\lambda$ and variance $1/\lambda^2$ and the waits add. The exponential is the case $k=1$. Waiting for $k_1$ events and then $k_2$ more is waiting for $k_1+k_2$, so independent gammas with the same rate add their shapes.

**Fact 3: the shape can be any positive number, and the integrals are free.** Replace $(k-1)!$ by $\Gamma(k)$ and the density is a valid distribution for every positive shape $a$ and rate $b$:

$$
\text{Gamma}(t;a,b)=\frac{b^a}{\Gamma(a)}\,t^{a-1}e^{-bt},\qquad\int_0^\infty t^{a-1}e^{-bt}\,dt=\frac{\Gamma(a)}{b^a}.
$$

The second formula is the point. A power of $t$ times an exponential in $t$ integrates in closed form, and multiplying two such kernels gives another one. Whenever a gamma is averaged against something of the same kind, the integral costs nothing.

That is why the gamma recurs. A Poisson count with a gamma rate gives the negative binomial. The share of two gammas is a beta. A normal with a gamma precision is Student's $t$, and sums of squared normals are gamma. A Poisson number of gamma-sized claims is the Tweedie. Each of those sections uses only the three facts above.

<details class="answer"><summary>Worked example: how long until the fifth order?</summary>

At 4 orders a minute the fifth order takes $5/4=1.25$ minutes on average, with a standard deviation of $\sqrt5/4=0.56$ minutes. The chance it arrives within the first minute is the chance of at least 5 orders in one minute: $1-\Pr(\text{Poisson}(4)\le4)=0.371$. The chance it takes more than two minutes is $\Pr(\text{Poisson}(8)\le4)=0.0996$, about one in ten.

</details>


## A probability you do not know: the beta

**Teaser.** An A/B test shows 12 conversions in 80 visits. Before running more traffic, what should we believe about the true rate?

The exponential and the gamma describe time. A different question is about a probability itself. Drop $n$ points uniformly on $[0,1]$ and look at the $k$-th from the left. For it to sit near $x$, exactly $k-1$ points must fall below $x$, one at $x$, and $n-k$ above:

$$
\htmlClass{eq-der}{f(x)}=\frac{n!}{(k-1)!\,(n-k)!}\,x^{k-1}(1-x)^{n-k}.
$$

This is the **beta** law with parameters $k$ and $n-k+1$. The same family appears from the other side. If a coin has unknown chance $x$ of yes, our prior is flat, and we see $k$ yeses and $n-k$ noes, then the posterior for $x$ is proportional to $x^k(1-x)^{n-k}$, a beta with parameters $k+1$ and $n-k+1$: the same family as the order statistic, shifted by one in the first parameter by the flat prior. A beta is a distribution over a probability.

<figure class="mfig ds-figure" id="ds-beta" aria-label="The k-th smallest of nine uniform points and its beta density">
<div class="ds-plot"></div>
<figcaption><b>The $k$-th smallest of nine points.</b> The dots are nine uniform draws, the $k$-th highlighted. The bars are where that same order statistic lands over 2,000 repeats; the curve is the beta density. Drag the handle to choose $k$.</figcaption>
</figure>


<details class="answer"><summary>Worked example: what is the rate after 12 conversions in 80 visits?</summary>

A flat prior (assumption 4: the unknown rate gets a law of its own) and 12 yeses with 68 noes give a beta posterior with parameters 13 and 69. Its mean is $13/82=0.159$ and its central 90% interval runs from $0.098$ to $0.229$. If the business needs the rate to be above 10%, the posterior probability is $0.94$; above 20%, $0.15$. The interval is wide and honest about it: 80 visits cannot say more.

</details>

**Shares of gammas: a beta again, and the Dirichlet.** *Teaser: customers choose among three plans, and the shares vary from week to week.* Take independent gamma variables and divide each by their sum. For two of them, $X\sim\text{Gamma}(a,1)$ and $Y\sim\text{Gamma}(b,1)$, the share $X/(X+Y)$ is a $\text{Beta}(a,b)$ and is independent of the total $X+Y\sim\text{Gamma}(a+b,1)$. With several gammas of shapes $\alpha_1,\dots,\alpha_m$ the same change of variables gives the **Dirichlet** law, a distribution over shares that add to one; each share on its own is a $\text{Beta}(\alpha_i,\alpha_0-\alpha_i)$ with $\alpha_0=\sum_j\alpha_j$. It is to the multinomial what the beta is to the binomial.

<details class="more"><summary>Derivation: the share of two gammas</summary>

The joint density is $f(x,y)=x^{a-1}y^{b-1}e^{-(x+y)}/(\Gamma(a)\Gamma(b))$. Let $s=x+y$ and $u=x/s$, so $x=us$, $y=(1-u)s$ and the Jacobian is $s$. Then $f(u,s)=\frac{u^{a-1}(1-u)^{b-1}}{\Gamma(a)\Gamma(b)}\,s^{a+b-1}e^{-s}$. It factorises into a function of $u$ alone, which is a $\text{Beta}(a,b)$ density after normalising, and a function of $s$ alone, a $\text{Gamma}(a+b,1)$ density.

</details>

<details class="answer"><summary>Worked example: how much does one plan's share swing?</summary>

Take three plans with Dirichlet parameters $(6,3,1)$: average shares $0.6$, $0.3$ and $0.1$, with total concentration 10. The third plan's share is a $\text{Beta}(1,9)$. The chance that it exceeds a 20% share in a given week is $0.8^9=0.134$, even though its average is 10%. The first plan's share has a standard deviation of $\sqrt{0.6\times0.4/11}=0.148$.

</details>

## Marginalisation: hide a parameter and a new distribution appears

Composing built things up from the same ticks. Marginalising takes a distribution that depends on a parameter $\theta$ and, when $\theta$ is itself random and unseen, averages it out:

$$
p(x)=\int p(x\mid\theta)\,p(\theta)\,d\theta .
$$

The result is a new distribution whose shape neither factor had. This is assumption 4 at work, and each one below follows from a story in which something is hidden.

**Thinning.** *Teaser: visitors arrive at 40 an hour and 3% buy; how many buyers per hour?* A Poisson number of visitors arrives, with mean $\lambda$, and each buys independently with chance $p$. The number of buyers given $n$ visitors is binomial; averaging over the unseen $n$ gives a Poisson with mean $\lambda p$. Marginalising a binomial over a Poisson count gives a Poisson again.

<details class="more"><summary>Derivation: thinning a Poisson</summary>

With $m=n-k$ visitors who do not buy, $\sum_{n\ge k}\frac{\lambda^ne^{-\lambda}}{n!}\binom nk p^k(1-p)^{n-k}=e^{-\lambda}\frac{(\lambda p)^k}{k!}\sum_{m\ge0}\frac{(\lambda(1-p))^m}{m!}=e^{-\lambda}\frac{(\lambda p)^k}{k!}e^{\lambda(1-p)}=\frac{(\lambda p)^ke^{-\lambda p}}{k!}.$

</details>

<details class="answer"><summary>Worked example: buyers per hour</summary>

Visitors at 40 an hour and a 3% buy rate give a Poisson with mean $40\times0.03=1.2$ buyers an hour. The chance of an hour with no buyer is $e^{-1.2}=0.301$, and of three or more buyers $0.121$.

</details>

**An unknown rate: the negative binomial.** *Teaser: daily orders swing more than a Poisson allows, so a forecast with a Poisson interval keeps being wrong. Why?* Suppose each day's rate $\Lambda$ is itself a draw from a gamma law with shape $r$ and mean $m$, and the day's count is Poisson given the rate. Then

$$
\Pr(N=k)=\int_0^\infty \frac{\lambda^ke^{-\lambda}}{k!}\,\htmlClass{eq-hid}{f_\Lambda(\lambda)}\,d\lambda
=\frac{\Gamma(k+r)}{k!\,\Gamma(r)}\Big(\frac{r}{r+m}\Big)^r\Big(\frac{m}{r+m}\Big)^k .
$$

This is the **negative binomial**. The mean is still $m$, but the variance is $m+m^2/r$: more than the Poisson's $m$. The extra spread is exactly the variance of the hidden rate. A count that looks overdispersed relative to a Poisson is what a Poisson with a wandering rate looks like.

<figure class="mfig ds-figure ds-wide" id="ds-negbin" aria-label="A gamma-distributed rate and the negative binomial count it produces, against a Poisson with the same mean">
<div class="ds-plot"></div>
<figcaption><b>A wandering rate fattens the count.</b> Left: the distribution of the rate, with mean 4. Right: the counts it produces (bars) against a Poisson with the same mean (dots). Drag the handle to make the rate more or less stable: a large $r$ pins it down.</figcaption>
</figure>

For $r=2$ and mean 4 the variance is $4+16/2=12$, three times the Poisson's. For $r=5$ it is $7.2$. The shop's forecast was too narrow because it assumed a rate that was the same every day.

<details class="answer"><summary>Worked example: a forecast interval for daily orders</summary>

Suppose the shop averages 40 orders a day and the daily rate varies with $r=8$, a relative spread of $1/\sqrt8=35\%$ in the rate. The variance is $40+40^2/8=240$, a standard deviation of $15.5$ orders, against $6.3$ for a Poisson. The central 95% range of daily orders is 28 to 53 under a Poisson and 15 to 75 under the negative binomial. A day above 60 orders has probability $0.0012$ under the Poisson, so it would be called a fluke, and $0.101$ under the negative binomial, about one day in ten.

</details>

**An unknown probability: the beta-binomial.** *Teaser: campaigns each reach 100 visitors, with average conversion 10%, and the number of buyers varies more from campaign to campaign than a binomial allows.* If each campaign has its own conversion chance, drawn from a beta law with parameters $a,b$, the count of buyers is a binomial averaged over a beta: the **beta-binomial**. With $a=2$ and $b=18$ (mean 0.1) and 100 visitors, the mean count is 10 and the variance is $51.4$, against $9$ for a single fixed chance: a factor of $1+(n-1)/(a+b+1)=5.7$. Campaigns differ, so counts across them are more spread out than coins would be.

<details class="answer"><summary>Worked example: how often does a campaign overperform?</summary>

With 100 visitors and an average conversion of 10%, a campaign with 20 or more buyers has probability $0.002$ if every campaign shared one chance, and $0.106$ under the beta-binomial with $a=2$, $b=18$: about one campaign in ten. What looks like a rare lucky campaign under a single chance is ordinary once campaigns differ.

</details>

## Adding up many parts gives the normal

**Teaser.** What will today's revenue be, and how far off could the forecast be?

Daily revenue is the sum of many separate orders. The binomial was already one such sum, and its bars become a bell as $n$ grows. The reason is general: a sum of many independent contributions, none dominant, has a distribution that forgets the shape of its parts and keeps only their mean and variance. After centring and scaling it approaches the **normal** law, with density $e^{-x^2/2}/\sqrt{2\pi}$.

<details class="more"><summary>Derivation: why only the mean and variance survive</summary>

Let the parts have mean $m$, variance $v$ and $j$-th cumulant $\kappa_j$ (the first is $m$, the second $v$). Cumulants add over independent parts, and the $j$-th cumulant scales with the $j$-th power of the scale. For the centred and scaled sum $Z_n=(S_n-nm)/\sqrt{nv}$ the first cumulant is 0, the second is 1, and for $j\ge3$ the $j$-th is $n\kappa_j/(nv)^{j/2}=\kappa_j/(v^{j/2}n^{j/2-1})$, which tends to 0. A distribution whose cumulants beyond the second all vanish is the normal; this is the standard route to the central limit theorem, and it needs those cumulants to exist. For $j=3$ it also gives the skewness $\kappa_3/(v^{3/2}\sqrt n)$, which is $2/\sqrt n$ for exponentials.

</details>

The shape of the parts is forgotten slowly. Orders are skewed, like exponentials. The sum of $n$ unit exponentials is a gamma with shape $n$, and its skewness is $2/\sqrt n$: $2$ for one order, $0.63$ for ten, $0.2$ for a hundred. A normal is symmetric, so the approximation is crude at ten orders and decent at a hundred.

<figure class="mfig ds-figure" id="ds-clt" aria-label="Standardised sums of exponential waits approaching a normal density as more are added">
<div class="ds-plot"></div>
<figcaption><b>Sums forget the shape of their parts.</b> The bars are standardised sums of $n$ skewed exponential waits; the curve is the normal. Drag the handle to add more terms.</figcaption>
</figure>


<details class="answer"><summary>Worked example: revenue forecast for a day</summary>

Take 400 orders a day on average, each worth 50 dollars on average with a standard deviation of 30 dollars. The expected revenue is $400\times50=20{,}000$ dollars. The number of orders is itself random, so the variance is $400\,(30^2+50^2)=1.36\times10^6$: the variance of a Poisson number $N$ of orders with sizes $X$ is $\mathbb E[N]\operatorname{Var}X+\operatorname{Var}N\,(\mathbb EX)^2$, and both $\mathbb E[N]$ and $\operatorname{Var}N$ equal 400. That is a standard deviation of about 1,166 dollars, twice the 600 dollars a fixed count of 400 orders would give. The normal 95% range is 17,700 to 22,300 dollars. With hundreds of orders the normal is reasonable; on a quiet day with ten, the skew from the previous paragraph would show.

</details>

## Multiplying many factors gives the lognormal

**Teaser.** Order values and session lengths are heavily skewed, and the average customer is worth more than the typical customer. Why?

Many quantities grow by multiplication. A customer's spend is the result of a series of factors: a promotion, a repeat visit, a bigger basket. A product of positive factors becomes a sum on the log scale, $\log(F_1\cdots F_n)=\sum\log F_i$, and a sum of many small terms is approximately normal. So the product is approximately the exponential of a normal: the **lognormal**, with density $\frac1x$ times a normal density in $\log x$. It is skewed to the right and has a long tail.

The mean and the typical value part ways. Take factors drawn evenly between 0.5 and 1.5: each has mean 1, so the product also has mean 1. But the log of one factor has mean $-0.045$ and variance $0.095$. The median of the product after 16 factors is $e^{16\times(-0.045)}=0.485$: a typical customer ends up at half the starting value, while the average stays at 1 because a few customers end up far above.

<figure class="mfig ds-figure ds-wide" id="ds-lognormal" aria-label="Products of random factors against the lognormal density, on the original and the logarithmic scale">
<div class="ds-plot"></div>
<figcaption><b>Products are skewed, their logs are not.</b> Left: the products of $n$ random factors, on the original scale, with the lognormal density. Right: the same values on the log scale, with a normal density. Drag the handle to multiply more factors.</figcaption>
</figure>


<details class="answer"><summary>Worked example: typical against average order</summary>

Let order values be lognormal with a typical (median) order of 40 dollars and a log-scale spread $\sigma=0.9$. The mean is $40\,e^{0.9^2/2}=60$ dollars: 1.5 times the median, pulled up by the tail. Orders above 200 dollars have probability $0.037$, and the top 1% of orders start at 325 dollars. Reporting only the mean of 60 dollars would describe almost no actual order: report the median and a high percentile.

</details>

## Not knowing the noise widens the tails: Student's t

**Teaser.** An experiment measures uplift on 8 customers per arm. How wide should the confidence interval be, when the noise level itself has to be estimated from the same eight?

A normal measurement has a spread $\sigma$. When $\sigma$ is not known it has to be estimated, and the estimate is itself random. The gamma gives its law in two steps, and a third step divides one by the other.

**Step 1: a squared normal is a gamma.** Let $Z$ be standard normal. The chance that $Z^2\le y$ is the chance that $-\sqrt y\le Z\le\sqrt y$, so the density of $Z^2$ at $y$ is $2\varphi(\sqrt y)\cdot\frac1{2\sqrt y}=\dfrac{y^{-1/2}e^{-y/2}}{\sqrt{2\pi}}$. Since $\Gamma(\tfrac12)=\sqrt\pi$, this is exactly the gamma density with shape $\tfrac12$ and rate $\tfrac12$. Shapes add for equal rates, so a sum of $\nu$ independent squared standard normals is $\text{Gamma}(\tfrac\nu2,\tfrac12)$, known as the chi-squared law with $\nu$ degrees of freedom. This is the story the half-integer shapes were missing.

<figure class="mfig ds-figure" id="ds-chi" aria-label="Sums of squared standard normals against the gamma density with shape nu over two and rate one half">
<div class="ds-plot"></div>
<figcaption><b>Squares of normals add up to a gamma.</b> The bars are 2,000 sums of $\nu$ squared standard normals; the curve is the gamma density with shape $\nu/2$ and rate $\tfrac12$. Drag the handle to add more squares.</figcaption>
</figure>

**Step 2: an estimated variance is a gamma with mean one.** Suppose the variance is estimated from $\nu$ independent standardised errors, so that $s^2=\sigma^2\cdot\frac1\nu\sum Z_i^2$. Dividing a gamma by $\nu$ multiplies its rate by $\nu$ (Fact 2), so $V=s^2/\sigma^2$ is $\text{Gamma}(\tfrac\nu2,\tfrac\nu2)$. Its mean is 1 and its variance is $2/\nu$: with eight observations the estimated variance is routinely far from the true one. For one group of $n$ observations centred on their own mean, $\nu=n-1$.

**Step 3: divide.** The statistic $(\bar x-\mu)/(s/\sqrt n)$ is a standard normal $Z$ divided by $\sqrt V$, and for normal data the sample mean and the sample variance are independent. Given $V=\tau$ the statistic is normal with precision $\tau$, so averaging over the gamma gives

$$
\htmlClass{eq-der}{p(x)}=\int_0^\infty \htmlClass{eq-obs}{\mathcal N\!\big(x;0,1/\tau\big)}\ \htmlClass{eq-hid}{\text{Gamma}\big(\tau;\tfrac\nu2,\tfrac\nu2\big)}\,d\tau=\frac{\Gamma(\frac{\nu+1}2)}{\sqrt{\nu\pi}\,\Gamma(\frac\nu2)}\Big(1+\frac{x^2}{\nu}\Big)^{-\frac{\nu+1}{2}},
$$

the **Student $t$** law with $\nu$ degrees of freedom. Its tails are heavier than the normal's, and they thin out toward the normal as $\nu$ grows.

<details class="more"><summary>Derivation: why the integral closes</summary>

Multiply the two factors. The normal contributes $\sqrt{\tau/2\pi}\,e^{-\tau x^2/2}$ and the gamma contributes $\frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)}\tau^{\nu/2-1}e^{-\nu\tau/2}$. In $\tau$ the product is $\tau^{(\nu+1)/2-1}e^{-\tau(\nu+x^2)/2}$ times constants, which is a gamma kernel with shape $\frac{\nu+1}2$ and rate $\frac{\nu+x^2}2$. By Fact 3 its integral is $\Gamma(\frac{\nu+1}2)\big(\frac{2}{\nu+x^2}\big)^{(\nu+1)/2}$. Collecting the constants gives $\frac{\Gamma((\nu+1)/2)}{\sqrt\pi\,\Gamma(\nu/2)}\nu^{\nu/2}(\nu+x^2)^{-(\nu+1)/2}$, and $\nu^{\nu/2}(\nu+x^2)^{-(\nu+1)/2}=\nu^{-1/2}(1+x^2/\nu)^{-(\nu+1)/2}$.

</details>

Read the integral as a story. Each measurement comes from a normal whose precision was drawn from the gamma. When the precision happens to be small, the normal is wide, and that rare wide normal supplies the far tail. The same construction, used on purpose as a model of noise whose level varies, is a common way to make a regression robust to outliers.

<figure class="mfig ds-figure" id="ds-tmix" aria-label="Eleven normal densities with precisions drawn from a gamma, and the Student t density they average to">
<div class="ds-plot"></div>
<figcaption><b>One $t$ is many normals.</b> Each thin curve is a normal at one typical precision, from the 5th to the 95th percentile of the gamma; the heavy curve is the $t$ they average to. Drag the handle to change the degrees of freedom $\nu$: with few of them the precisions spread widely.</figcaption>
</figure>

<figure class="mfig ds-figure" id="ds-t" aria-label="The Student t density against the normal, for a chosen number of degrees of freedom">
<div class="ds-plot"></div>
<figcaption><b>Not knowing the noise costs width.</b> The solid curve is the $t$; the dots are the normal. Drag the handle to change the degrees of freedom $\nu$. The tail beyond $\pm3$ is shaded.</figcaption>
</figure>

With $\nu=3$ the chance of landing beyond $\pm3$ is $0.058$; for the normal it is $0.0027$, twenty-one times smaller. An interval for an uplift measured on eight customers has to be wider than the normal one says, and the $t$ law says by how much.

<details class="answer"><summary>Worked example: a 95% interval from eight customers per arm</summary>

Say the pooled standard deviation is 20 and the measured uplift is $+12$. The standard error is $20\sqrt{2/8}=10$. The normal 95% interval is $12\pm1.96\times10$, or $(-7.6,\ 31.6)$. The variance is pooled from two groups of eight, so it has $2\times8-2=14$ degrees of freedom, and the 97.5% point of the $t$ with 14 degrees is $2.145$: the interval is $12\pm21.4$, or $(-9.4,\ 33.4)$, 9% wider. Both include zero here, so the data do not settle that the uplift is positive, and the $t$ interval says how much less sure to be.

</details>

## A lump at zero and a skewed rest: the Tweedie

**Teaser.** Revenue per customer per month: most customers spend nothing, and the rest spread widely. What distribution has an exact zero and a skewed positive part?

Now put the two operations together. A random number $N$ of orders arrives, Poisson with mean $\lambda$, and each brings a random size $X_i$, independent of how many orders there were (assumption 5), gamma with shape $\alpha$ and scale $\theta$ (so a size has mean $\alpha\theta$). The total is $Y=X_1+\dots+X_N$. Two things follow.

If no order arrives, the total is exactly zero, which has probability $e^{-\lambda}$. If $n\ge1$ arrive, the sum of $n$ gamma sizes is gamma with shape $n\alpha$, so for $y>0$

$$
f(y)=\sum_{n=1}^\infty \frac{\lambda^ne^{-\lambda}}{n!}\ \htmlClass{eq-der}{\text{Gamma}(y;\,n\alpha,\ \theta)} .
$$

This is a **compound Poisson–gamma** law: a lump of probability at zero and a smooth density on the positive numbers. Its mean is $\lambda\alpha\theta$, and its variance, by adding the variance of the count to the variance of the sizes, is $\lambda\alpha(\alpha+1)\theta^2$.

The Tweedie family is this law written in terms of a mean $\mu$, a dispersion $\phi$ and a power $p$ so that the variance is a power of the mean, $\operatorname{Var}=\phi\mu^p$ (here $p$ is the power, not a Bernoulli chance). The map is

$$
\lambda=\frac{\mu^{2-p}}{\phi\,(2-p)},\qquad
\alpha=\frac{2-p}{p-1},\qquad
\theta=\phi\,(p-1)\,\mu^{p-1},
$$

with $1<p<2$. Shape $\alpha$ is large when $p$ is near 1, so sizes are almost constant and the total behaves like a scaled Poisson; $\alpha$ is small when $p$ is near 2, so many tiny sizes add up to something close to a gamma.

<figure class="mfig ds-figure ds-wide" id="ds-tweedie" aria-label="The compound Poisson-gamma density and its probability of exactly zero, for a chosen power">
<div class="ds-plot"></div>
<figcaption><b>A lump at zero and a curve beyond it.</b> The bar at the left is the chance that no event arrived; the curve is the density of the total when some did. The mean is 1.5 and the dispersion 2 throughout. Drag the handle to change the power $p$.</figcaption>
</figure>

With these settings the chance of an exact zero is $0.449$ at $p=1.1$, $0.294$ at $p=1.5$ and $0.0055$ at $p=1.9$: the zeros disappear as $p$ moves toward 2. That is a consequence of the construction, not an extra choice. The mean, dispersion and power fix the rate $\lambda$, and the rate fixes the probability of zero as $e^{-\lambda}$. A Tweedie model cannot be told that zeros are more or less common than that formula says, while the positive part stays the same.

When zeros have their own cause, the other story is a gate. A customer is a never-buyer with probability $g$, and otherwise orders as a Poisson with mean $\lambda$. Marginalising over the gate gives the **zero-inflated Poisson**: $\Pr(0)=g+(1-g)e^{-\lambda}$, and $(1-g)\lambda^k e^{-\lambda}/k!$ for $k\ge1$. With $g=0.3$ and $\lambda=2$, the chance of zero is $0.395$, where a plain Poisson with the same $\lambda$ gives $0.135$. The gate sets the number of zeros on its own, which the Tweedie cannot.

<details class="answer"><summary>Worked example: spend per customer, with its zeros</summary>

Take a Tweedie with mean $1.5$, dispersion 2 and power $p=1.5$ (in tens of dollars per month). The rate is $\lambda=1.22$ orders, so a share $e^{-1.22}=0.294$ of customers spend exactly nothing. The mean among those who do spend is $1.5/(1-0.294)=2.12$, larger than the overall mean because the zeros are left out. The chance of spending more than 5 is $0.062$. For comparison, if 30% of the base were never-buyers and the rest ordered as a Poisson with mean 2, the zero-inflated model gives $0.395$ zeros, where a plain Poisson with that mean would give $0.135$.

</details>

<details class="more"><summary>Derivation: why the Tweedie variance is a power of the mean</summary>

By the law of total variance, $\operatorname{Var}Y=\mathbb E[N]\operatorname{Var}X+\operatorname{Var}N\,(\mathbb EX)^2=\lambda\alpha\theta^2+\lambda\alpha^2\theta^2=\lambda\alpha(\alpha+1)\theta^2$. With $\alpha=(2-p)/(p-1)$ we have $\alpha+1=1/(p-1)$, so $\operatorname{Var}Y=\lambda\alpha\theta\cdot\theta/(p-1)=\mu\theta/(p-1)$. Substituting $\theta=\phi(p-1)\mu^{p-1}$ gives $\operatorname{Var}Y=\phi\mu^p$.

</details>

## When the clock has a memory: the Weibull

**Teaser.** Do subscribers cancel at a steady rate, or does the risk of cancelling change with how long they have been subscribed?

The exponential's constant chance per tick is assumption 3. Relax it: let the chance of an event in the next tick, given that none has happened yet, depend on the time $t$ that has passed. Call it the **hazard** $h(t)$. The survival function is then the product of "no event" over all ticks so far, which in the limit is

$$
\Pr(T>t)=\exp\!\Big(-\!\int_0^t h(s)\,ds\Big).
$$

A constant hazard gives back the exponential. Relaxing assumption 3 is the third move, alongside composing and marginalising. The simplest change is a hazard that is a power of time, $h(t)=k\lambda(\lambda t)^{k-1}$, which gives the **Weibull** law, $\Pr(T>t)=e^{-(\lambda t)^k}$. The shape $k$ reads directly as a story about risk: with $k<1$ the hazard falls (the longer you have survived, the safer you are), with $k=1$ it is constant, and with $k>1$ it rises (things wear out).

<figure class="mfig ds-figure ds-wide" id="ds-weibull" aria-label="The hazard and the density of the Weibull for a chosen shape">
<div class="ds-plot"></div>
<figcaption><b>The shape is the story of risk.</b> Top: the rate at which people cancel right now, among those still subscribed. Bottom: the resulting waiting-time density. Drag the handle to change the shape $k$; at $k=1$ the hazard is flat and the density is the exponential.</figcaption>
</figure>

For the shop, a fitted shape below 1 means most cancellations happen early and loyalty accumulates; a shape above 1 means tenure raises the risk, which is the case for contracts that expire. Which it is cannot be settled by thought: it is a number to estimate from the cancellation times, and a shape near 1 says the exponential was enough.

<details class="answer"><summary>Worked example: does cancellation risk fall with tenure?</summary>

Suppose cancellation times fit a Weibull with shape $k=0.7$ and rate $\lambda=0.1$ per month. The hazard is $h(t)=k\lambda(\lambda t)^{k-1}$. In the first month it is $0.14$ per month; at twelve months it is $0.066$. Early risk is $2.1$ times later risk, as a shape below 1 promises. The chance of still being subscribed at twelve months is $e^{-(\lambda\cdot 12)^{k}}=e^{-1.2^{0.7}}=0.321$. A constant-hazard model with the same survival at twelve months would put the hazard at $0.095$ per month throughout, overstating the risk for veterans and understating it for newcomers.

</details>

## The family tree

Every distribution above is one move from another, and the moves are the three named at the start. Click a distribution in the map to see how it is built, and the shop question that goes with it.

<figure class="mfig ds-figure ds-wide ds-graph" id="ds-graph" aria-label="A map of how the distributions are built from one another">
<div class="ds-plot"></div>
<figcaption><b>Everything starts at a coin and a uniform draw.</b> Solid black arrows compose: add, wait, multiply, take a limit; the Weibull arrow relaxes an assumption. Copper arrows marginalise: average over a hidden quantity. Columns are only layout. Click a distribution, or focus the map and use the arrow keys, to light up how it is built and where it leads.</figcaption>
</figure>

## What the story cannot tell you

A histogram does not say which story produced it. The negative binomial is the answer to a Poisson with a wandering rate, and it is also the answer to events that arrive in clumps, and the two histograms are the same. When a count is overdispersed, the negative binomial fits it, but fitting it is not evidence that a rate was hidden.

The assumptions can also fail. If the chance per tick changes with the time of day, the process is no longer stationary, the gaps are no longer exponential, and a single $\lambda$ is a poor summary; the Weibull handles a change with tenure, not a change with the clock. If one ticket raises the chance of the next, ticks are not independent, and the counts are overdispersed for a reason a Poisson with a better rate will not repair. The normal needs many small, nearly independent terms; a revenue day dominated by one huge order is not that. The Tweedie assumes the sizes come from a gamma and are independent of the count, and it ties the zeros to the rest, which the zero-inflated model does not.

The opening asked what a 12% conversion in a small test means: with 80 visits the true rate sits somewhere between about 10% and 23%, which is a beta, not a coin. Thirty-one buyers out of 200 was a one-in-twelve event, a quiet ticket desk owed nothing, and a customer's spend came with a heavy right tail and, often, a lump at zero. Each answer was a different construction from the same ticks, and the useful habit for the next question is to ask what the process is, what is hidden, and which move the answer needs.

## Appendix: an A/B test uses four of these distributions

**Teaser.** The shop converts 12% of visitors on the current page and believes a redesign lifts that to 13.5%. How many visitors does it need, and when may it call the result?

Each arm of the test is a row of coin flips, one per visitor, so the number of buyers among $n$ visitors is binomial. The true rate behind it is hidden, and after $x$ buyers its posterior is a $\text{Beta}(x+1,\,n-x+1)$. Two readings of the same data follow.

- **The belief reading.** Put the two betas side by side and integrate: the chance that $B$'s rate exceeds $A$'s is $\int f_B(r)\,F_A(r)\,dr$, where $F_A$ is the cumulative distribution of $A$'s beta.
- **The surprise reading.** For large $n$ each observed rate is approximately normal, so their difference is too, with standard error $\sqrt{\bar p(1-\bar p)\cdot 2/n}$, where $\bar p$ is the pooled rate. The distance of the gap from zero, in standard errors, is $z$, and the $p$-value is the chance of a gap at least that far from zero if the two arms were identical.

<figure class="mfig ds-figure ds-wide" id="ds-ab" aria-label="Two beta posteriors for the conversion rates of arm A and arm B, with the chance that B is better and the z-test result">
<div class="ds-plot"></div>
<figcaption><b>Same data, two readings.</b> Arm A converts 12% and arm B converts 13.5% at every size. The curves are the beta posteriors for each arm's true rate; the lines below give the chance that B beats A and the z-test. Drag the handle to change the visitors per arm.</figcaption>
</figure>

The two numbers answer different questions. At 2,000 visitors per arm, B leads by 1.5 points, and the chance that B's true rate is higher is 92%. The $p$-value is 0.155, because a gap that size appears in about one test in six between identical pages. The first sentence is about the rate; the second is about how unusual the data would be under no difference. Neither is wrong, and neither alone says what to do. The next section asks how long the test must run before the second one clears the usual bar.

<details class="answer"><summary>Worked example: how many visitors does the test need?</summary>

To detect a move from 12% to 13.5% with a 5% false-win rate ($z=1.96$) and an 80% chance of catching the real lift ($z=0.84$), each arm needs $n=2(1.96+0.84)^2\,\bar p(1-\bar p)/\delta^2$ visitors, where $\delta=0.015$ and $\bar p=0.1275$. That is about 7,800 per arm. With 1,000 visitors a day split evenly, each arm gets 500 a day, so the test needs about 16 days. Halving the lift to detect, from 1.5 points to 0.75, roughly quadruples the visitors to about 30,000 per arm.

</details>

<details class="more"><summary>Derivation: where the sample size comes from</summary>

The gap between the observed rates has standard error $s=\sqrt{2\bar p(1-\bar p)/n}$. To call a gap real at the 5% level, the observed gap must exceed $1.96\,s$. If the true gap is $\delta$, the observed gap is approximately normal around $\delta$ with spread $s$, and it exceeds $1.96\,s$ with probability 80% when $\delta=(1.96+0.84)\,s$. Solving $\delta=2.8\sqrt{2\bar p(1-\bar p)/n}$ for $n$ gives the formula. The factor $1/\delta^2$ is why small lifts are expensive.

</details>

### Looking every day is a different test

**Teaser.** The dashboard updates daily. On day six the $p$-value dips below 0.05. Should the shop stop and ship?

The 5% promise holds for one look at a size chosen in advance. Over time the $z$ value of two identical pages wanders like a random walk, and each extra look is another chance for the walk to cross $\pm1.96$ by luck.

<figure class="mfig ds-figure ds-wide" id="ds-peek" aria-label="Forty simulated tests of two identical pages, and the share declared a winner for different numbers of looks">
<div class="ds-plot"></div>
<figcaption><b>Every look is another chance.</b> Left: forty tests in which both pages convert 12% (the dashed lines are the 5% cut-offs); the dark paths crossed the line at one of the looks, marked on the axis, and stopped being honest. Right: of 2,000 such tests, the share declared a winner for each number of looks, against the promised 5% (dashed). Drag the handle to change how often you look.</figcaption>
</figure>

With one look the share is close to the promised 5%. With five looks it is about 15%, with ten about 20%, and with twenty about 26%: one test in four ships a page that does nothing. A threshold on the chance that $B$ is better is not exempt, because a stopping rule of that kind also rewards luck, and its error rate is worth checking by simulation. The remedies are to fix the visitors per arm in advance and read the result once, or to use a sequential design that raises the bar at each look so the overall promise holds.

<details class="answer"><summary>Worked example: what does daily peeking cost across many tests?</summary>

A team runs 100 tests where the new page does nothing. Reading each once at the planned size, about 5 of them show a false win. Checking twenty times and stopping at the first crossing, about 26 do. That is 21 extra pages shipped for no gain, and each one also carries the redesign's maintenance cost.

</details>

### Counts and money: change the process, keep the recipe

**Teaser.** The shop also cares about orders per visitor and revenue per visitor. Do these need a different test?

They need a different likelihood and the same recipe. A visitor who can order several times produces a count, so the orders in $n$ visitors are Poisson with a hidden rate, and a gamma belief about that rate becomes a $\text{Gamma}(a+\text{orders},\;b+n)$ belief after the data. Comparing the two arms is again a question about two hidden rates. Revenue per visitor has a lump at zero and a skewed rest, which is the Tweedie of the previous sections, and its average moves slowly, so the test needs many more visitors than a conversion test with the same relative lift. In each case the questions are the ones the essay keeps asking: what is the process, what is hidden, and which move turns the data into a belief about it.

## Further reading

- [Blitzstein and Hwang — *Introduction to Probability*](https://projects.iq.harvard.edu/stat110/home). The distributions told as stories; I did not check the edition or the specific chapters here.
- [The Tweedie distribution](https://aggregate.readthedocs.io/en/latest/5_technical_guides/5_x_tweedie.html), a technical guide from the `aggregate` documentation: the compound Poisson–gamma form for $1<p<2$ with its point mass at zero.
- [SAS — The GENMOD procedure, Tweedie distribution](https://support.sas.com/documentation/cdl/en/statug/68162/HTML/default/statug_genmod_details28.htm). The map from the Tweedie parameters to the Poisson rate, gamma shape and gamma scale.
- Smyth — *Regression analysis of quantity data with exact zeros*, 1996. The standard reference for Tweedie regression with exact zeros; I confirmed that it exists and where it was published but did not read it.