Closed formula for D4 constant calculation? (Moving range chart constant)

Question

I need to build a Moving Range Shewhart control chart given a series of observations. In short, I have to calculate the central line and the upper and lower limits as follows

cl = MR
lcl = MR - D3 * MR
ucl = MR + D4 * MR

Where MR is the average moving range and D3, D4 are unbiasing constants. Both D3 and D4 depend on the size of the moving range (n).

My textbook, as most of the books on this subject do, suggests to calculate the constants just by looking at the tables however this is really unpractical when using a software such as R or Python and you need to test multiple moving range sizes, therefore I would like to define a function to calculate D3 and D4 given the moving range size. In order to do this, I need to know the closed formula for D3 and D4 which is not explained anywhere, apparently, or I'm missing something (most likely).

Is there a formula for calculating the D3 and D4 constant for unbiasing the Moving Range chart?

Some more explanation on why these constant are needed would be helpful, the book doesn't say much on the subject.

Answer 1

Repeat the question:
So you want a form for the D4 coefficient for a Shewhart moving range control chart?

Answer:
Here is a table (link), and you can make a really fast lookup for windows within the range. It is piece-wise constant interpolation, unless you have non-integer sample sizes.

Here is a picture in case the link breaks.

I get an okay fit for $ 3 \le n \le 25$ using a linear model.

x <- seq(from=2,to=25,by=1)
y <- c(3.267,       2.574,       2.282,       2.114,
       2.004,       1.924,       1.864,       1.816,
       1.777,       1.744,       1.717,       1.693,
       1.672,       1.653,       1.637,       1.622,
       1.608,       1.597,       1.585,       1.575,
       1.566,       1.557,       1.548,       1.541)

est <- lm (I(log(y ))~ 1 + I(log(x)) +I(log(x)^2)+I(log(x)^(0.345)) )
summary(est)

Where the result is:

$log(D4) = 3.0326040 + 0.2940527*log(n)-0.0063287*(log(n)^2) - 2.3257758*(log(n)^{0.345}) $

The summary values for the fit were:

Call:
lm(formula = I(log(y)) ~ 1 + I(log(x)) + I(log(x)^2) + I(log(x)^(0.345)))

Residuals:
       Min         1Q     Median         3Q        Max 
-4.024e-04 -1.388e-04  1.084e-05  1.242e-04  3.469e-04 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.0326040  0.0058494  518.45  < 2e-16 ***
I(log(x))          0.2940527  0.0037063   79.34  < 2e-16 ***
I(log(x)^2)       -0.0063287  0.0003954  -16.00 7.24e-13 ***
I(log(x)^(0.345)) -2.3257758  0.0091925 -253.01  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0002038 on 20 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 6.188e+06 on 3 and 20 DF,  p-value: < 2.2e-16

It is an un-adjusted $R^2$ of darn near 6-9's. I had to use vegan to get to that. Every parameter has a p-value smaller than 1e-13, on 24 samples. When extrapolation occurs, the function appears for a while to retain "physics". If you look at the first differences of the raw data, there is a very linear trend at progressively higher values of n. It also does not make physical sense that there should be a negative, or even a very small value for D4.

Here is the code to get the R^2

library(vegan)
RsquareAdj(est)

Here is the result. Yes, that is 5-9's.

$r.squared
[1] 0.9999989

$adj.r.squared
[1] 0.9999988

This was how I came to the 0.345 power. It is the midpoint between which the adjusted $R^2$ is constantly 0.9999988. One unit either side of the edge, the value drops to 0.99999987.

Here is the code to get the plot (and yes I am using "y" instead of "d4"):

n <- seq(from=3, to=50,by=1)
y2 <- exp(3.0326040 + 0.2940527*log(n)-0.0063287*(log(n)^2) - 2.3257758*(log(n)^0.345))

plot(n,y2,type="b",col="Green",xlab="n",ylab="D4")
points(x,y,pch=2,col="Red")
legend(x = 35,y = 2.5,legend = c("fit","data"),
       col=c("Green","Red"),
       pch=c(1,2),
       lty=c(1,-2))

Here is the result

I don't like extrapolating, especially without a decent understanding of what is actually going on, but if I was forced to do something then this might be where I would start.

Answer 2

The definition of $D_3$ and $D_4$, from "ASTM MNL7, Manual on Presentation of Data and Control Chart Analysis" is:

$D_3=1-3\frac{d_3}{d_2}$

The equation in the standard should be, $D_3=\max\left(0,1-3\frac{d_3}{d_2} \right)$

$D_4=1+3\frac{d_3}{d_2}$

The factors $d_2$ and $d_3$ are calculated with:

$d_2=\int_{-\infty }^{\infty}\left [ 1-\left ( 1-\alpha_1 ^n-\alpha_1^n\right ) \right ]dx_1$ where $\alpha_1=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x_1}e^{-\left (x^2/2 \right )}dx$ and $n=$sample size.

$d_3=\sqrt{2\int_{-\infty}^\infty\int_{-\infty}^{x_1}\left [ 1-\alpha_1^n-\left ( 1-\alpha_n \right )^n+\left ( \alpha_1-\alpha_n \right )^n \right ]dx_ndx_1-d_2^2}$ where $\alpha_1=\frac{1}{\sqrt{2\pi}\int_{-\infty}^{x_1}e^{-\left (x^2/2 \right )}dx}$ and $\alpha_n=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x_n}e^{-\left (x^2/2 \right )}dx$

As a result, interoplation of values based upon already established tables or simple using a lookup function in an array of table values would be the most efficient method for calculating the desired control limits. Looking at a few different tables (D.C. Mongomery's Introduction to Statistical Quailty Control, ANSI/ASQC B3-1996, ASTM E2587-12, MNL7, and ISO 7870-2:2013), three decimal places for this factor seem to be sufficient. Juran's Quality Handbook: Fifth Edition only has values out to two decimal places.

Montgomery also refers to Grant, E.L., and R.S. Leavenworth (1980). Statistical Quality Control, 5th ed. in using $D$ values based upon probability limits of a distribution and that those results are typically very agreable with the $D_3$ and $D_4$ values. They may provide another option for a closed formula, however Mongomery points to a table of values for the proability-based limits as well.

Answer 3

I know this is a closed answer, but just to extend a little bit of the discussion, and because it's also one of my favorite subjects I would like to add an article I found.

The article it's called: Tables of Range and Studentized Range by H. Leon Harter https://www.jstor.org/stable/2237810?seq=1.

Page 1124 presents the following equation that was used in order to compute the moments of the range statistics.

$$ E\left(W^{k}\right) = n\left(n-1\right)\int_{-\infty}^\infty \Bigg\{ \int_{0}^\infty W^{k} \Big[ \Phi(X+W) - \Phi(X) \Big]^{n-2} \phi(X+W)dW \Bigg\}\phi(X)dx $$

With

$$ \phi(X) = (2\pi)^{-1/2}e^{-X^{2}/2} \qquad \Phi(X) = \int_{0}^X \phi(X)dX $$

So, this is the closed-form formula to compute constants $d_{2}$ and $d_{3}$ that Harter used in 1960. The article itself, on page 1130, has the most complete table I've ever seen, it has values for the constant from $n=2 \ldots 100$ and with 10 decimals figures. Constant $d_{4}$ is computed as @Tavrock has already mentioned.

I've also posted a question on how to solve Harter's integral Harter's Integral, if someone knows how to tackle this problem any comment would be well received.

Answer 4

The constants $d_2$ and $d_3$ can be found by numerical integration. The equations for all values are provided by @Tavrock.

In R, $d_2$ can be found with the following code, with $n$ the sample size:

integrand <- function(x1) { 
    a1 = pnorm(x1)
    1 - a1^n - (1 - a1)^n
}
d2 = integrate(integrand, lower = -Inf, upper = Inf)$value

For $d_3$, a double integration is necessary. Once $d_2$ is calculated, it can be done using the following code:

InnerFunc = function(x1, xn) { 
    a1 = pnorm(x1); an = pnorm(xn)
    1 - a1^n - (1 - an)^n + (a1 - an)^n 
}
InnerIntegral = function(x1) {
    sapply(x1, function(x1) { 
        integrate(function(xn) InnerFunc(x1, xn), -Inf, x1)$value
    }) 
}
integral = integrate(InnerIntegral, -Inf, Inf)$value
d3 = sqrt( 2 * integral - d2^2)

From $d_2$ and $d_3$, we can readily find $D_1$, $D_2$, $D_3$, and $D_4$.

Here are the values I computed, for $n$ ranging from 26 to 100:

n	$d_2$	1/$d_2$	$d_3$	$D_1$	$D_2$	$D_3$	$D_4$
26	3.9643	0.2523	0.7050	1.8494	6.0793	0.4665	1.5335
27	3.9965	0.2502	0.7017	1.8914	6.1016	0.4733	1.5267
28	4.0274	0.2483	0.6986	1.9318	6.1231	0.4797	1.5203
29	4.0570	0.2465	0.6955	1.9704	6.1437	0.4857	1.5143
30	4.0855	0.2448	0.6927	2.0075	6.1635	0.4914	1.5086
31	4.1129	0.2431	0.6899	2.0432	6.1826	0.4968	1.5032
32	4.1393	0.2416	0.6872	2.0776	6.2011	0.5019	1.4981
33	4.1648	0.2401	0.6847	2.1107	6.2189	0.5068	1.4932
34	4.1894	0.2387	0.6822	2.1427	6.2361	0.5115	1.4885
35	4.2132	0.2373	0.6799	2.1736	6.2528	0.5159	1.4841
36	4.2362	0.2361	0.6776	2.2035	6.2690	0.5202	1.4798
37	4.2586	0.2348	0.6754	2.2324	6.2847	0.5242	1.4758
38	4.2802	0.2336	0.6732	2.2604	6.2999	0.5281	1.4719
39	4.3012	0.2325	0.6712	2.2876	6.3147	0.5319	1.4681
40	4.3216	0.2314	0.6692	2.3140	6.3291	0.5355	1.4645
41	4.3414	0.2303	0.6673	2.3396	6.3431	0.5389	1.4611
42	4.3606	0.2293	0.6654	2.3645	6.3568	0.5422	1.4578
43	4.3794	0.2283	0.6636	2.3887	6.3700	0.5454	1.4546
44	4.3976	0.2274	0.6618	2.4123	6.3830	0.5485	1.4515
45	4.4154	0.2265	0.6601	2.4352	6.3956	0.5515	1.4485
46	4.4328	0.2256	0.6584	2.4576	6.4080	0.5544	1.4456
47	4.4497	0.2247	0.6568	2.4794	6.4200	0.5572	1.4428
48	4.4662	0.2239	0.6552	2.5007	6.4318	0.5599	1.4401
49	4.4824	0.2231	0.6536	2.5214	6.4433	0.5625	1.4375
50	4.4981	0.2223	0.6521	2.5417	6.4546	0.5651	1.4349
51	4.5136	0.2216	0.6507	2.5615	6.4656	0.5675	1.4325
52	4.5286	0.2208	0.6492	2.5809	6.4764	0.5699	1.4301
53	4.5434	0.2201	0.6479	2.5998	6.4870	0.5722	1.4278
54	4.5578	0.2194	0.6465	2.6183	6.4973	0.5745	1.4255
55	4.5720	0.2187	0.6452	2.6365	6.5075	0.5767	1.4233
56	4.5858	0.2181	0.6439	2.6542	6.5174	0.5788	1.4212
57	4.5994	0.2174	0.6426	2.6716	6.5272	0.5809	1.4191
58	4.6127	0.2168	0.6413	2.6887	6.5367	0.5829	1.4171
59	4.6258	0.2162	0.6401	2.7054	6.5461	0.5848	1.4152
60	4.6386	0.2156	0.6389	2.7217	6.5554	0.5868	1.4132
61	4.6511	0.2150	0.6378	2.7378	6.5645	0.5886	1.4114
62	4.6635	0.2144	0.6366	2.7536	6.5734	0.5905	1.4095
63	4.6756	0.2139	0.6355	2.7690	6.5821	0.5922	1.4078
64	4.6875	0.2133	0.6344	2.7842	6.5907	0.5940	1.4060
65	4.6992	0.2128	0.6333	2.7991	6.5992	0.5957	1.4043
66	4.7106	0.2123	0.6323	2.8138	6.6075	0.5973	1.4027
67	4.7219	0.2118	0.6313	2.8282	6.6157	0.5989	1.4011
68	4.7330	0.2113	0.6302	2.8423	6.6238	0.6005	1.3995
69	4.7440	0.2108	0.6292	2.8562	6.6317	0.6021	1.3979
70	4.7547	0.2103	0.6283	2.8699	6.6395	0.6036	1.3964
71	4.7653	0.2099	0.6273	2.8834	6.6472	0.6051	1.3949
72	4.7757	0.2094	0.6264	2.8966	6.6548	0.6065	1.3935
73	4.7860	0.2089	0.6254	2.9096	6.6623	0.6080	1.3920
74	4.7960	0.2085	0.6245	2.9225	6.6696	0.6093	1.3907
75	4.8060	0.2081	0.6236	2.9351	6.6769	0.6107	1.3893
76	4.8158	0.2077	0.6227	2.9475	6.6840	0.6121	1.3879
77	4.8254	0.2072	0.6219	2.9598	6.6911	0.6134	1.3866
78	4.8349	0.2068	0.6210	2.9718	6.6980	0.6147	1.3853
79	4.8443	0.2064	0.6202	2.9837	6.7049	0.6159	1.3841
80	4.8535	0.2060	0.6194	2.9954	6.7117	0.6172	1.3828
81	4.8627	0.2056	0.6186	3.0070	6.7183	0.6184	1.3816
82	4.8716	0.2053	0.6178	3.0184	6.7249	0.6196	1.3804
83	4.8805	0.2049	0.6170	3.0296	6.7314	0.6208	1.3792
84	4.8893	0.2045	0.6162	3.0407	6.7378	0.6219	1.3781
85	4.8979	0.2042	0.6154	3.0516	6.7442	0.6230	1.3770
86	4.9064	0.2038	0.6147	3.0624	6.7504	0.6242	1.3758
87	4.9148	0.2035	0.6139	3.0730	6.7566	0.6253	1.3747
88	4.9231	0.2031	0.6132	3.0835	6.7627	0.6263	1.3737
89	4.9313	0.2028	0.6125	3.0939	6.7688	0.6274	1.3726
90	4.9394	0.2025	0.6118	3.1041	6.7747	0.6284	1.3716
91	4.9474	0.2021	0.6111	3.1142	6.7806	0.6295	1.3705
92	4.9553	0.2018	0.6104	3.1241	6.7864	0.6305	1.3695
93	4.9631	0.2015	0.6097	3.1340	6.7922	0.6315	1.3685
94	4.9708	0.2012	0.6090	3.1437	6.7979	0.6324	1.3676
95	4.9784	0.2009	0.6084	3.1533	6.8035	0.6334	1.3666
96	4.9859	0.2006	0.6077	3.1628	6.8091	0.6343	1.3657
97	4.9934	0.2003	0.6071	3.1722	6.8146	0.6353	1.3647
98	5.0007	0.2000	0.6064	3.1814	6.8200	0.6362	1.3638
99	5.0080	0.1997	0.6058	3.1906	6.8254	0.6371	1.3629
100	5.0152	0.1994	0.6052	3.1996	6.8307	0.6380	1.3620

For $n$ between 26 and 30, these results agree with Kenith Grey's simulations.

We can compare the numerical results to the approximation provided by @EngrStudent: