%% The exponentiallyâ€weighted control charts for geometrically moving averages
% Authors: Saeed Maghsoodloo, Daniel F Silva
% Quality and Reliability Engineering International 39.7 (2023), 2706-2719

%% Inputs
% w - omega: the weight of the present observation
% kvec: values of k, number of standard deviations in alternate hypothesis
% Lvec: values of L, parameter for width of the control limit
format longE
w=0.2;
w1=1-w;
n=6;

kvec = [0.00, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2];
Lvec = [2.70, 2.80 , 2.90 , 3.00 , 3.10];

knum = length(kvec);
Lnum = length(Lvec);

ARLmat= zeros(knum,Lnum);

%% Section 6 - The average run-length when sigma calculated from subgroups
% BetaLimit: Type II error probability for large t. Eq (22)
% g, fun, PrE_fun : intermediate calculations to obtain BetaLimit
g=sqrt(w/(1+w1));
fun = @(Z1,Z2) 1/(2.*pi.*sqrt(1-w1.^2)).*exp(-(Z1.^2-2.*w1.*Z1.*Z2+Z2.^2)/(2.*(1-w1.^2))); %integrand
PrE_fun =@(Z) 1/sqrt(2.*pi).*exp(-Z.^2/2); %integrand for Pr(E_{t-1})
BetaLimit =@(k,L) integral2(fun,-L ,L ,-L-k*sqrt(n)/g, L-k*sqrt(n)/g )/integral(PrE_fun,-L,L);

%% Main code
% Calculates a table of average run lengths for sigma calculated from
% sample subgroups for multiple values of k and L provided by the user.
% (Table 2).
for k = 1:knum
for L = 1:Lnum
ARLmat(k,L)=1./(1-BetaLimit(kvec(k),Lvec(L)));
end
end

ARLmat