Probability

1. Demand is a random variable

A random variable assigns a number to each possible outcome of an uncertain situation. Next week's demand $D$ is a random variable: before the week begins, we do not know its value. Its distribution describes which values are possible and how likely they are. In some weeks demand may be 150 pairs, in others 340. A probability model describes this uncertainty. It does not predict the exact value of the next observation.

Outcomes, events and probabilities

An ice-cream vendor is planning tomorrow's sales. An outcome is one complete possible description of what happens. An event is a collection of outcomes, such as “it rains.” A probability is a number from 0 to 1 measuring how likely an event is. In this simple model, tomorrow is either sunny or rainy: the two outcomes are mutually exclusive (they cannot both happen) and exhaustive (one must happen). Their probabilities therefore add to 1.

With two sources of uncertainty, an outcome specifies both stock availability and weather. The stockout happens with probability $0.05$, and no stockout happens with probability $0.95$. For this example we assume independence: learning whether chocolate is available does not change the weather probabilities. For events $A$ and $B$,

The symbol $\cap$ means “and.” The expression $P(B\mid A)$ means “the probability of $B$, given that $A$ happened.” Two events are independent if and only if $P(A\cap B)=P(A)\,P(B)$. Between symbols, the multiplication sign is often omitted: $P(A)P(B)$ means $P(A)\times P(B)$. Between numbers, we always write it.

To find the probability of a joint outcome, multiply the probabilities along its path. To find the probability of an event that contains several outcomes, add the probabilities of its mutually exclusive paths. For example, with cloudy weather included, the probability of a stockout and sun is

Sun occurs on two mutually exclusive paths, one with chocolate available and one with a stockout. Adding the two path probabilities gives

If stockouts were more likely on sunny days, we would use the relevant conditional probabilities on the branches. Multiplying the two unconditional probabilities would then be unjustified.

A random variable and its realisations

For Amélie's sunglasses, we also distinguish uncertain demand from an observed value. We write random variables in capital letters: $D$ is next week's demand, a quantity not yet known. Once a week has passed, the value we observe is called a realisation of the random variable, and we write it in lower case. Amélie's spreadsheet contains the 52 realisations $d_1 = 108,\ d_2 = 112,\ \dots,\ d_{52} = 112$. The random variable $D$ has a distribution, and each realisation $d_t$ is a number drawn from it. Every probability statement in this course, such as $P(D\le250)$ or $E[D]$, is about the capital letter, and every dataset is a collection of lower-case values.

A random variable can be discrete (it can only take specific values, like the number of customers who walk in: 0, 1, 2, …) or continuous (it can take any value in a range, like the weight of a shipment). Weekly demand for whole pairs of sunglasses is discrete. Its observed values range from 108 to 367, so many values are possible, and we can approximate it by a continuous random variable with a smooth density curve.

Sample statistics estimate population parameters

We assume that the random variable $D$ has a true distribution, called the population distribution. Population parameters summarise it, for example the expectation $\mu = E[D]$ and the standard deviation $\sigma$. In practice we observe only a sample, and a finite sample generally cannot tell us the exact values of these parameters. We therefore compute sample counterparts and use them as estimates.

From her 52 realisations, Amélie computes two sample statistics: the sample mean $\bar d$ and the sample standard deviation $s_d$: $$\begin{aligned} \bar d &= \frac{1}{52}\sum_{t=1}^{52} d_t = 230.5\text{ pairs}, \\[8pt] s_d &= \sqrt{\frac{1}{51}\sum_{t=1}^{52}(d_t-\bar d)^2} \\ &= 74.8\text{ pairs}. \end{aligned}$$ These numbers are computed from the data. A different sample of 52 weeks would give slightly different values. They are distinct from the population mean $\mu = E[D]$ and the population standard deviation $\sigma$ of the distribution of $D$, which are properties of the random variable itself. Section 2 makes this distinction precise.

To compute probabilities, we need a model with explicit assumptions. The histogram is roughly bell-shaped and symmetric. We therefore assume that weekly demand follows a normal distribution, $$D \sim N(\mu, \sigma^2),$$ with mean $\mu$ and variance $\sigma^2$. We estimate the unknown parameters with the sample statistics: $\hat\mu = \bar d = 230.5$ and $\hat\sigma = s_d = 74.8$. The hat denotes an estimate. The histogram describes the observed weeks. The normal curve is a model that we choose because it fits the data reasonably well and is easy to compute with. Every later calculation in the course that uses this model, such as service levels, safety stocks and newsvendor quantities, depends on this assumption. For that reason, we check the model against the data, as in the figure above.

The normal distribution

The normal distribution is a continuous distribution that approximates many real-world distributions. One way to describe it is through its probability density function, or density, $f(x)$ (Section 4 explains densities in more detail):

You do not need to memorise this density. Remember its shape and its two parameters: $\mu$ locates the centre, and $\sigma$ controls the spread. Increasing $\mu$ shifts the curve to the right without changing its spread. Increasing $\sigma$ makes the curve wider and lower, and its total area remains 1. Areas under the curve give probabilities, not the height at a point. We write $N(\mu,\sigma^2)$, with the variance as the second parameter.

The standard normal has mean 0 and standard deviation 1. Subtracting the mean and dividing by the standard deviation converts any normal variable to this scale:

The variable $Z$ has mean 0 and standard deviation 1. A value $Z=2$ means “two standard deviations above the mean.” We write $\phi$ for the standard normal density and $\Phi$ for its cumulative distribution function, introduced in Section 4. A normal model allows negative values, so it is only an approximation for demand, which cannot be negative. Here the model probability of negative demand is about 0.1%.

2. Expectation: a probability-weighted average

That number is the expectation (or mean) $E[X]$, the probability-weighted average of all possible values. For a discrete random variable, the probability mass function (pmf) gives the probability $P(X=x)$ of each possible value $x$, and these probabilities sum to 1. To compute the expectation, multiply each value by its probability, then add: $$E[X] = \sum_x x \, P(X = x).$$ Here $x\,P(X=x)$ means $x\times P(X=x)$. With the probabilities $0.05$, $0.15$, $0.35$, $0.30$ and $0.15$ for the values 0 to 4, which are the initial values of the sliders below,

For a continuous random variable with density $f(x)$, an integral adds up the contributions over a continuous range, and the density supplies the weights: $$E[X] = \int_{-\infty}^{\infty} x\, f(x)\, dx.$$ The expectation is the balance point of the distribution: a cardboard cut-out of the bar chart would balance at $E[X]$. It is also the long-run average. For independent repeated draws from the same distribution with a finite expected absolute value, the observed average converges to $E[X]$. The expectation need not be a possible outcome. A mean of 2.35 units per day does not mean that any day has a fractional order. The expectation $E[X]$ is defined by the distribution alone. To compute it, we need to know the probabilities. Here we know them because we set them with the sliders below.

A sample mean is itself random

We now connect the expectation to data, using the distinction from Section 1. Suppose we do not set the probabilities but instead observe a random sample: $n$ draws $X_1, X_2, \dots, X_n$ from the distribution that are independent and identically distributed (abbreviated i.i.d.). From the sample we compute the sample mean $$\bar X_n = \frac{1}{n}\sum_{i=1}^{n} X_i.$$ Before the draws are observed, $\bar X_n$ is a random variable. If we observe the process for $n$ days twice, the two sample means will usually differ. The law of large numbers connects the sample mean to the expectation. For i.i.d. draws with a finite expected absolute value, as $n \to \infty$, $$\bar X_n \;\longrightarrow\; E[X].$$ In other words, as the sample grows, the probability that the sample mean is far from the expectation tends to zero. For any fixed number of days $n$, however, the observed sample mean is only an estimate of $E[X]$ and contains sampling variation. Amélie's 52-week sample mean $\bar d = 230.5$ from Section 1 is such an estimate of the unknown population mean $\mu = E[D]$. A different year of data would give a different estimate.

Linearity of expectation

Expectation is linear: for constants $a$ and $b$, $$E[aX + b] = a\,E[X] + b.$$ For example, if each unit sells for a fixed price $p$, the expected revenue from selling all demand is $E[pD]=p\,E[D]$. Multiplying demand by a constant multiplies its expectation by the same constant, and adding a constant adds it to the expectation.

More generally, expectation distributes over sums. For random variables with finite expected absolute values and constants $a,b$,

Independence is not required. Expected total demand is the sum of expected demands even when the demands move together.

3. Variance and standard deviation: measuring spread

The expectation $E[D]$ gives the centre of a distribution, but it does not say how far a typical outcome lies from that centre. For that we use the variance, the expected squared deviation from the mean: $$Var(D) = E\big[(D - \mu)^2\big].$$ We square the deviations so that positive and negative deviations do not cancel out. Without squaring, the expected deviation is always zero, because by linearity $E[D-\mu]=E[D]-\mu=0$. Squaring changes the units: the variance of demand measured in pairs is measured in pairs². To return to the original units, we take the square root, which gives the standard deviation: $$\sigma=\sqrt{Var(D)}.$$ The standard deviation measures the typical scale of deviations. It is not the same as their average absolute size. As with the mean, we keep the population and sample versions apart. $Var(D)$ and $\sigma$ are properties of the distribution. The value $s_d = 74.8$ pairs from Section 1 is a sample standard deviation, computed from the 52 observations and used as an estimate $\hat\sigma$ of the model parameter $\sigma$.

Estimating spread from a sample

For $n\ge2$ observations $x_1,\ldots,x_n$, first calculate their sample mean. Then measure each observation's distance from that mean:

Divide by $n-1$, not $n$, when estimating variance from an i.i.d. sample. We have estimated the mean from these same observations, so the deviations add to zero. Knowing $n-1$ of them therefore determines the last one. Dividing by $n$ would underestimate the population variance on average. The $n-1$ correction makes $s^2$ an unbiased estimator of the variance: for an i.i.d. sample with finite variance, $E[s^2]=\sigma^2$.

The square root $s$ is the usual sample standard deviation and is in the observations' original units. Strictly speaking, taking the square root does not preserve unbiasedness, so $s$ itself is not generally unbiased for $\sigma$. For Amélie's 52 weeks, $n-1=51$, which gives the $s_d=74.8$ pairs used above.

Shifting and scaling a random variable

Variance follows its own rule, and it is not linear like expectation: $$Var(aX + b) = a^2\, Var(X).$$ Adding a constant $b$ shifts the distribution to the left or right without changing its shape, so the variance does not change. Multiplying by $a$ multiplies the standard deviation by $|a|$ and the variance by $a^2$. To see why, write $\mu=E[X]$. By linearity, $E[aX+b]=a\mu+b$, so

The fourth line uses linearity of expectation to move the constant $a^2$ outside.

The standard deviation of a sum of independent variables

To calculate the uncertainty in total demand over several periods, we use the following rule. If $X_1,\dots,X_n$ are independent and identically distributed (i.i.d.) with mean $\mu$ and standard deviation $\sigma$, their sum has mean $n\mu$ and standard deviation $$SD(X_1 + X_2 + \cdots + X_n) = \sigma\sqrt{n}.$$ The mean of the sum grows in proportion to $n$, while its standard deviation grows in proportion to $\sqrt n$, because independent deviations partly cancel out. In Section 5, we use this rule to calculate the uncertainty in demand during a delivery lead time.

Averages become more precise with a larger sample size

The sum grows with the number of observations, and the average divides that sum by $n$. Consider an i.i.d. sample with mean $\mu$ and finite variance $\sigma^2$. By linearity of expectation, the expected sample mean is

For the variance, the rule $Var(aX)=a^2Var(X)$ with $a=1/n$ and the sum rule above give

Taking the square root gives the standard deviation of the sample mean:

The standard deviation of a sample mean is called its standard error. With four times as many independent observations, the standard error is halved, because $\sqrt{4n}=2\sqrt n$. Individual observations remain just as variable, but their average is more precise. Using the sunglasses sample standard deviation gives an estimated standard error of $74.8/\sqrt{52}\approx10.4$ pairs, if the weeks can be treated as i.i.d.

4. The CDF: cumulative probabilities

For a continuous random variable, the smooth idealisation of the histogram is the probability density function (pdf), written $f(x)$. The density describes the probability close to a point: for a small width $\delta$, $$P(x \le D \le x + \delta) \approx f(x)\,\delta.$$ The probability of an interval is the area under the curve over that interval: $$P(a \le D \le b) = \int_a^b f(x)\,dx.$$ The value $f(x)$ itself is a density, not a probability, and it can exceed 1 for narrow distributions. In a continuous model, a single exact value has probability zero, and probabilities apply to intervals. Actual demand is a whole number of pairs, so a single value can have positive probability. The smooth model is an approximation.

The pdf shows where the probability is concentrated. Amélie's question, however, is cumulative: for a given stock level, she wants the probability that she has enough stock. The cumulative distribution function (CDF) gives this probability for every stock level $x$: $$F(x) = P(D \le x).$$ The CDF never decreases, which means it is weakly increasing. It approaches 0 for small $x$ and approaches 1 for large $x$. The pdf and the CDF describe the same distribution in two different ways. The pdf gives the probability close to each point, and the CDF gives the accumulated probability at or below each value. Either one determines the other: $$F(x) = \int_{-\infty}^{x} f(t)\,dt, \qquad f(x) = F'(x).$$ The CDF is the accumulated area under the pdf, and the pdf is the slope of the CDF where the derivative exists.

In the figure below, the smooth curves come from the model we assumed in Section 1, $D \sim N(\mu, \sigma^2)$ with the estimates $\hat\mu = 230.5$ and $\hat\sigma = 74.8$. For this model, $F(x)=\Phi\big((x-\mu)/\sigma\big)$. The step-shaped curve is the empirical CDF, computed directly from the 52 observed weeks without a distributional assumption. At each $x$, it is the fraction of weeks with demand at or below $x$.

The CDF gives the stockout risk for each stock level. If Amélie stocks $x$ pairs, $F(x)$ is the probability that she has enough stock, and $1-F(x)$ is the probability of a stockout. For example, the model gives $F(250)\approx0.603$, so stocking 250 pairs leads to a stockout with probability of about $1-0.603=0.397$. In the data, 31 of the 52 weeks had demand of at most 250 pairs, and $31/52\approx0.596$.

5. Quantiles and service levels: reading the CDF backwards

The answer is a quantile. For a target probability $\alpha$, the $\alpha$-quantile $F^{-1}(\alpha)$ is the smallest value $x$ at which the CDF reaches at least $\alpha$: $$F^{-1}(\alpha)=\inf\{x:F(x)\ge\alpha\},\qquad 0<\alpha<1.$$ Here $\inf$ (infimum) denotes the smallest value of the set. In supply chain applications, we also call $\alpha$ the service level. For the continuous and strictly increasing normal CDF, $P(D\le F^{-1}(\alpha))=\alpha$. For discrete demand, the CDF jumps, so the covered probability can be greater than $\alpha$. For example, with the initial spare-part probabilities of Section 2, $F(2)=0.55$ and $F(3)=0.85$, so $F^{-1}(0.80)=3$, and stocking 3 units covers demand with probability $0.85$.

The same idea is used in finance. The Value-at-Risk (VaR) at level $\alpha$ is the $\alpha$-quantile of a loss variable, $\text{VaR}_\alpha=F_{\text{loss}}^{-1}(\alpha)$. It is a loss threshold that is exceeded with probability at most $1-\alpha$. A service quantile is a stock threshold that covers demand with probability at least $\alpha$. Neither quantile tells us how large the outcomes beyond the threshold will be.

We apply quantiles to an $(r,Q)$ inventory system, named after its two parameters: the reorder point $r$ and the order quantity $Q$. The firm continuously monitors its inventory position, which is on-hand stock plus on-order stock minus backorders. Whenever the inventory position falls to $r$, the firm places a replenishment order of fixed size $Q$, which arrives after the lead time $L$ (here measured in weeks). Between ordering and receiving, the firm must cover demand from its available stock. In this session we calculate $r$ for a chosen service target. The order quantity $Q$ is a separate decision.

We assume a fixed lead time of $L$ whole weeks and independent weekly demands $D_t\sim N(\mu,\sigma^2)$. The demand during the lead time is the sum $D_L=\sum_{t=1}^{L}D_t$. A sum of independent normal variables is also normal. Its mean follows from linearity of expectation, and its variance from the sum rule of Section 3: $$\begin{aligned} \mu_L &=E[D_L]=\sum_{t=1}^{L}E[D_t]=L\mu, \\[4pt] \sigma_L^2 &=Var(D_L)=\sum_{t=1}^{L}Var(D_t)=L\sigma^2, \\[4pt] \sigma_L &=\sqrt{L\sigma^2}=\sigma\sqrt L, \\[4pt] D_L &\sim N(\mu_L,\sigma_L^2). \end{aligned}$$ Here $\mu_L$ is the product $L\times\mu$. The expected demand over $L$ weeks is $L$ times the expected weekly demand, while the standard deviation grows only by the factor $\sqrt{L}$.

From a normal quantile to the reorder point

A Type-1 service level $\alpha$ is the probability of covering demand during the lead time. We therefore choose $r$ as the $\alpha$-quantile of $D_L$, such that $$P(D_L\le r)=\alpha.$$ Section 1 introduced a useful property of normal variables: subtracting the mean and dividing by the standard deviation gives a standard normal variable. For lead-time demand, this is

We subtract $\mu_L$ from both sides of the inequality $D_L\le r$, then divide both sides by $\sigma_L$. Because $\sigma_L>0$, the inequality keeps its direction. These steps describe the same event on a new scale, so its probability does not change:

We set this probability equal to $\alpha$. The standard normal quantile $z_\alpha=\Phi^{-1}(\alpha)$ is the value such that $P(Z\le z_\alpha)=\alpha$ for $Z\sim N(0,1)$. Applying $\Phi^{-1}$ to both sides gives the standardised threshold $z_\alpha$. We then solve for $r$:

The third line multiplies both sides by $\sigma_L$, and the fourth adds $\mu_L$ to both sides. More generally, the $\alpha$-quantile of $N(\mu,\sigma^2)$ is $\mu+\sigma\Phi^{-1}(\alpha)$: the mean plus the standard normal quantile times the standard deviation.

The reorder point $r$ is the inventory position at which we place an order. The safety stock $ss$ is the stock held above expected lead-time demand to cover demand uncertainty:

For example, with $L=2$ weeks and $\alpha=0.95$, the standard normal quantile is $z_{0.95}\approx1.645$, and

With the exact value of $z_{0.95}$, the model quantile is about 634.998 pairs. For whole pairs and an “at least $\alpha$” target, we round the exact model quantile up, which gives 635 pairs. Rounding intermediate values, such as $z_{0.95}$ to 1.645, can change the last digit, so we keep full precision until the final answer.

With no safety stock, $r=\mu_L$. A normal variable is equally likely to be above or below its mean, so this gives $\alpha=0.5$, and stock runs out in about half of the replenishment cycles.

6. Counting arrivals: the Poisson distribution

The distribution is named after the French mathematician Siméon Denis Poisson. Because poisson means “fish” in French, fishing gives a memory aid: the number of fish caught in one hour. If catches occur independently at a constant average rate of $\lambda$ per hour, a Poisson model describes the number caught. The fishing example is a memory aid, not the historical origin of the name.

So far we have modelled the quantity demanded each week. Strictly speaking, that was also a count (whole pairs of sunglasses), but a large count, in the hundreds. For this reason we approximated it by a continuous random variable with a normal distribution. Other random quantities are small counts, for which that approximation is poor. Examples are the number of orders that arrive in a period, the number of trucks that arrive today, and the number of your 25 SKUs that receive an order this week. For such quantities we model the count directly. We assume that arrivals occur independently of one another and at a constant average rate $\lambda$ per period. Then the number of arrivals $N$ in a fixed period follows the Poisson distribution: $$P(N = k) = \frac{e^{-\lambda}\lambda^k}{k!}, \qquad k = 0, 1, 2, \dots$$ The parameter $\lambda$ is the expected number of arrivals in the period, and $k!$ means $1\times2\times\cdots\times k$, with $0!=1$. The mean and the variance of a Poisson random variable are both equal to $\lambda$: $$E[N]=Var(N)=\lambda,\qquad SD(N)=\sqrt\lambda.$$ For example, with $\lambda=3$,

Comparing the sample mean and the sample variance is a useful first check of a Poisson model. Similar values are consistent with a Poisson model, but they do not prove the model or its assumptions.

We now check the assumptions against data. Snyder & Shen's battery dataset records orders for 25 products over 26 periods. For each period, we count how many of the 25 products received at least one order. Each product counts once, even if it received several orders. This gives 26 counts, each between 0 and 25. A Poisson distribution allows counts above 25, so it can only be an approximation here. If the products ordered independently with the same probability, the count would follow a binomial distribution, and the Poisson distribution approximates it well when the individual ordering probabilities are small.

These tools describe the risk, but they do not choose an order quantity. In Session 2, Amélie chooses an order quantity by trading off the cost of holding too much stock against the cost of holding too little. This is an optimization problem.