# Inspeção por Amostragem - Atributos
#
# Escolha do Plano de Amostragem

# Walter Ponge-Ferreira
# São Paulo, 05.06.2018

setwd("~/Documents/Data/Pasta Tecnica/R/PME3463 Aula 07")
rm(list=ls())

# PRP - Producer Risk Point
Alpha <- 0.05 # 5%
P0 <- 0.01 # 1%

# CRP - Consumer Risk Point
Beta <- 0.10 # 10%
P1 <- 0.05 # 5%

# Distribuição Binomial
# X ~ bin(n,p)
#

# Busca do plano de amostragem por tentativa e erro
# n = ?
# Ac = ?

# Primeira Tentativa
n <- 120
Ac <- 5
p <- seq(0,0.1,by=0.001)
Pa <- pbinom(q = Ac, size = n, prob = p)
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
     main=leg,
     xlab = "Proporção de defeituosas",
     ylab = "Probabilidade de Aceitação do Lote")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")

# Segunda Tentativa
n <- 120
Ac <- 3
p <- seq(0,0.1,by=0.001)
Pa <- pbinom(q = Ac, size = n, prob = p)
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
     main="Curva Característica para n = 120 e Ac = 3",
     xlab = "Proporção de defeituosas",
     ylab = "Probabilidade de Aceitação do Lote")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")

# Terceira Tentativa
n <- 120
Ac <- 2
p <- seq(0,0.1,by=0.001)
Pa <- pbinom(q = Ac, size = n, prob = p)
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
     main= leg,
     xlab = "Proporção de defeituosas",
     ylab = "Probabilidade de Aceitação do Lote")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")

# Quarta Tentativa
n <- 150
Ac <- 3
p <- seq(0,0.1,by=0.001)
Pa <- pbinom(q = Ac, size = n, prob = p)
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
    main = leg,
    xlab = "Proporção de defeituosas",
    ylab = "Probabilidade de Aceitação do Lote")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")

# Quinta Tentativa
n <- 130
Ac <- 3
p <- seq(0,0.1,by=0.001)
Pa <- pbinom(q = Ac, size = n, prob = p)
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
     main = leg,
     xlab = "Proporção de defeituosas",
     ylab = "Probabilidade de Aceitação do Lote")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")

# Busca da solução ótima

library(AcceptanceSampling)

(meuplano <- OC2c(n = n, c = Ac, pd = seq(0, 0.1, 0.001)))
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
plot(meuplano, xlim = c(0, 0.10), type = "l", lwd = 2,  las = 1,
     main = leg)
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")
text(P1, Beta, "CRP", pos = 3, col = "red")
assess(meuplano, PRP = c(P0,1-Alpha), CRP = c(P1,Beta))

(melhorparametro <- find.plan(PRP = c(0.01,0.95), CRP = c(0.05,0.10)))


melhorplano <- OC2c(n = melhorparametro$n, c = melhorparametro$c, r = melhorparametro$r, pd = seq(0, 0.1, 0.001))
leg <- paste("Curva Característica para n = ",melhorparametro$n," e Ac = ",melhorparametro$c)
plot(melhorplano, xlim = c(0, 0.10), type = "l", lwd= 2, las = 1,
     main = leg)
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")
text(P1, Beta, "CRP", pos = 3, col = "red")
assess(melhorplano, PRP = c(P0,1-Alpha), CRP = c(P1,Beta))

# Avalia Solução Ótima

p <- seq(0, 0.1, by = 0.001)
Pa <- pbinom(q = Ac, size = n, prob = p)
leg <- paste("Curva Característica para n = ",melhorparametro$n," e Ac = ",melhorparametro$c)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
     main = leg,
     xlab = "Proporção de Defeituosas",
     ylab = "Probabilidade de Aceitação")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")
text(P1, Beta, "CRP", pos = 3, col = "red")

# Qualidade Média Resultante - QMR
AOQ <- Pa * p # N >> n
plot(AOQ ~ p, type = "l", lwd = 2, las = 1, xlab = "Proporção de Defeituosas")
grid()
(AOQL <- max(AOQ))
abline(h = AOQL, col = "red", lty = 2)
abline(v = P0, col = "green", lty = 2)
text(8*P0,AOQL,paste("AOQL = ",round(AOQL*100, digits = 2),"%"), pos = 1, col = "red")
text(P0, 0, paste("AQL = ",round(P0*100, digits = 2),"%"), pos = 4, col = "green")
title(main = "Qualidade Média Resultante - QMR")

# ------------------------------------------------
# Plano de Amostragem pela NBR 5426
#

# Palno de Amostragem Simples
# Tamanho do Lote
N <- 300
# Nível de Inspeção Normal (II)
# Da tabela 1: Código literal -> H

# Limite de Qualidade Aceitavel
# AQL = 1%
# Da tabela 2:
# n = 50
# Ac = 1
# Re = 2
n <- 50
Ac <- 1
Re <- 2

p <- seq(0, 0.1, by = 0.001)
leg <- paste("Curva Característica para n = ",n," e Ac = ",Ac)
Pa <- pbinom(q = Ac, size = n, prob = p)
plot(Pa ~ p, type = "l", lwd = 2, las = 1,
     main = leg,
     xlab = "Proporção de Defeituosas",
     ylab = "Probabilidade de Aceitação",
     sub = "Norma NBR 5426 - Plano Simples - Nivel Normal")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")
text(P1, Beta, "CRP", pos = 3, col = "red")

AOQ <- Pa * p
plot(AOQ ~ p, type = "l", lwd = 2, las = 1, xlab = "Proporção de Defeituosas")
grid()
(AOQL <- max(AOQ))
abline(h = AOQL, col = "red", lty = 2)
abline(v = P0, col = "green", lty = 2)
text(8*P0,AOQL,paste("AOQL = ",round(AOQL*100, digits = 2),"%"), pos = 1, col = "red")
text(P0, 0, paste("AQL = ",round(P0*100, digits = 2),"%"), pos = 4, col = "green")
title(main = "Qualidade Média Resultante - QMR",
      sub = "Norma NBR 5426 - Plano Simples - Nível Normal")

# Avalia Plano recomendado pela NBR 5426

planoNBR <- OC2c(n, Ac, Re, pd = seq(0, 0.1, 0.001))
plot(planoNBR, xlim = c(0, 0.10), type = "l", lwd = 2, las = 1,
     main = "CCO para Nível Normal - Código Literal H",
     sub = "Norma NBR 5426 - Plano Simples - Nivel Normal")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")
text(P1, Beta, "CRP", pos = 3, col = "red")

assess(planoNBR, PRP = c(P0,1-Alpha), CRP = c(P1,Beta))

AOQ <- planoNBR@paccept * planoNBR@pd
plot(AOQ ~ planoNBR@pd, xlim = c(0,0.10),type = "l", lwd = 2, las = 1,
     xlab = "Proporção de Defeituosas",
     ylab = "AQL")
grid()
(AOQL <- max(AOQ))
abline(h = AOQL, col = "red", lty = 2)
abline(v = P0, col = "green", lty = 2)
text(8*P0,AOQL,paste("AOQL = ",round(AOQL*100, digits = 2),"%"), pos = 1, col = "red")
text(P0, 0, paste("AQL = ",round(P0*100, digits = 2),"%"), pos = 4, col = "green")
title(main = "Qualidade Média Resultante - QMR",
      sub = "Norma NBR 5426 - Plano Simples - Nível Normal")

# Palno de Amostragem Duplo
# Tamanho do Lote
N <- 300
# Nível de Inspeção Normal (II)
# Da tabela 1: Código literal -> H

# Limite de Qualidade Aceitavel
# AQL = 1%
# Da tabela 5:
# n1 = 31
# Ac1 = 0
# Re1 = 2
# n2 = 32
# Ac2 = 1
# Re2 = 2
n <- c(32,32)
Ac <- c(0,1)
Re <- c(2,2)

# Avalia Plano recomendado pela NBR 5426

planoNBR <- OC2c(n, c = Ac, r = Re, type = "binomial", pd = seq(0, 0.1, 0.001))
plot(planoNBR, xlim = c(0, 0.10), type = "l", lwd = 2, las = 1,
     main = "CCO para Nível Normal - Código Literal H",
     sub = "Norma NBR 5426 - Plano Duplo - Nível Normal")
grid()
points(c(P0,P1),c(1-Alpha,Beta), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")
text(P1, Beta, "CRP", pos = 3, col = "red")

assess(planoNBR, PRP = c(P0,1-Alpha), CRP = c(P1,Beta))

AOQ <- planoNBR@paccept * planoNBR@pd
plot(AOQ ~ planoNBR@pd, xlim = c(0,0.10),type = "l", lwd = 2, las = 1,
     xlab = "Proporção de Defeituosas",
     ylab = "AOQ")
grid()
(AOQL <- max(AOQ))
abline(h = AOQL, col = "red", lty = 2)
abline(v = P0, col = "green", lty = 2)
text(8*P0,AOQL,paste("AOQL = ",round(AOQL*100, digits = 2),"%"), pos = 1, col = "red")
text(P0, 0, paste("AQL = ",round(P0*100, digits = 2),"%"), pos = 4, col = "green")
title(main = "Qualidade Média Resultante - QMR",
      sub = "Norma NBR 5426 - Plano Duplo - Nível Normal")

# Palno de Amostragem Multiplo
# Tamanho do Lote
N <- 300
# Nível de Inspeção Normal (II)
# Da tabela 1: Código literal -> H

# Limite de Qualidade Aceitavel
# AQL = 6,5% # Mudança de especificação: vide observação abaixo na definição do Ac
# Da tabela 8: (7 amostras)
# n = 13; 13; 13; 13; 13; 13; 13
# Ac = 0;  1;  3;  5;  7; 10; 13
# Re = 4;  6;  8; 10; 11; 12; 14
P0 <- 0.065
Alpha = 0.05

n <- rep(13,7)
Ac <- c(0,1,3,5,7,10,13) # Atenção: função não aceita valores abaixo de 0!
Re <- c(4,6,8,10,11,12,14)

# Avalia Plano recomendado pela NBR 5426

planoNBR <- OC2c(n, c = Ac, r = Re, type = "binomial", pd = seq(0, 10*P0, P0/10))
plot(planoNBR, xlim = c(0, 10*P0), type = "l", lwd = 2, las = 1,
     main = "CCO para Nível Normal - Código Literal H",
     sub = "Norma NBR 5426 - Plano Múltipla - Nível Normal")
grid()
legend("topright", title = "Plano de Amostragem", paste("AQL = ",round(P0*100, digits = 2),"%"))
points(c(P0),c(1-Alpha), pch = 16, col = "red")
text(P0, (1-Alpha), "PRP", pos = 1, col = "red")

AOQ <- planoNBR@paccept * planoNBR@pd
plot(AOQ ~ planoNBR@pd, xlim = c(0,10*P0),type = "l", lwd = 2, las = 1,
     xlab = "Proporção de Defeituosas",
     ylab = "AOQ")
grid()
legend("right", title = "Plano de Amostragem",paste("AQL = ",round(P0*100, digits = 2),"%"))
(AOQL <- max(AOQ))
abline(h = AOQL, col = "red", lty = 2)
abline(v = P0, col = "green", lty = 2)
text(8*P0,AOQL,paste("AOQL = ",round(AOQL*100, digits = 2),"%"), pos = 1, col = "red")
text(P0, 0, paste("AQL = ",round(P0*100, digits = 2),"%"), pos = 4, col = "green")
title(main = "Qualidade Média Resultante - QMR",
      sub = "Norma NBR 5426 - Plano Múltiplo - Nível Normal")

# Fim do Arquivo