Fitting the Demand, and Solving the Buy

Manuel D. Rossetti, PhD P.E.

Agenda

  • Look, summarize, hypothesize, estimate, test
  • The one statistic that decides the family before any test runs
  • What a goodness of fit test does, and does not, settle
  • When a history is mostly zeros
  • The belt, solved: 475 against a mean of 465
  • Why the optimum is flat, and what the choice of family costs

The Distribution Carries More Than the Formula

Every model from here to the end of the course 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.

The distributor has 104 weeks of demand for the belt, recorded after the deck went out of production.

One difference from the general procedure matters here: demand for a spare part is a count, so every tool is a discrete one.

Read the Data as Data

val df = DataFrame.readCSV(
    File("data/ch7-weekly-demand.csv"),
    colTypes = mapOf("week" to ColType.Int, "demand" to ColType.Int)
)
val history: IntArray = df.getColumn("demand").toList().map { it as Int }.toIntArray()

Notice that the column types are named rather than left to be inferred.

A column of counts silently read as text, or as real numbers, fails later and somewhere else.

Look Before You Model

A summary statistic answers a question you already thought to ask.

A picture shows you the thing you did not.

Three plots answer three different questions.

val h = Histogram(Histogram.createBreakPoints(0.0, 12, 10.0), "Weekly demand")
h.collect(history.toDoubles())
h.histogramPlot().showInBrowser()

ObservationsPlot(history.toDoubles()).showInBrowser()
ACFPlot(history.toDoubles()).showInBrowser()

What Shape Is It?

Weekly demand over 104 weeks, in bins of ten belts.

Reading the Histogram

A peak well to the left.

A long thin tail reaching past 110.

No suggestion of symmetry anywhere in it.

That single picture has already ruled out the normal distribution.

The Weeks in the Order They Occurred

Is It the Same Process Throughout?

Plot the weeks in the order they occurred. Look for a trend, a level shift, or a season.

There is none of those: the series wanders about a constant level for two years.

That matters, because it means the season really is over, which the problem statement claimed and this plot lets a reader check rather than accept.

Autocorrelation Through Twenty Lags

Are the Weeks Independent?

The autocorrelation plot. Every lag falls inside the band, including the first, whose correlation is 0.1423 against a band of about

\[ \pm 1.96/\sqrt{104} = \pm 0.19 \]

Autocorrelated demand breaks the variance of a sum.

We are about to add fourteen of these weeks together and treat their variances as adding. This is the plot that licenses that arithmetic.

Taken together: the weeks are independent, identically distributed, and strongly right-skewed.

Summarize It

Statistic Value Statistic Value
Weeks observed 104 Skewness 1.2630
Mean 33.2404 Smallest week 0
Variance 613.0387 First quartile 16.0
Standard deviation 24.7596 Median 26.5
Coefficient of variation 0.7449 Third quartile 40.0
Variance to mean ratio 18.4426 Largest week 113
Weeks with no demand 2

Read three things off this before going further.

One: The Median Sits Well Below the Mean

A median of 26.5 against a mean of 33.24 says the distribution is right-skewed: most weeks are quiet and a few are very large.

The skewness of 1.2630 is the same statement with a number attached.

The histogram is the same statement with a picture attached.

Three routes, one finding.

Two: Demand Is a Count

The weeks are whole belts. Two of them are zero. The quantity has no meaning between whole numbers.

So the families to consider are the discrete ones, and not the continuous ones.

Three: For a Count, the Diagnostic Is the Variance to Mean Ratio

Not the coefficient of variation, because the ratio is what each discrete family fixes or bounds.

\[ \mathit{VMR} = \frac{613.0387}{33.2404} = 18.44 \]

Step 1, Hypothesize

The ratio does most of this step by itself.

The Poisson fixes \(\mathit{VMR} = 1\). The binomial forces \(\mathit{VMR} < 1\).

A ratio of 18.44 is near neither, and no choice of parameter repairs that, because in both families the variance is decided the moment the mean is estimated.

What remains is the negative binomial, which admits any \(\mathit{VMR} > 1\), and which is the Poisson with its constant-rate assumption relaxed.

That is a plausible story for a spare part: the rate at which mower decks need belts is not the same in every week.

Step 2, Estimate

val modeler = PMFModeler(history)
for (result in modeler.estimateParameters(modeler.defaultEstimators)) {
    println(if (result.success) "${result.parameters}" else "FAILED: ${result.message}")
}
RV Type = NegativeBinomial
Double Parameters {probOfSuccess=0.05422232292974839, numSuccesses=1.9057024844439698}
RV Type = Poisson
Double Parameters {mean=33.24038461538462}
FAILED: Cannot match moments when sample average <= estimated variance
RV Type = Binomial
Double Parameters {probOfSuccess=0.29133420825260514, numTrials=114.09708737864078}

The most useful line in that output is the failure.

The Software Performed Step 1

The moment-matching binomial estimator refuses to return an answer, and the reason it gives is the variance to mean ratio in words.

The binomial that did return an answer did so by matching the largest observation instead of the moments.

It reports 114.1 trials: a fractional number of trials, for a quantity that has no trials in it.

An estimate being returned is not the same thing as a family being appropriate.

Confirm the Estimate by Hand

Matching the negative binomial to the sample mean and variance:

\[ \hat{p} = \frac{\bar{x}}{s^{2}} = \frac{33.2404}{613.0387} = 0.054222 \]

\[ \hat{r} = \frac{\bar{x}^{2}}{s^{2}-\bar{x}} = \frac{33.2404^{2}}{613.0387-33.2404} = 1.9057 \]

Notice that \(1/\hat{p} = 18.44\) is the variance to mean ratio we started with.

So the fitted \(\hat{p}\) carries no information the summary did not already have.

That is true of the method of moments generally: the fit is a restatement of two summary numbers in the parameters of a family.

Step 3, Test

A chi-squared goodness of fit test, on bins of approximately equal probability.

Model \(\chi^{2}\) d.f. \(p\)-value At the 0.05 level
Negative binomial, \(\hat{p}=0.0542\), \(\hat{r}=1.9057\) 23.35 17 0.1382 do not reject
Poisson, \(\hat{\lambda}=33.2404\) 571.98 16 0.0000 reject

The Poisson is not rejected narrowly. Its test statistic is more than twenty times the negative binomial’s.

The table is the confirmation, not the decision.

The decision was made by the variance to mean ratio. A test that contradicted it would be telling us something is wrong with the data.

The Fitted Model Against the Data

Reading the Comparison

The empirical mass function is spiky because it has to be: 104 observations spread over 54 distinct values put a mass of \(1/104\) or \(2/104\) on most of them.

Which is exactly why the test bins before it tests.

What to read here is whether the fitted curve runs through the middle of the empirical forest.

It does.

What the Test Does Not Settle

A goodness of fit test asks whether the data are inconsistent with a family. Failing to reject is not proof.

With 104 observations over twenty bins, several bins carry an expected count below five, which is where the chi-squared approximation is least trustworthy.

The test and the ratio narrow the field. The choice among what survives is yours to defend.

Where two survivors disagree enough to change the decision: compute the answer under both and report the spread.

And What It Did Not Do

It did not deseasonalize, because there is nothing to remove and the season is over.

It did not forecast, which is a different job from fitting a distribution.

It did not correct for censoring, because the distributor does not record lost sales.

That last one is the limitation from last session, and no amount of care here repairs it.

When the History Is Mostly Zeros

The belt sold in 102 of its 104 weeks. Many items are not like that.

Classify a history by two statistics: the average interval between demands \(\mathit{ADI}\), and the squared coefficient of variation of the non-zero sizes \(\mathit{CV}^{2}\).

\(\mathit{ADI} > 1.32\) is intermittent. Also \(\mathit{CV}^{2} > 0.49\) is lumpy.

Those cut-offs are not conventions. They are the points at which one forecasting method overtakes another, which is why they are so oddly specific.

For the belt, \(\mathit{ADI} = 104/102 = 1.02\). Comfortably under, which is the check that licensed everything above.

Two Things Go Wrong Otherwise

A single family has to explain two different things. Period demand for an intermittent item is a mixture: a point mass at zero, and a distribution of sizes. Fitting one family forces it to account for how often demand happens and how large it is at once.

The coefficient of variation goes past where the normal is usable. A normal model at \(CV = 1.0\) puts a sixth of its probability below zero, and intermittent items routinely run past 1.0.

That is not a small error in a tail. It is a model asserting that demand is frequently negative.

Forecasting Intermittent Demand Is a Different Question

Worth naming precisely enough that you can go and read it.

Exponentially smoothing an intermittent series directly is wrong: the smoothing is dragged toward zero by every empty period, so the estimate is lowest just after a demand occurs, which is precisely when the next one is least likely.

Croston’s method smooths two series instead, the size when a demand happens and the interval between demands, and combines them into a rate.

The obvious combination is biased: the ratio of two independent smoothed estimates does not estimate the ratio of their means. Syntetos and Boylan gave the correction.

Fitting asks what the distribution of demand is. Forecasting asks what next period’s demand will be. This chapter asks only the first.

The Belt, Solved

Let \(X\) represent total demand over the fourteen remaining weeks.

Weekly demands are independent, so the mean and variance of a sum add:

\[ E[X] = 14(33.2404) = 465.37, \qquad \mathit{Var}[X] = 14(613.0387) = 8582.54 \]

So \(\sigma = 92.64\) belts, and the coefficient of variation is \(92.64/465.37 = 0.1991\).

Two Approximations, Both Earned

Treating the total as continuous. One belt is \(1/92.64\), or 1.1%, of the standard deviation of the total. Rounding cannot decide the answer.

The weekly data, where one belt is 4% of the standard deviation, would not have earned the same treatment.

Using the normal. We rejected the normal for weekly demand, at \(CV = 0.7449\). The fourteen-week total has \(CV = 0.1991\), and a normal at that coefficient puts essentially nothing below zero.

Safe here and not safe there, and the reason is that a sum of many independent terms tends toward the normal whatever the terms look like.

The Quantity

\[ \frac{c_{u}}{c_{u}+c_{o}} = \frac{45}{83} = 0.542169 \]

Standardizing, \(Q^{*} = \mu + z\sigma\) with \(z = \Phi^{-1}(0.542169) = 0.1059\):

\[ Q^{*} = 465.37 + 0.1059(92.64) = 475.18 \]

The distributor orders 475 belts.

The Consequences

The standard normal loss function at \(z = 0.1059\) is 0.3482, so the expected shortage is \(92.64(0.3482) = 32.26\) belts.

The expected leftover is \(32.26 + 475.18 - 465.37 = 42.07\) belts.

The distributor expects to scrap about 42 belts and to turn away about 32 dealers.

\[ E[C(Q^{*})] = 38(42.07) + 45(32.26) = 1598.71 + 1451.72 = \$3{,}050.43 \]

The Check

The other form of the expected cost must agree:

\[ (45+38)(32.26) + 38(475.18 - 465.37) = 2677.62 + 372.81 = \$3{,}050.43 \]

The two agree, which confirms both the loss function and the quantity computed from it.

Worth performing every time. It catches a missing \(\sigma\) in the loss function immediately.

In business terms: expected profit is \(45(465.37) - 3050.43 = \$17{,}891\). The distributor commits $23,750 to buy 475 belts and expects to clear about $17,900.

The Optimum Is Flat

Ordering the mean rounded to whole belts, 465, costs $3,068.88 against the optimum’s $3,050.43.

That is 0.6% worse.

The EOQ is flat in the same way and for the same reason: the objective is smooth and its derivative vanishes at the optimum.

So the value of this calculation is not the last belt.

It is knowing that 475 is right and 700 is not, and 700 costs $8,930, nearly three times the optimum.

What the Choice of Family Costs

A sum of independent negative binomials sharing \(p\) is negative binomial, so the fitted weekly model gives the fourteen-week distribution exactly.

Model used to choose \(Q\) Normal Lognormal Gamma Exact negative binomial
\(Q\), rounded to whole belts 475 466 469 469
Expected cost under the exact model $3,061.28 $3,057.19 $3,055.38 $3,055.38
Penalty against the exact answer 0.19% 0.06% 0.00% 0.00%

Read that as a bound, not a ranking. The worst of the three continuous families costs 0.19% more than choosing correctly.

That spread is small because \(CV = 0.1991\). It would not be small at the weekly 0.7449, which is another reason to fit carefully.

The Normal’s Defect Was Real and Did Not Bite

A normal with this mean and standard deviation puts \(\Phi(-5.02)\) below zero.

About one in four million.

That is a statement about this aggregation, not about the normal in general.

The same family applied to the weekly data would have put about 9% of its mass below zero.

What We Did Today

The variance to mean ratio chose the family, and the goodness of fit test confirmed it rather than deciding it.

A fit by moments restates two summary numbers in the parameters of a family, and carries no more information than they do.

The normal was wrong weekly and safe over fourteen weeks, and the reason is the coefficient of variation of the sum.

The optimum is flat. The value of the calculation is not the last belt. It is knowing that 475 is right and 700 is not.

The choice of continuous family cost at most 0.19%, and it would not have been that cheap at the weekly coefficient of variation.

For Next Time

Read before next session: Discrete Demand, Extensions, and Estimating the Answer by Simulation.

Everything today leaned on a continuous approximation that this item earned. Next session is about the items that do not:

  • The staircase rule, which is not an inverse
  • A case where rounding the continuous answer buys one unit too many, at $3,800
  • A stockout penalty, and stock already on hand, and what an unbounded answer means
  • Computing the same answer by sampling, and why you do that on a problem you have already solved

One belt was 1.1% of the standard deviation of the total. One transformer is going to be 40% of it, and that changes the arithmetic.

⌂ Index