qcc::oc.curves.xbar

用例

该函数用于绘制 xbar 图的 OC 曲线,用例如下:

• 将数据按样本分组.数据原先为向量,并带有标记其样本号的变量. qcc.groups() 将数据按其样本号分为一个个子组,并组成矩阵.

library(qcc)

## Loading required package: MASS
## Package 'qcc', version 2.2
## Type 'citation("qcc")' for citing this R package in publications.

data(pistonrings)
head(pistonrings)

##   diameter sample trial
## 1    74.03      1  TRUE
## 2    74.00      1  TRUE
## 3    74.02      1  TRUE
## 4    73.99      1  TRUE
## 5    74.01      1  TRUE
## 6    74.00      2  TRUE

diameter <- qcc.groups(pistonrings$diameter, pistonrings$sample)
head(diameter)

##    [,1]  [,2]  [,3]  [,4]  [,5]
## 1 74.03 74.00 74.02 73.99 74.01
## 2 74.00 73.99 74.00 74.01 74.00
## 3 73.99 74.02 74.02 74.00 74.00
## 4 74.00 74.00 73.99 74.02 74.01
## 5 73.99 74.01 74.02 73.99 74.01
## 6 74.01 73.99 74.00 73.98 73.99

• 建立控制图对象,但不将图像绘制出,而是将其传递给 oc.curves(). 由其类型, oc.curves() 会自动调用 oc.curves.xbar(), 并绘制 OC 曲线.

cc <- qcc(diameter, type = "xbar", nsigmas = 3, plot = FALSE)
beta <- oc.curves(cc, identify = TRUE)

函数功能

其核心功能为计算 beta 的值,而 beta 应该由正态分布某些点的概率密度得出.具体见注释.

function(object, n, c = seq(0, 5, length = 101), nsigmas = object$nsigmas, identify = FALSE, 
    restore.par = TRUE) {
    # 判断输入对象类型
    type <- object$type
    if (!(object$type == "xbar")) 
        stop("not a `qcc' object of type \"xbar\".")


    # 提取组内样本量序列（qcc()允许每组样本量不同），排除样本量不等的情况
    size <- unique(object$sizes)
    if (length(size) > 1) 
        stop("Operating characteristic curves available only for equal sample sizes!")


    # 查看给定的样本量变化范围， 如果未给定，则用c(1, 5, 10, 15,
    # 20)代替（当然当前样本量必须要计算）
    if (missing(n)) 
        n <- unique(c(size, c(1, 5, 10, 15, 20)))


    # beta 计算 （Montgonery eq:6.19） 默认计算0.05为跨度的 sigma
    # 101个，从0sigma到5sigma
    beta <- matrix(NA, length(n), length(c))
    for (i in 1:length(n)) beta[i, ] <- pnorm(nsigmas - c * sqrt(n[i])) - pnorm(-nsigmas - 
        c * sqrt(n[i]))
    rownames(beta) <- paste("n=", n, sep = "")
    colnames(beta) <- c


    # 绘图
    oldpar <- par(bg = qcc.options("bg.margin"), cex = qcc.options("cex"), no.readonly = TRUE)
    if (restore.par) 
        on.exit(par(oldpar))
    plot(c, beta[1, ], type = "n", ylim = c(0, 1), xlim = c(0, max(c)), xlab = "Process shift (std.dev)", 
        ylab = "Prob. type II error ", main = paste("OC curves for", object$type, 
            "chart"))
    rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = qcc.options("bg.figure"))
    for (i in 1:length(n)) lines(c, beta[i, ], type = "l", lty = i)


    # identify = TRUE 时可以交互式的在 OC 曲线图上选点，见下图。
    beta <- t(beta)
    names(dimnames(beta)) <- c("shift (std.dev)", "sample size")
    if (identify) {
        cs <- rep(c, length(n))
        betas <- as.vector(beta)
        labels <- paste("c=", formatC(cs, 2, flag = "-"), ": beta=", formatC(betas, 
            4, flag = "-"), ", ARL=", formatC(1/(1 - betas), 2, flag = "-"), 
            sep = "")
        i <- identify(cs, betas, labels, pos = 4, offset = 0.2)
        apply(as.matrix(labels[i$ind]), 1, cat, "\n")
    } else {
        legend(max(c), 1, legend = paste("n =", n), bg = qcc.options("bg.figure"), 
            lty = 1:length(n), xjust = 1, yjust = 1)
    }
    invisible(beta)
}