Appendix C — Common Probability Distributions
After reading this chapter you should be able to:
- identify the distributions used to model demand and lead time
- compute the loss functions that appear in inventory cost expressions
- explain why a variance of a shortage needs a loss function one order higher
C.1 Discrete Distributions
Demand arrives in whole units, so a discrete distribution is the natural default and a continuous one is an approximation made for convenience. This section gives the four families the book uses, what each assumes, and the single diagnostic that usually decides between them. Section C.3 gives their loss functions.
The diagnostic is the variance to mean ratio. Let \(X\) represent a demand over some interval, and let
\[ \mathit{VMR} = \frac{\mathit{Var}[X]}{E[X]} \tag{C.1}\]
Recall from Section 2.2.6 that a demand history is characterized before it is modeled. Equation C.1 is the characterization that selects among the families below, because each of them fixes the ratio or bounds it on one side. Notice that \(\mathit{VMR}\) is not the variability coefficient of Equation 5.10 or the squared coefficient of variation of Equation 2.3. Those are dimensionless; this one carries the units of \(X\) and is therefore a statement about a particular time bucket.
C.1.1 The Poisson Distribution
\[ g(x) = \frac{\lambda^{x}e^{-\lambda}}{x!}, \quad x = 0, 1, 2, \ldots, \qquad E[X] = \mathit{Var}[X] = \lambda \tag{C.2}\]
The Poisson is what demand looks like when customers arrive independently at a constant average rate and each takes one unit. That is the standard model for a slow-moving part, and it is the distribution to reach for first when the item is demanded singly.
Notice what Equation C.2 fixes. The mean and the variance are the same parameter, so \(\mathit{VMR} = 1\) exactly. Thus, the Poisson is not a family you fit to a variance; it is a family whose variance is decided the moment you estimate its mean. A demand history whose \(\mathit{VMR}\) is far from one is telling you that the Poisson is the wrong family, and no choice of \(\lambda\) will repair it.
C.1.2 The Geometric and Negative Binomial Distributions
The negative binomial is the family to reach for when \(\mathit{VMR} > 1\), which is the usual case for real demand. Let \(\beta = (1-p)/p\). Then
\[ g(x) = \binom{x+r-1}{r-1}p^{r}(1-p)^{x}, \quad x = 0, 1, 2, \ldots, \qquad E[X] = r\beta, \qquad \mathit{Var}[X] = \frac{r\beta}{p} \tag{C.3}\]
so that \(\mathit{VMR} = 1/p > 1\) for every admissible \(p\). The geometric is the special case \(r = 1\), with \(E[X] = \beta\) and \(\mathit{Var}[X] = \beta/p\).
Two readings of Equation C.3 are useful and they are different. Counting failures before the \(r\)th success is the textbook reading and explains the binomial coefficient. The reading that matters here is that the negative binomial is a Poisson whose rate is itself random, drawn from a gamma distribution. That is, it is what you get when demand is Poisson but the rate varies from interval to interval, which is the situation an inventory analyst meets. Thus, the negative binomial is the Poisson with the constant-rate assumption relaxed, and the extra variance is the price of relaxing it.
A warning about \(p\). Two conventions are in circulation. Equation C.3 makes \(p\) the probability of success, so the mean is \(r(1-p)/p\). Much of the fitting literature, including the selection rules of Section C.4, makes \(p\) the probability of failure, so the mean is \(rp/(1-p)\). The two are exchanged by replacing \(p\) with \(1-p\), and a fit that comes out inverted is almost always this and not an error in the data. You should check which convention a routine uses before passing it a parameter.
C.1.3 The Binomial Distribution
\[ g(x) = \binom{n}{x}p^{x}(1-p)^{n-x}, \quad x = 0, 1, \ldots, n, \qquad E[X] = np, \qquad \mathit{Var}[X] = np(1-p) \tag{C.4}\]
Here \(\mathit{VMR} = 1-p < 1\), so the binomial is the family for demand that is less variable than Poisson. That is rarer than the other direction, and when it appears it usually has a structural cause: a fixed population of \(n\) machines each of which may or may not need the part this period, or a contract that caps the quantity. Notice that the support is bounded above, which is the strongest statement any of these families makes and should be believed only when something in the situation actually enforces it.
C.1.4 The Discrete Empirical Distribution
Let \(S = \{x_{1}, \ldots, x_{k}\}\) represent the distinct values observed and \(p_{j}\) their relative frequencies. The empirical distribution is the data used as its own model, with no family assumed.
It is the right choice more often than its simplicity suggests, and Chapter 7 uses it. Thus, when the history is long enough to populate the range and the decision is a single quantile, fitting a family buys nothing and risks a poor fit. Two costs come with it. The support is bounded by what happened to be observed, so a quantile beyond the largest observation cannot be produced at all, and a long tail is exactly where inventory decisions live. And the estimate inherits every irregularity of the sample, including the ones that are noise.
The rule of thumb. Use the empirical distribution when the question is inside the data and a fitted family when it is outside. A service level of 90% on two hundred observations is inside. A service level of 99% is not.
C.2 Continuous Distributions
A continuous distribution is an approximation to a demand that arrives in whole units, and it is made for two reasons. The arithmetic is easier, and once demand is aggregated over a lead time or a season the granularity stops mattering. This section gives the four families the book uses and, for each, the condition under which the approximation is safe. Section C.3 gives their loss functions.
The diagnostic here is the coefficient of variation. Let
\[ c_{X} = \frac{\sqrt{\mathit{Var}[X]}}{E[X]}, \qquad c_{X}^{2} = \frac{\mathit{Var}[X]}{(E[X])^{2}} \tag{C.5}\]
Unlike the variance to mean ratio of Equation C.1, \(c_{X}\) is dimensionless, so it compares across items and across time buckets. Each family below either fixes \(c_{X}\) or restricts it, and that is what selects among them.
C.2.1 The Exponential Distribution
\[ g(x) = \lambda e^{-\lambda x}, \quad x \ge 0, \qquad E[X] = \frac{1}{\lambda}, \qquad \mathit{Var}[X] = \frac{1}{\lambda^{2}} \tag{C.6}\]
Thus, \(c_{X} = 1\) exactly. The exponential is the continuous counterpart of the Poisson: it is the time between Poisson arrivals, and it appears in this book as a lead time rather than as a demand. Like the Poisson it has no free parameter for variability, so it is a hypothesis to be tested rather than a family to be fitted to a spread.
Notice also that it is memoryless. A lead time modeled as exponential is one for which an order already three weeks late is no more likely to arrive tomorrow than one placed today. That is a strong claim about a supplier, and it is usually false.
C.2.2 The Gamma and Erlang Distributions
\[ g(x) = \frac{x^{\alpha-1}e^{-x/\beta}}{\beta^{\alpha}\Gamma(\alpha)}, \quad x \ge 0, \qquad E[X] = \alpha\beta, \qquad \mathit{Var}[X] = \alpha\beta^{2} \tag{C.7}\]
with shape \(\alpha > 0\) and scale \(\beta > 0\). Then \(c_{X} = 1/\sqrt{\alpha}\), so the shape parameter alone sets the variability and the scale sets the size. The gamma is the flexible workhorse of this appendix: it covers every \(c_{X} > 0\), it lives on the non-negative half line where demand lives, and it reduces to the exponential at \(\alpha = 1\).
The Erlang is the gamma with integer \(\alpha\), which is the sum of \(\alpha\) independent exponentials. That buys a closed form for the distribution function and nothing else, and Table C.5 shows the loss functions are the gamma’s unchanged. Its importance here is structural rather than computational: because an Erlang is a sum of exponentials, mixtures of Erlangs can be built to match any mean and variance, which is the procedure of Section C.4.2.
Notice what Equation C.7 does not have. There is no closed form for \(G(x)\), so a quantile requires a numerical routine. That is rarely an obstacle, and it is not an argument for the normal.
C.2.3 The Normal Distribution
The normal is the most used and the least defensible of the four. It is used because Section B.2 makes it the limit of a sum, so a demand aggregated over many periods or many customers tends toward it, and because its loss functions are the tabulated ones of Table C.2.
Its defect is the one every inventory text mentions and few quantify. The support is the whole real line, so a normal demand model assigns positive probability to negative demand, and the amount is governed entirely by \(c_{X}\).
| \(c_{X}\) | \(P\{X < 0\}\) under a normal model |
|---|---|
| 0.2 | 0.0000 |
| 0.4 | 0.0062 |
| 0.6 | 0.0478 |
| 1.0 | 0.1587 |
| 1.5 | 0.2525 |
Table C.1 is the rule in practice. Thus, a normal model is safe for demand whose coefficient of variation is below roughly 0.3, questionable to about 0.5, and wrong above it. Recall from Example 2.4 that intermittent items routinely run past \(c_{X} = 1\), which is precisely where the normal fails, and why Chapter 7 does not reach for it there.
C.2.4 The Lognormal Distribution
If \(\ln X\) is normal with mean \(\mu\) and variance \(\sigma^{2}\), then \(X\) is lognormal. It is the natural repair of the normal’s defect: the support is strictly positive, the shape is right-skewed as demand data usually is, and it accommodates any \(c_{X}\).
The parameterization is the trap. A fitted demand distribution reports the mean and variance of \(X\), while the density is written in terms of the mean and variance of \(\ln X\). The conversion is
\[ \sigma^{2} = \ln\left(\frac{\mathit{Var}[X]}{(E[X])^{2}} + 1\right) = \ln\left(1 + c_{X}^{2}\right), \qquad \mu = \ln(E[X]) - \frac{\sigma^{2}}{2} \tag{C.8}\]
and the lognormal rows of Table C.5 and Table C.6 are written in terms of \(E[X]\) and \(\mathit{Var}[X]\) for that reason. Notice that \(\sigma^{2}\) depends on the data only through \(c_{X}\), so the shape of a lognormal is set by the coefficient of variation as the gamma’s is by \(\alpha\).
Choosing between the gamma and the lognormal is usually not worth a long argument. Both live on the positive half line, both take any \(c_{X}\), and on demand data they rarely separate by enough to change a policy. The lognormal has the heavier right tail, so it gives the larger answer at a high service level, and that difference grows with \(c_{X}\). Where the two disagree enough to matter, the disagreement is the finding, and Section C.4 says what to do about it.
C.3 Loss Functions
Two quantities appear in almost every stochastic inventory cost expression in this book, and neither is an ordinary expectation. How much demand goes unmet when \(b\) units are on hand, and how much stock is left over when \(b\) units were provided. Both are expectations of a truncated difference, and the functions that give them are called loss functions.
Let \(X\) represent a random variable, usually a demand over some interval. Let \(b\) represent a stock level the demand is compared against. The notation \((X-b)^{+} = \max(X-b, 0)\) is used throughout, so \((X-b)^{+}\) is the amount by which demand exceeds \(b\), and zero when it does not.
Four functions carry almost all of the work, and a fifth is needed once, in Section C.3.6. Let \(G(b) = P\{X \le b\}\) represent the cumulative distribution function of \(X\). Then
\[ G^{0}(b) = 1 - G(b), \qquad G^{1}(b) = E\left[(X-b)^{+}\right] \tag{C.9}\]
define the complementary distribution function and the first order loss function. The second order loss function is defined differently for the two cases, and the difference is not cosmetic:
\[ G^{2}(b) = \begin{cases} \tfrac{1}{2}E\left[(X-b)^{+}(X-b-1)^{+}\right] & X \text{ discrete} \\[4pt] \tfrac{1}{2}E\left[\left((X-b)^{+}\right)^{2}\right] & X \text{ continuous} \end{cases} \tag{C.10}\]
Read \(G^{0}\) carefully in the discrete case. Because \(G(b) = P\{X \le b\}\), the complement is \(G^{0}(b) = P\{X > b\}\), which is \(P\{X \ge b+1\}\) and not \(P\{X \ge b\}\). Every formula in this section depends on that convention, and using \(P\{X \ge b\}\) instead produces answers that are correct only where the two happen to coincide. You should check which convention a source uses before taking a formula from it.
The first order function is the one with the obvious meaning, and the second order function has a use of its own. Chapter 8 needs the expected backorder level of a policy that orders \(Q\) at a time, and that level is a difference of second order loss functions evaluated at \(r\) and \(r+Q\). Chapter 7 needs only the first. The variance of that backorder level needs one order higher again, which is what Section C.3.6 supplies.
C.3.1 Discrete Random Variables
Let \(g(x) = P\{X = x\}\) represent the probability mass function. Writing the expectation as a sum over the range gives
\[ G^{1}(b) = \sum_{x > b}(x-b)g(x) = \sum_{x \ge b}G^{0}(x) \tag{C.11}\]
\[ G^{2}(b) = \tfrac{1}{2}\sum_{x > b}(x-b)(x-b-1)g(x) = \sum_{x > b}G^{1}(x) \tag{C.12}\]
The second form of each is the one to notice. A loss function is a tail sum of the function one order below it, which is what makes both computable by a single downward pass over the support.
For a non-negative \(X\) the pass can start from a moment. Notice that \(G^{1}(0) = E[X]\) and that \(G^{2}(0) = \tfrac{1}{2}E[X(X-1)]\), the latter being the second binomial moment. Thus, for \(b \ge 0\),
\[ G^{1}(b) = E[X] - \sum_{0 \le x < b}G^{0}(x), \qquad G^{2}(b) = \tfrac{1}{2}E[X(X-1)] - \sum_{0 < x \le b}G^{1}(x) \tag{C.13}\]
Equation C.13 is how these are computed in practice. That is, one accumulating loop produces \(G^{0}\), \(G^{1}\) and \(G^{2}\) at every \(b\) from a single evaluation of the mass function at each point.
C.3.2 The Discrete Empirical Distribution
A demand history that is used directly, without fitting, is a discrete empirical distribution. Let \(S = \{x_{1}, x_{2}, \ldots, x_{k}\}\) represent its ordered set of distinct values and let \(g(x_{j}) = p_{j}\) represent the mass at each. Inside the support the sums of Equation C.11 apply unchanged. Outside it they need care, because the support is finite in both directions.
\[ G^{1}(b) = \begin{cases} E[X] - b & b \le x_{1} \\ \sum_{x > b,\; x \in S}(x-b)g(x) & x_{1} < b < x_{k} \\ 0 & b \ge x_{k} \end{cases} \tag{C.14}\]
Notice the first line. Below the smallest possible value every unit of demand exceeds \(b\), so the truncation never binds and \(G^{1}(b)\) is the plain expectation \(E[X] - b\). It is not constant at \(G^{1}(x_{1})\), which is the value it takes only at \(b = x_{1}\) itself. The same reasoning gives
\[ G^{2}(b) = \tfrac{1}{2}\left(E\left[X^{2}\right] - (2b+1)E[X] + b(b+1)\right) \qquad\text{for } b \le x_{1} \tag{C.15}\]
and zero for \(b \ge x_{k}\). For example, take \(S = \{2, 5, 9, 14\}\) with masses \(0.20\), \(0.30\), \(0.35\) and \(0.15\), so that \(E[X] = 7.15\) units. At \(b = 0\) the first order loss function is \(7.15 - 0 = 7.15\) units, not the \(7.15 - 2 = 5.15\) that treating it as constant below the support would give. The two agree at \(b = 2\) and nowhere else.
C.3.3 Continuous Random Variables
Let \(g(x)\) represent the probability density function. The sums become integrals and the tail relations survive:
\[ G^{1}(b) = \int_{b}^{\infty}(x-b)g(x)\,dx = \int_{b}^{\infty}G^{0}(x)\,dx, \qquad G^{2}(b) = \tfrac{1}{2}\int_{b}^{\infty}(x-b)^{2}g(x)\,dx = \int_{b}^{\infty}G^{1}(x)\,dx \tag{C.16}\]
For a non-negative \(X\) and \(b \le 0\) the truncation again never binds, so \(G^{1}(b) = E[X] - b\) and \(G^{2}(b) = \tfrac{1}{2}\left(E[X^{2}] - 2bE[X] + b^{2}\right)\). Thus, \(G^{1}(0) = E[X]\) and \(G^{2}(0) = \tfrac{1}{2}E[X^{2}]\).
The closed forms for the continuous families are all obtained the same way, through the partial expectation functions. Let
\[ H^{1}(b) = \int_{b}^{\infty}x\,g(x)\,dx, \qquad H^{2}(b) = \int_{b}^{\infty}x^{2}g(x)\,dx \tag{C.17}\]
Expanding \((x-b)\) and \((x-b)^{2}\) inside the integrals of Equation C.16 and collecting terms gives
\[ G^{1}(b) = H^{1}(b) - bG^{0}(b), \qquad G^{2}(b) = \tfrac{1}{2}\left(H^{2}(b) - 2bH^{1}(b) + b^{2}G^{0}(b)\right) \tag{C.18}\]
That is, a loss function for a new distribution follows from its partial expectations, and those are often the easier integrals. The lognormal entries of Table C.5 and Table C.6 are derived exactly this way.
C.3.4 The Standard Normal Case
The normal distribution is used often enough that its loss functions carry their own notation. Let \(\phi(z)\) and \(\Phi(z)\) represent the standard normal density and distribution function, and let \(\Phi^{0}(z) = 1 - \Phi(z)\). Then
\[ \Phi^{1}(z) = \phi(z) - z\Phi^{0}(z), \qquad \Phi^{2}(z) = \tfrac{1}{2}\left[\left(z^{2}+1\right)\Phi^{0}(z) - z\phi(z)\right] \tag{C.19}\]
and for \(X \sim N(\mu, \sigma^{2})\) with \(z = (x-\mu)/\sigma\),
\[ G^{0}(x) = \Phi^{0}(z), \qquad G^{1}(x) = \sigma\Phi^{1}(z), \qquad G^{2}(x) = \sigma^{2}\Phi^{2}(z) \tag{C.20}\]
Notice the units. \(\Phi^{1}\) and \(\Phi^{2}\) are dimensionless, and the scaling by \(\sigma\) and \(\sigma^{2}\) is what puts \(G^{1}\) in units and \(G^{2}\) in units squared. A loss function computed without that scaling is one of the more common errors in this material, and the units are what catch it.
Table C.2 gives the three functions at the values of \(z\) that arise most often.
| \(z\) | \(\phi(z)\) | \(\Phi^{0}(z)\) | \(\Phi^{1}(z)\) | \(\Phi^{2}(z)\) | \(\Phi^{3}(z)\) |
|---|---|---|---|---|---|
| \(-2.00\) | 0.0540 | 0.9772 | 2.0085 | 2.4971 | 2.3342 |
| \(-1.00\) | 0.2420 | 0.8413 | 1.0833 | 0.9623 | 0.6819 |
| \(0.00\) | 0.3989 | 0.5000 | 0.3989 | 0.2500 | 0.1330 |
| \(0.50\) | 0.3521 | 0.3085 | 0.1978 | 0.1048 | 0.0485 |
| \(1.00\) | 0.2420 | 0.1587 | 0.0833 | 0.0377 | 0.0152 |
| \(1.28\) | 0.1758 | 0.1003 | 0.0475 | 0.0197 | 0.0074 |
| \(1.50\) | 0.1295 | 0.0668 | 0.0293 | 0.0114 | 0.0041 |
| \(1.645\) | 0.1031 | 0.0500 | 0.0209 | 0.0078 | 0.0027 |
| \(2.00\) | 0.0540 | 0.0228 | 0.0085 | 0.0029 | 0.0009 |
| \(2.33\) | 0.0264 | 0.0099 | 0.0034 | 0.0010 | 0.0003 |
| \(3.00\) | 0.0044 | 0.0013 | 0.0004 | 0.0001 | 0.0000 |
For example, take demand over a lead time as \(N(50, 12^{2})\) units and a stock level of 62 units. Then \(z = (62-50)/12 = 1.00\), so from Table C.2 the probability of a shortage is \(\Phi^{0}(1.00) = 0.1587\) and the expected shortage is \(\sigma\Phi^{1}(1.00) = 12(0.0833) = 0.9996\) units, i.e. about one unit short on average.
That answer is worth checking a second way, and the partial expectation gives one. Reading \(\phi(1.00) = 0.2420\) from the same table,
\[ H^{1}(62) = \mu\Phi^{0}(z) + \sigma\phi(z) = 50(0.1587) + 12(0.2420) = 10.8390 \text{ units} \]
so from Equation C.18, \(G^{1}(62) = 10.8390 - 62(0.1587) = 0.9996\) units. The two routes agree. Carried to full precision both give \(0.999786\) units, and the difference is the four decimal places Table C.2 is printed to.
That check is worth performing whenever a loss function is computed for the first time in a new setting, because it catches a missing \(\sigma\) immediately. Without the scaling the first route returns 0.0833 units while the second still returns 0.9996, and a disagreement by a factor of \(\sigma\) is unmistakable.
C.3.5 Loss Functions for Common Distributions
The closed forms are collected below, one table per order so that each formula has the width of the page to itself. Throughout, \(g\) is the mass or density function and \(G^{0}\) is the complementary distribution function, so every loss function is expressed in terms of quantities the distribution already provides. The parameters and moments are not repeated here; they are in Section C.1 and Section C.2 with the family that owns them.
| Distribution | \(G^{1}(x)\) |
|---|---|
| Poisson, mean \(\lambda\) | \(-(x-\lambda)G^{0}(x) + \lambda g(x)\) |
| Geometric | \(\beta(1-p)^{x}\) |
| Negative binomial, size \(r\) | \(-(x-r\beta)G^{0}(x) + (x+r)\beta g(x)\) |
| Binomial | No simpler form than Equation C.11 |
| Distribution | \(G^{2}(x)\) |
|---|---|
| Poisson, mean \(\lambda\) | \(\tfrac{1}{2}\left\{\left[(x-\lambda)^{2}+x\right]G^{0}(x) - \lambda(x-\lambda)g(x)\right\}\) |
| Geometric | \(\beta^{2}(1-p)^{x}\) |
| Negative binomial, size \(r\) | \(\tfrac{1}{2}\left\{\left[r(r+1)\beta^{2} - 2r\beta x + x(x+1)\right]G^{0}(x) + \left[(r+1)\beta - x\right](x+r)\beta g(x)\right\}\) |
| Binomial | No simpler form than Equation C.12 |
| Distribution | \(G^{1}(x)\) |
|---|---|
| Exponential, rate \(\lambda\) | \(\dfrac{1}{\lambda}e^{-\lambda x}\) |
| Gamma, shape \(\alpha\) and scale \(\beta\) | \(\beta\left[\left(\alpha - \tfrac{x}{\beta}\right)G^{0}(x) + xg(x)\right]\) |
| Erlang | The gamma form, unchanged |
| Normal | \(\sigma\Phi^{1}(z)\), with \(z = (x-\mu)/\sigma\) and \(\Phi^{1}\) from Equation C.19 |
| Lognormal | \(E[X]\,\Phi(\sigma - z) - x\Phi^{0}(z)\), with \(z = (\ln x - \mu)/\sigma\) |
| Distribution | \(G^{2}(x)\) |
|---|---|
| Exponential, rate \(\lambda\) | \(\dfrac{1}{\lambda^{2}}e^{-\lambda x}\) |
| Gamma, shape \(\alpha\) and scale \(\beta\) | \(\tfrac{1}{2}\beta^{2}\left\{\left[\left(\alpha - \tfrac{x}{\beta}\right)^{2} + \alpha\right]G^{0}(x) + \left(\alpha - \tfrac{x}{\beta} + 1\right)xg(x)\right\}\) |
| Erlang | The gamma form, unchanged |
| Normal | \(\sigma^{2}\Phi^{2}(z)\), with \(\Phi^{2}\) from Equation C.19 |
| Lognormal | \(\tfrac{1}{2}\left[\left(\mathit{Var}[X] + (E[X])^{2}\right)\Phi(2\sigma - z) - 2xE[X]\,\Phi(\sigma - z) + x^{2}\Phi^{0}(z)\right]\) |
The lognormal rows deserve a word, because they are the ones most often misquoted. They are Equation C.18 applied to the partial expectations \(H^{1}(x) = E[X]\Phi(\sigma - z)\) and \(H^{2}(x) = \left(\mathit{Var}[X] + (E[X])^{2}\right)\Phi(2\sigma - z)\). It is easy to quote one of those partial expectations as though it were the loss function itself, and the resulting figure is too large by roughly a factor of three at stock levels near the mean.
C.3.6 The Third Order Loss Function
Everything above computes an average. Chapter 8 also needs the variability of a shortage, and that is where a third function enters.
Start with what the first two already give. For a fixed stock level \(b\), the shortage \((X-b)^{+}\) is a random variable whose mean is \(G^{1}(b)\), and its second moment follows directly from Equation C.10:
\[ E\left[\left((X-b)^{+}\right)^{2}\right] = \begin{cases} 2G^{2}(b) + G^{1}(b) & X \text{ discrete} \\[4pt] 2G^{2}(b) & X \text{ continuous} \end{cases} \tag{C.21}\]
so \(\mathit{Var}\left[(X-b)^{+}\right]\) needs nothing beyond \(G^{2}\). Notice the extra \(G^{1}\) in the discrete case. It is there because \((X-b)^{+}(X-b-1)^{+} = \left((X-b)^{+}\right)^{2} - (X-b)^{+}\) when \(X\) and \(b\) are whole numbers, which is the same convention difference Equation C.10 records.
The stock level is not always fixed, and that is what forces a third order function. In the policies of Chapter 8 the quantity the demand is compared against is itself random, and the shortage has to be averaged over its distribution. Averaging \(G^{1}\) over a range produces a difference of \(G^{2}\), by the tail relations of Equation C.12 and Equation C.16. Averaging Equation C.21 over the same range produces a difference of the function one order higher again, and that function is not yet defined.
Thus, define the loss function of order \(n\) by extending Equation C.10 in the obvious way:
\[ G^{n}(b) = \begin{cases} \dfrac{1}{n!}E\left[(X-b)^{+}(X-b-1)^{+}\cdots(X-b-n+1)^{+}\right] & X \text{ discrete} \\[8pt] \dfrac{1}{n!}E\left[\left((X-b)^{+}\right)^{n}\right] & X \text{ continuous} \end{cases} \tag{C.22}\]
and in both cases each order is the tail accumulation of the one below it:
\[ G^{n}(b) = \sum_{x > b}G^{n-1}(x) \qquad\text{or}\qquad G^{n}(b) = \int_{b}^{\infty}G^{n-1}(x)\,dx \tag{C.23}\]
The book needs \(n = 3\) and no higher. Setting \(n = 3\) in Equation C.22 gives the third order loss function, and Equation C.23 gives the identity that makes it useful:
\[ \sum_{j = b+1}^{c} G^{2}(j) = G^{3}(b) - G^{3}(c) \tag{C.24}\]
That is, a sum of second order loss functions across a range of stock levels collapses to a difference of third order loss functions at its endpoints, in the way a sum of first order functions collapses to a difference of second order ones. Chapter 8 reports the mean backorder level with the second identity and its variance with Equation C.24.
C.3.7 Closed Forms for the Third Order
Table C.7 gives \(G^{3}\) for the families whose lower orders appear in Table C.4 and Table C.6. The standard normal values are in the last column of Table C.2.
| Distribution | \(G^{3}(x)\) |
|---|---|
| Geometric | \(\beta^{3}(1-p)^{x}\) |
| Poisson, negative binomial | Computed by the recursion of Equation C.26 |
| Exponential | \(\dfrac{1}{\lambda^{3}}e^{-\lambda x}\) |
| Gamma, shape \(\alpha\) and scale \(\beta\) | \(\tfrac{1}{6}\left[\alpha^{(3)}\beta^{3}G^{0}_{\alpha+3}(x) - 3x\,\alpha^{(2)}\beta^{2}G^{0}_{\alpha+2}(x) + 3x^{2}\alpha\beta\,G^{0}_{\alpha+1}(x) - x^{3}G^{0}_{\alpha}(x)\right]\) |
| Normal | \(\sigma^{3}\Phi^{3}(z)\), with \(\Phi^{3}(z) = \tfrac{1}{6}\left[(z^{2}+2)\phi(z) - z(z^{2}+3)\Phi^{0}(z)\right]\) |
| Lognormal | \(\tfrac{1}{6}\left[H^{3}(x) - 3xH^{2}(x) + 3x^{2}H^{1}(x) - x^{3}G^{0}(x)\right]\), with \(H^{3}(x) = E[X^{3}]\,\Phi(3\sigma - z)\) |
Two of those rows generalize a pattern already visible. The geometric has \(G^{1} = \beta(1-p)^{x}\) and \(G^{2} = \beta^{2}(1-p)^{x}\), and the exponential has \(1/\lambda\) and \(1/\lambda^{2}\) against the same exponential factor, so in both families each order contributes one more factor of the mean and nothing else. Thus, for those two distributions \(G^{n}\) is known at every order without any further work.
The gamma and lognormal rows are Equation C.18 carried one order up. Expanding \((u-x)^{3}\) inside the integral of Equation C.22 and collecting terms gives
\[ G^{3}(x) = \tfrac{1}{6}\left[H^{3}(x) - 3xH^{2}(x) + 3x^{2}H^{1}(x) - x^{3}G^{0}(x)\right], \qquad H^{3}(x) = \int_{x}^{\infty}u^{3}g(u)\,du \tag{C.25}\]
so a third order loss function for any continuous family follows from its third partial expectation. For the gamma, \(u\,g_{\alpha,\beta}(u) = \alpha\beta\,g_{\alpha+1,\beta}(u)\) applied three times gives \(H^{3}(x) = \alpha^{(3)}\beta^{3}G^{0}_{\alpha+3}(x)\), which is the gamma row.
For the discrete families the practice of Equation C.13 extends unchanged. Notice that \(G^{3}(0) = \tfrac{1}{6}E[X(X-1)(X-2)]\), the third binomial moment, so for \(b \ge 0\)
\[ G^{3}(b) = \tfrac{1}{6}E[X(X-1)(X-2)] - \sum_{0 < x \le b}G^{2}(x) \tag{C.26}\]
and the same accumulating loop that produces \(G^{0}\), \(G^{1}\) and \(G^{2}\) produces \(G^{3}\) at no extra cost in evaluations of the mass function. The third binomial moment is \(\lambda^{3}/6\) for the Poisson and \(r(r+1)(r+2)\beta^{3}/6\) for the negative binomial.
For example, take demand over a lead time as \(N(50, 12^{2})\) units and a stock level of 62 units, the case worked in Section C.3.4. Then \(z = 1.00\) and from Table C.2, \(\Phi^{3}(1.00) = 0.0152\), so
\[ G^{3}(62) = \sigma^{3}\Phi^{3}(z) = 12^{3}(0.0152) = 1728(0.0152) = 26.27 \text{ units}^{3} \]
Carried to full precision that is 26.29 units cubed, and the difference is again the four decimal places Table C.2 is printed to.
Notice the units. \(G^{1}\) is in units, \(G^{2}\) in units squared and \(G^{3}\) in units cubed, and the scaling by \(\sigma\), \(\sigma^{2}\) and \(\sigma^{3}\) is what produces them. That check catches a missing power of \(\sigma\) in the same way Section C.3.4 catches a missing \(\sigma\).
C.3.8 Mixtures
Chapter 7 describes demand that is intermittent, and a mixture is one way to represent it. Let \(X\) be distributed as a mixture of \(X_1\) and \(X_2\) with weights \(p_{1}\) and \(p_{2} = 1-p_{1}\), so that
\[ g_{X}(x) = p_{1}g_{X_{1}}(x) + p_{2}g_{X_{2}}(x), \qquad G_{X}(x) = p_{1}G_{X_{1}}(x) + p_{2}G_{X_{2}}(x) \tag{C.27}\]
Because the loss functions are expectations and the expectation operator is linear, they mix the same way:
\[ G_{X}^{1}(x) = p_{1}G_{X_{1}}^{1}(x) + p_{2}G_{X_{2}}^{1}(x), \qquad G_{X}^{2}(x) = p_{1}G_{X_{1}}^{2}(x) + p_{2}G_{X_{2}}^{2}(x) \tag{C.28}\]
and so do the raw moments, \(E[X^{n}] = p_{1}E[X_{1}^{n}] + p_{2}E[X_{2}^{n}]\). The variance does not, and the extra term is the reason mixtures are useful for intermittent demand:
\[ \mathit{Var}[X] = p_{1}\mathit{Var}[X_{1}] + p_{2}\mathit{Var}[X_{2}] + p_{1}p_{2}\left(E[X_{1}] - E[X_{2}]\right)^{2} \tag{C.29}\]
Notice that last term. Two components with identical variances still produce a mixture more variable than either, and the excess grows with the square of the distance between their means. Thus, a mixture can reach a variability the component families cannot, which is what an intermittent demand series demands of a model.
What a mixture does not inherit is an inverse. There is no closed form for \(G_{X}^{-1}(\cdot)\), which the newsvendor solution of Chapter 7 needs. However, the quantile is bracketed by the component quantiles. For any \(p \in [0,1]\), with \(x_{1} = G_{X_{1}}^{-1}(p)\) and \(x_{2} = G_{X_{2}}^{-1}(p)\),
\[ \min\{x_{1}, x_{2}\} \le G_{X}^{-1}(p) \le \max\{x_{1}, x_{2}\} \tag{C.30}\]
Thus, the search has a starting interval that is free. For a discrete mixture, step up from \(\min\{x_{1}, x_{2}\}\); for a continuous one, use the bracket as the initial interval of a bisection search. Either way the bracket is what turns an unbounded search into a bounded one.
C.3.9 Computing Them
The KSL implements these functions directly, so none of the closed forms in Section C.3.5 has to be coded by hand. A distribution that provides them implements LossFunctionDistributionIfc, which adds two methods:
fun firstOrderLossFunction(x: Double): Double
fun secondOrderLossFunction(x: Double): DoubleNine distributions implement it, and they are the ones tabulated above: Poisson, Geometric, NegativeBinomial, Binomial and DEmpiricalCDF for the discrete case, and Exponential, Gamma, Normal and Lognormal for the continuous case. ShiftedLossFunctionDistribution wraps any of them when the support has to be moved off the origin.
val ltd = Normal(mean = 50.0, variance = 144.0)
println(ltd.complementaryCDF(62.0)) // 0.158655, the shortage probability
println(ltd.firstOrderLossFunction(62.0)) // 0.999786, the expected shortage in units
println(ltd.secondOrderLossFunction(62.0)) // 5.424464The Erlang is the Gamma with an integer shape parameter, and needs no separate class.
Check any implementation before trusting it, including this one. Two identities make that cheap, and both come from earlier in this section. For a non-negative random variable,
\[ G^{1}(0) = E[X], \qquad G^{2}(b) \ge 0 ext{ for every } b \tag{C.31}\]
The first is Equation C.13 and its continuous counterpart read at \(b = 0\). The second holds because \(G^{2}\) is the expectation of a product of two non-negative quantities. Thus, an implementation whose \(G^{1}(0)\) differs from the distribution’s own mean has its parameters crossed somewhere, and one that returns a negative \(G^{2}\) is wrong outright. That check is worth performing every time a loss function is used from a new source, because it is one line and it catches the two errors that actually happen: a rate used where a mean belongs, and a truncation applied at the wrong end.
C.4 Choosing a Demand Distribution
Every model from Chapter 7 onward takes a demand distribution as an input, and the distribution carries more of the answer than the policy formula built on top of it does. This section covers the two routes to one. Fitting uses a sample of demand. Moment matching uses only a mean and a variance, which is the situation whenever the quantity of interest was derived rather than observed.
C.4.1 Fitting to a Sample
The sequence is the same one Simulation Modeling using the KSL gives at length in its distribution modeling appendix, compressed here to what an inventory analyst needs.
- Characterize before hypothesizing. Plot the series in time order, then as a histogram. Compute the mean, the variance, and the diagnostic of Equation C.1 or Equation C.5. Count the zeros.
- Hypothesize from the characterization, not from habit. A variance to mean ratio near one points at the Poisson, above one at the negative binomial, below one at the binomial. A coefficient of variation past 0.5 rules out the normal by Table C.1.
- Estimate the parameters, by maximum likelihood or by matching moments.
- Test the fit, with a chi-squared test for a discrete family or a Kolmogorov-Smirnov test for a continuous one.
- Look at the fit, as a density over a histogram and as a quantile plot. The tail is what an inventory decision uses, and a test statistic averages over the whole distribution.
The KSL does the last four steps. PDFModeler for continuous data and PMFModeler for discrete data estimate parameters across a set of families, score the candidates, and rank them:
val results = PDFModeler(data).estimateAndEvaluateScores()Three cautions are worth more than the mechanics.
Failing to reject is not proof. A goodness-of-fit test asks whether the data are inconsistent with a family. With two hundred observations it will fail to reject several families at once, and with twenty thousand it will reject all of them. Thus, the test narrows the field and does not choose.
A sales history is not a demand history. Sales are censored at whatever was on the shelf, so fitting to sales understates both the mean and the variance, and understates them most in the periods when the item mattered. Chapter 7 returns to this.
The reporting bucket changes the answer. Recall from Example 2.4 that the same transactions bucketed daily rather than monthly show more zeros and a longer interval between demands, which can move an item from one demand class to another. Thus, fit at the bucket the decision is made in, and state which one it was.
C.4.2 Matching the First Two Moments
Often there is no sample. The demand during a lead time, which Chapter 8 needs, is a random sum whose mean and variance follow from the moments of demand and of lead time by Section B.5, and no observations of it exist at all. What is available is \(\hat{\mu}\) and \(\hat{\sigma}^{2}\), and what is wanted is a distribution.
For example, take demand averaging 5 units a week with a variance of 16, and a lead time averaging 10 weeks with a variance of 25. Then
\[ E[D(L)] = 10(5) = 50 \text{ units}, \qquad \mathit{Var}[D(L)] = 10(16) + 25(5)^{2} = 785 \tag{C.32}\]
so \(c^{2} = 785/2500 = 0.3140\) and \(\mathit{VMR} = 15.70\) units. Notice that the lead time variance enters multiplied by the square of the demand rate, which is why an unreliable supplier dominates a variable demand in almost every real instance.
The gamma match is the default. The gamma can match any mean and variance, and inverting Equation C.7 gives the parameters directly:
\[ \hat{\beta} = \frac{\hat{\sigma}^{2}}{\hat{\mu}}, \qquad \hat{\alpha} = \frac{\hat{\mu}^{2}}{\hat{\sigma}^{2}} \tag{C.33}\]
On Equation C.32 that gives \(\hat{\beta} = 785/50 = 15.70\) and \(\hat{\alpha} = 2500/785 = 3.1847\). The check is immediate: \(\alpha\beta = 50.0000\) and \(\alpha\beta^{2} = 785.0000\). That check is worth performing every time, because a matching procedure that does not reproduce its own targets is the most common failure in this material.
A mixture of two Erlangs is the alternative, and it is the one to use when the distribution must be discrete in spirit or when a phase-type representation is wanted. Algorithm C.1 returns \(\mathrm{Erlang}(k_{1}, \lambda_{1}, k_{2}, \lambda_{2}, p_{1}, p_{2})\) matching \(E[X]\) and \(\mathit{Var}[X]\).
disc, which is exactly zero when \(c^{2}\) is the reciprocal of an integer and can go slightly negative in floating point, so it is clamped rather than trusted. At that same point \(p_{1}\) evaluates to one, and to a shade above one in floating point, so it is clamped as well. Neither clamp changes an answer; both stop the procedure returning something that is not a probability.
erlangMixture(E[X], Var[X]):
c2 = Var[X] / (E[X] * E[X])
if c2 < 1: // both components are Erlang
k1 = floor(1 / c2)
k2 = k1 + 1
disc = max(0, k2*(1 + c2) - k2*k2*c2) // exactly 0 when c2 = 1/k
p1 = min(1, (k2*c2 - sqrt(disc)) / (1 + c2))
lambda1 = (k2 - p1) / E[X]
lambda2 = lambda1
else: // both components are exponential
k1 = 1
k2 = 1
lambda1 = (2 / E[X]) * (1 + sqrt((c2 - 1/2) / (c2 + 1)))
lambda2 = 4 / E[X] - lambda1
p1 = lambda1 * (lambda2 * E[X] - 1) / (lambda2 - lambda1)
p2 = 1 - p1
return (k1, lambda1, k2, lambda2, p1, p2)On Equation C.32, \(c^{2} = 0.3140\) falls in the first branch, giving \(k_{1} = 3\), \(k_{2} = 4\), \(p_{1} = 0.5893\) and a common rate of \(\lambda = 0.06821\) per unit. The mixture’s mean is 50.0000 and its variance 785.0000, matching to six figures.
One implementation note. When \(c^{2}\) is exactly \(1/k\) for an integer \(k\), the quantity under the square root in the first branch is exactly zero, and floating point arithmetic can make it very slightly negative. Thus, clamp it at zero rather than letting the square root fail.
Adan’s rule selects a family rather than assuming one, and it is the procedure to use when the quantity being represented is known to be discrete. Let
\[ a = c^{2} - \frac{1}{\mu} = \frac{\mathit{VMR} - 1}{\mu} \tag{C.34}\]
Notice that \(a\) is a rescaled statement about how the variance compares with the mean, so \(a = 0\) is the Poisson case. Table C.8 gives the rule, in the form of Adan et al. (1995) as extended by Rossetti and Ünlü (2011) to cover \(a < -1\), where no discrete family can match.
| Condition | Family |
|---|---|
| \(a < -1\) | Gamma, by Equation C.33 |
| \(-1 \le a < 0\) | A mixture of \(\mathrm{Bin}(n, p)\) and \(\mathrm{Bin}(n+1, p)\) |
| \(a = 0\) | Poisson with mean \(\mu\) |
| \(0 < a < 1\) | A mixture of \(\mathrm{NB}(r, p)\) and \(\mathrm{NB}(r+1, p)\) |
| \(a \ge 1\) | A mixture of \(\mathrm{Geom}(p_{1})\) and \(\mathrm{Geom}(p_{2})\) |
In every mixture below, \(q\) is the weight on the first component.
Mixed binomial, for \(-1 \le a < 0\):
\[ n = \left\lfloor \frac{-1}{a} \right\rfloor, \qquad q = \frac{1 + a(n+1) + \sqrt{-an(n+1) - n}}{1+a}, \qquad p = \frac{\mu}{n+1-q} \tag{C.35}\]
Mixed negative binomial, for \(0 < a < 1\):
\[ r = \left\lfloor \frac{1}{a} \right\rfloor, \qquad q = \frac{a(r+1) - \sqrt{(r+1)(1-ar)}}{1+a}, \qquad p = \frac{\mu}{r+1-q+\mu} \tag{C.36}\]
Mixed geometric, for \(a \ge 1\):
\[ q = \frac{1}{1 + a + \sqrt{a^{2}-1}}, \qquad p_{1} = \frac{\mu}{2q + \mu}, \qquad p_{2} = \frac{\mu}{2(1-q) + \mu} \tag{C.37}\]
In Equation C.36 and Equation C.37, \(p\) is the probability of a failure, so that \(\mathrm{NB}(r,p)\) has mean \(rp/(1-p)\) and \(\mathrm{Geom}(p)\) has mean \(p/(1-p)\). That is the opposite of the convention in Equation C.3, and it is the single most likely thing to go wrong when these are implemented. The check of Equation C.33 applies here too and catches it at once: compute the mixture’s mean and variance from its components and compare them with the targets.
On Equation C.32, \(a = 0.3140 - 1/50 = 0.2940\), which falls in the mixed negative binomial branch. Then \(r = \lfloor 1/0.2940 \rfloor = 3\), \(q = 0.37788\) and \(p = 0.93245\), and the resulting mixture has mean 50.0000 and variance 785.0000.
Notice that all three routes reproduce the same two moments and disagree about everything else. That is the lesson of moment matching: two numbers do not determine a distribution, and the choice of family is an assumption that the procedure hides rather than removes. Where the answer matters, compute it under two of them and report the spread.
C.4.3 When the Demand Is Intermittent
A series in which most periods have no demand at all is not served well by any of the above, because a single family has to explain both how often demand occurs and how large it is when it does. Two responses are available.
The first is a zero-modified family, which mixes a point mass at zero with one of the families of Section C.1. Rossetti and Ünlü (2011) recommend a zero-modified form of Adan’s rule when the probability of a zero period exceeds one half, and switching back to the ordinary rule below that.
The second is a compound representation, in which the time between demands and the size of a demand are modeled separately and the period demand is the sum of a random number of sizes. Section B.5 gives the moments. That representation is the more faithful one, and it is the one that connects to the demand classification of Section 2.2.6, where the interval between demands and the variability of the sizes are the two axes.
Both are beyond what Chapter 7 needs, and both matter from Chapter 8 onward, where the quantity being modeled is demand over a lead time and the zeros no longer conveniently aggregate away.