cusum
R2026bDetect small changes in mean using cumulative sum
Syntax
Description
[
returns the first index of the upper and lower cumulative sums of
iupper,ilower] = cusum(x)x that have drifted beyond five standard deviations
above and below a target mean, respectively. The minimum detectable mean shift
is one standard deviation. The function estimates the target mean and standard
deviations from the first 25 samples of x.
[
returns all the indices at which the upper and lower cumulative sums exceed the
control limit.iupper,ilower] = cusum(___,"all")
[ also returns the upper and lower
cumulative sums.iupper,ilower,uppersum,lowersum]
= cusum(___)
cusum(___) with no output arguments
plots the upper and lower cumulative sums normalized to one standard
deviation above and below the target mean.
Examples
Generate and plot a 100-sample random signal with a linear trend.
fnc = rand(1,100) + linspace(0,1,100); plot(fnc)
Use cusum with no output arguments to plot the upper and lower cumulative sums with default settings.
cusum(fnc)
Compute the mean and standard deviation of the first 25 samples. Apply cusum using these numbers as the target mean and the target standard deviation. Highlight the point where the cumulative sum drifts more than five standard deviations beyond the target mean. Set the minimum detectable mean shift to one standard deviation.
mfnc = mean(fnc(1:25)); sfnc = std(fnc(1:25)); cusum(fnc,5,1,mfnc,sfnc)
Generate a signal resembling motion about an axle that becomes unstable due to wear. Add white Gaussian noise of variance 1/9.
sz = 200; dr = airy(2,linspace(-14.9371,1.2,sz)); rd = dr + sin(2*pi*(1:sz)/5) + randn(1,sz)/3;
Plot the growing background drift and the resulting signal.
plot(dr) hold on plot(rd,".-") hold off
Find the mean and standard deviation if the drift is not present and there is no noise. Plot the ideal noiseless signal and its stable background.
id = 0.3*sin(2*pi*(1:sz)/20); st = id + sin(2*pi*(1:sz)/5); mf = mean(st)
mf = -3.8212e-16
sf = std(st)
sf = 0.7401
plot(id) hold on plot(st,".-") hold off
Use the CUSUM Control Chart to pinpoint the onset of instability. Assume that the system becomes unstable when the signal is three standard deviations beyond its ideal behavior. Specify a minimum detectable shift of one standard deviation.
cusum(rd,3,1,mf,sf)
Increase the minimum detectable shift to 1.2 standard deviations. Return all instances of unwanted drift.
cusum(rd,3,1.2,mf,sf,"all")
Every hole in golf has an associated "par" that indicates the expected number of strokes needed to sink the ball. Skilled players usually complete each hole with a number of strokes very close to par. It is necessary to play several holes and let scores accumulate before a clear winner emerges in a match.
Ben, Jen, and Ken play a full round, which consists of 18 holes. The course has an assortment of par-3, par-4, and par-5 holes. At the end of the game, the players tabulate their scores.
hole = 1:18;
par = [4 3 5 3 4 5 3 4 4 4 5 3 5 4 4 4 3 4];
nms = {"Ben";"Jen";"Ken"};
Ben = [4 3 4 2 3 5 2 3 3 4 3 2 3 3 3 3 2 3];
Jen = [4 3 4 3 4 4 3 4 4 4 5 3 4 4 5 5 3 3];
Ken = [4 3 4 3 5 5 4 4 4 4 5 3 5 4 5 4 3 5];
T = table(hole',par',Ben',Jen',Ken', ...
VariableNames=["hole";"par";nms])T = 18×5 table
hole par Ben Jen Ken
____ ___ ___ ___ ___
1 4 4 4 4
2 3 3 3 3
3 5 4 4 4
4 3 2 3 3
5 4 3 4 5
6 5 5 4 5
7 3 2 3 4
8 4 3 4 4
9 4 3 4 4
10 4 4 4 4
11 5 3 5 5
12 3 2 3 3
13 5 3 4 5
14 4 3 4 4
15 4 3 5 5
16 4 3 5 4
⋮
The winner of the round is the player whose lower cumulative sum drifts the most below par at the end. Compute the sums for the three players to determine the winner. Make every shift in mean detectable by setting a small threshold.
[~,~,~,Bensum] = cusum(Ben-par,1,1e-4,0);
[~,~,~,Jensum] = cusum(Jen-par,1,1e-4,0);
[~,~,~,Kensum] = cusum(Ken-par,1,1e-4,0);
plot([Bensum;Jensum;Kensum]')
legend(nms,Location="best")
Ben wins the round. Simulate their next game by adding or subtracting a stroke per hole at random.
Ben = Ben+randi(3,1,18)-2;
Jen = Jen+randi(3,1,18)-2;
Ken = Ken+randi(3,1,18)-2;
[~,~,~,Bensum] = cusum(Ben-par,1,1e-4,0);
[~,~,~,Jensum] = cusum(Jen-par,1,1e-4,0);
[~,~,~,Kensum] = cusum(Ken-par,1,1e-4,0);
plot([Bensum;Jensum;Kensum]')
legend(nms,Location="best")
