Q:

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

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.

A:

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.

A:

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.

