Skip to content
Snippets Groups Projects
Commit a656187d authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Adiciona script de graf de controle para S e bar_X.

parent 12a6050f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 239 additions and 0 deletions
scripts/03_S_distr.R 0 → 100644
+ 239
− 0
View file @ a656187d
#=======================================================================
# Gráfico de controle para bar_X e S.
#-----------------------------------------------------------------------
# Definições da sessão.
# Instalar o pacote IQCC.
# install.package("IQCC")
suppressMessages(library(IQCC))
library(latticeExtra)
# Trellis graphical style.
mycol <- c("#E41A1C", "#377EB8", "#4DAF4A",
"#984EA3", "#FF7F00", "#FFFF33")
ps <- list(box.rectangle = list(col = 1, fill = c("gray70")),
box.umbrella = list(col = 1, lty = 1),
dot.symbol = list(col = 1, pch = 19),
dot.line = list(col = "gray50", lty = 3),
plot.symbol = list(col = 1, cex = 0.8),
plot.line = list(col = 1),
plot.polygon = list(col = "gray95"),
superpose.line = list(col = mycol, lty = 1),
superpose.symbol = list(col = mycol, pch = 1),
superpose.polygon = list(col = mycol),
strip.background = list(col = c("gray80", "gray50")))
trellis.par.set(ps)
library(reshape)
#-----------------------------------------------------------------------
# Distribuição amostral de R e S. Qual é o melhor estimador?
# Diferentes tamanhos de amostra.
n <- c(3, 5, 7, 10, 15, 25)
# Para cada tamanho de amostra, são geradas 500 amostras e calculados R
# e S.
L <- lapply(n,
FUN = function(n) {
res <- replicate(500, {
c4 <- c4(n)
d2 <- d2(n)
x <- rnorm(n)
r <- diff(range(x))
s <- sd(x)
c(r = r/d2, s = s/c4)
})
res <- cbind(n = n, as.data.frame(t(res)))
return(res)
})
str(L)
L <- do.call(rbind, L)
library(reshape)
M <- melt(data = L, id.vars = "n",
variable_name = "estim")
str(M)
# Qual é o estimador mais preciso ou eficiente?
ecdfplot(~value | factor(n), groups = estim, data = M,
scales = "free",
auto.key = TRUE,
as.table = TRUE)
str(M)
# names(trellis.par.get())
# trellis.par.get()$superpose.line$col
with(M, tapply(value, list(estim, n), mean))
with(M, tapply(value, list(estim, n), sd))
aggregate(M$value[subscripts],
M$estim[subscripts],
FUN = function(x) c(mean(x), sd(x)))
a <- aggregate(cbind(y = M$value),
by = list(M$estim),
FUN = function(x) c(mean(x), sd(x)))
a
c(t(apply(a[, -1], 1, FUN = function(x) { x[1] + c(-1, 0, 1) * x[2]})))
xtabs(~estim + n, M)
densityplot(~value | factor(n), groups = estim, data = M,
# scales = "free",
auto.key = TRUE, as.table = TRUE,
plot.points = FALSE, ref = TRUE,
panel=function(x, subscripts, ...){
panel.densityplot(x, subscripts = subscripts, ...)
a <- aggregate(cbind(M$value[subscripts]),
list(estim = M$estim[subscripts]),
FUN = function(x) c(mean(x), sd(x)))
m <- c(t(apply(a[, -1], 1,
FUN = function(x) {
x[1] + c(-1, 0, 1) * x[2]
})))
k <- trellis.par.get(
"superpose.line")$col[1:nlevels(M$estim)]
panel.abline(v = m, col = k, lty = 2)
# panel.grid(v = 2)
})
bwplot(value ~ estim | factor(n), data = M, pch = "|",
as.table = TRUE, layout = c(NA, 1),
panel = function(x, y, subscripts, ...) {
panel.bwplot(x = x, y = y, subscripts = subscripts, ...)
m <- tapply(y, x, mean)
panel.points(x = 1:nlevels(x), y = m, pch = 19)
}) + layer(panel.abline(h = 1, lty = 2))
#-----------------------------------------------------------------------
# Gráfico de controle para n fixo.
# x <- read.table("clipboard")
# dput(x[, 1])
browseURL(paste0("http://bcs.wiley.com/he-bcs/Books?action=",
"chapter&bcsId=7327&itemId=1118146816&chapterId=80371"))
# arquivo: ch06.xlxs, sheet: Exercise Data, coluna: Ex6-15ID (AB).
x <- c(74.03, 74.002, 74.019, 73.992, 74.008, 73.995, 73.992, 74.001,
74.011, 74.004, 73.988, 74.024, 74.021, 74.005, 74.002, 74.002,
73.996, 73.993, 74.015, 74.009, 73.992, 74.007, 74.015, 73.989,
74.014, 74.009, 73.994, 73.997, 73.985, 73.993, 73.995, 74.006,
73.994, 74, 74.005, 73.985, 74.003, 73.993, 74.015, 73.988,
74.008, 73.995, 74.009, 74.005, 74.004, 73.998, 74, 73.99,
74.007, 73.995, 73.994, 73.998, 73.994, 73.995, 73.99, 74.004,
74, 74.007, 74, 73.996, 73.983, 74.002, 73.998, 73.997, 74.012,
74.006, 73.967, 73.994, 74, 73.984, 74.012, 74.014, 73.998,
73.999, 74.007, 74, 73.984, 74.005, 73.998, 73.996, 73.994,
74.012, 73.986, 74.005, 74.007, 74.006, 74.01, 74.018, 74.003,
74, 73.984, 74.002, 74.003, 74.005, 73.997, 74, 74.01, 74.013,
74.02, 74.003, 73.982, 74.001, 74.015, 74.005, 73.996, 74.004,
73.999, 73.99, 74.006, 74.009, 74.01, 73.989, 73.99, 74.009,
74.014, 74.015, 74.008, 73.993, 74, 74.01, 73.982, 73.984,
73.995, 74.017, 74.013)
length(x)
n <- rep(1:25, each = 5)
library(qcc)
da <- t(unstack(x = data.frame(i = n, x = x), form = x ~ i))
obj <- qcc(data = da, type = "xbar")
obj <- qcc(data = da, type = "S")
#-----------------------------------------------------------------------
# Gráfico de controle para n variável.
# n <- scan()
# dput(n)
n <- c(5, 3, 5, 5, 5, 4, 4, 5, 4, 5, 5, 5, 3, 5, 3, 5, 4, 5, 5, 3, 5, 5,
5, 5, 5)
# Cada amostra separada.
L <- do.call(rbind,
by(data = da, INDICES = da$i,
FUN = function(d) {
i <- d$i[1]
d[1:n[i], ]
}))
str(L)
# As estimativas de média e desvio-padrão por amostra.
A <- aggregate(x ~ i, data = L,
FUN = function(x) {
c(m = mean(x), s = sd(x))
})
str(A)
# Média das médias: bar_bar_X.
bbx <- sum(A$x[, "m"] * n)/sum(n)
# Médias dos desvios-padrões: bar_S.
bs <- sqrt(sum((n - 1) * A$x[, "s"]^2)/c(sum(n - 1)))
# As constantes que dependem do n.
A3 <- 3/(c4(n) * sqrt(n))
B3 <- 1 - 3 * sqrt(1 - c4(n)^2)/c4(n)
B4 <- 1 + 3 * sqrt(1 - c4(n)^2)/c4(n)
# Os limites para as duas medidas.
LC <- data.frame(i = 1:length(n),
LCm = bbx,
LICm = bbx - A3 * bs,
LSCm = bbx + A3 * bs,
LCs = bs,
LICs = pmax(0, B3 * bs),
LSCs = B4 * bs)
# bar_X chart.
matplot(x = LC$i - 0.5, y = LC[, c("LCm", "LICm", "LSCm")],
type = "n",
xlab = "Número da amostra",
ylab = expression(bar(x)))
lines(LC$i, LC$LCm)
with(LC, segments(x0 = i - 0.5, x1 = i + 0.5, y0 = LICm, y1 = LICm))
with(LC, segments(x0 = i - 0.5, x1 = i + 0.5, y0 = LSCm, y1 = LSCm))
points(x = LC$i, y = A$x[, "m"], type = "o")
# S chart.
matplot(x = LC$i - 0.5,
y = LC[, c("LCs", "LICs", "LSCs")],
type = "n",
xlab = "Número da amostra",
ylab = expression(s))
lines(LC$i, LC$LCs)
with(LC, segments(x0 = i - 0.5, x1 = i + 0.5, y0 = LICs, y1 = LICs))
with(LC, segments(x0 = i - 0.5, x1 = i + 0.5, y0 = LSCs, y1 = LSCs))
points(x = LC$i, y = A$x[, "s"], type = "o")
# DANGER: tem alguma coisa nesse código passo a passo, algum fator de
# escala, por exemplo, que faz com que esses resultados sejam
# ligeiramente diferentes do pacote qqc.
#-----------------------------------------------------------------------
# help(qcc, h = "html")
# qcc()
# help.search("rbind.fill")
sam <- by(L, L$i,
FUN = function(d) {
as.data.frame(matrix(d$x, nrow = 1))
})
db <- do.call(plyr::rbind.fill, sam)
qcc(data = db, type = "xbar")
qcc(data = db, type = "S")
#-----------------------------------------------------------------------
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment