#!/usr/bin/perl -w
use strict;
use Math::Random;
my %shiftSchedule = (
"first" => { "start" => 6.00, "end" => 14.00 },
"second" => { "start" => 14.00, "end" => 22.00 },
"third" => { "start" => 22.00, "end" => 6.00 }
);
my $shift = "third"; # shift to monitor
my $inspectionRate = 1 / 2; # every 1/2 hour
my $drift = 1.5; # drift to simulate
my $m = 25; # samples in control chart setup
my $n = 12; # observations per sample
my $target = 100.0; # quality characteristic target
my $hour;
my $i;
my $j;
my $minute;
my $observation;
my $setupM = $m;
print "timestamp sample observation phase\r\n";
for ($i = 1; $i <= $m; $i++) {
for ($j = 0; $j < $n; $j++) {
$observation = $target + random_normal();
printf " 0:00 %6d %7.3f setup\r\n", $i, $observation;
}
}
$m = $shiftSchedule{$shift}{"end"} - $shiftSchedule{$shift}{"start"};
if ($m < 0) {
$m += 24;
}
$m /= $inspectionRate;
for ($i = 1; $i <= $m; $i++) {
$hour = int($i * $inspectionRate + $shiftSchedule{$shift}{"start"});
if ($hour >= 24) {
$hour -= 24;
}
$minute = ($i & 0x1) ? (60 * $inspectionRate) : 0;
for ($j = 0; $j < $n; $j++) {
$observation = $target + random_normal();
if ($i >= (0.25 * $m)) {
if ($i < (0.75 * $m)) {
$observation += ($drift / (0.5 * $m)) * ($i - (0.25 * $m));
} else {
$observation += $drift;
}
}
printf " %2d:%02d %6d %7.3f monitoring\r\n", $hour, $minute, $setupM + $i, $observation;
}
}
observations <- read.table("observations.txt", TRUE)
require(qcc)
attach(observations)
# identify "observation" and "sample" columns as providing the
# observations and sample numbers
observation <- qcc.groups(observation, sample)
# number of observations per sample
n = ncol(observation)
# number of samples in phase I ("setup")
m = length(phase[phase == "setup"]) / n
# we do not plot the phase I ("setup") data, so just provide a filler
setupTimestamps <- rep.int(c(""), m)
# extract the timestamps to display for the phase II ("monitoring") data
monitoringTimestamps <- as.character(timestamp[phase == "monitoring"])
# reduce from one per observation to one per sample
monitoringTimestamps <- monitoringTimestamps[seq(1, length(monitoringTimestamps), n)]
# reduce to one per hour
for (i in 1:length(monitoringTimestamps)) {
minutes <- strsplit(monitoringTimestamps[[i]], ":")[[1]][[2]]
if (minutes != "00") monitoringTimestamps[[i]] <- ""
}
# plot xbar chart
obj <- qcc(data = observation[1:m,],
type = "xbar",
newdata = observation[-(1:m),],
labels = setupTimestamps,
newlabels = monitoringTimestamps,
axes.las = 3,
chart.all = FALSE,
title = "xbar chart for quality characteristic XXX",
xlab = "Sample",
ylab = "Mean value (units)")