Lista 1 - Controle Estatístico da qualidade
Arthur Machado; Willams Batista; Heitor Victor
31/10/2019
Pacotes utilizados
library(qcc)
library(nortest)
library(readr)
library(readxl)Exercícios
1)
Medidas de viscosidade de um polímero são feitas a cada 10 minutos por um viscômetro on-line. Trinta e seis observações são mostradas na tabela abaixo (leia de cima para baixo e da esquerda para a direita). O alvo da viscosidade para esse processo é \(\mu_0 = 3200\).
#lendo o banco de dados
data<-data.frame(c(3169, 3173, 3162, 3154, 3139, 3145, 3160, 3172, 3175,
3205, 3203, 3209, 3208, 3211, 3214, 3215, 3209, 3203,
3185, 3187, 3192, 3199, 3197, 3193, 3190, 3183, 3197,
3188, 3183, 3175, 3174, 3171, 3180, 3179, 3175, 3174))
colnames(data)<-paste(c("viscosidade"))a)
Estime o desvio padrão do processo;
sd_process<-sd(data$viscosidade)b)
Construa um gráfico CUSUM tabular para esse processo, usando valores de \(h = 8.01\) e \(k = 0.25\);
h<-8.01;
k<-0.25
mu0<-3200
sig0<-sd_process
mu1<-mu0 + k*sig0
q <- cusum(data$viscosidade, decision.interval = h, se.shift = k, center=mu0, std.dev=sig0)
summary(q)##
## Call:
## cusum(data = data$viscosidade, center = mu0, std.dev = sig0, decision.interval = h, se.shift = k)
##
## cusum chart for data$viscosidade
##
## Summary of group statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3139 3174 3184 3185 3200 3215
##
## Group sample size: 1
## Number of groups: 36
## Center of group statistics: 3200
## Standard deviation: 19.26062
##
## Decision interval (std.err.): 8.01
## Shift detection (std. err.): 0.25c)
Discuta a escolha de h e k na parte (b) em relação ao desempenho do gráfico CUSUM;
Temos que para essas especificações de h e k na questão b, temos um comprimento médio da sequência de 370 obsevações até o primeiro alarme falso ocorra, ou seja, nessas especificações esse gráfico CUSUM tabular fica pouco proveitoso no seu uso, pois apresenta o mesmo desempenho do gráfico de Shewhart.
Podemos provar isso, atráves dos códigos abaixo:
b = h + 1.166
#delta <- abs(mu1-mu0)/sig0
#Dado o processo estar sob controle, como delta=0 então CMS_mais = CMS_menos
delta <-0
Delta0 = -delta - k
Delta1 = delta - k
CMS0_mais = (exp(-2*Delta0*b) + 2*Delta0*b - 1)/(2*Delta0^2)
CMS0_menos = (exp(-2*Delta1*b) + 2*Delta1*b - 1)/(2*Delta1^2)
CMS0 = 1 / ((1/CMS0_mais) + (1/CMS0_menos))
CMS0## [1] 370.8386d)
Estabeleça e aplique um gráfico de controle MMEP usando \(λ = 0.1\) e \(L = 2.7\);
#MMEP - Usando o pacote qcc
lam<-0.1;
L<-2.7;
q <- ewma(data$viscosidade, lambda= lam, nsigmas= L, center=mu0, std.dev=sig0)
summary(q)##
## Call:
## ewma(data = data$viscosidade, center = mu0, std.dev = sig0, lambda = lam, nsigmas = L)
##
## ewma chart for data$viscosidade
##
## Summary of group statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3139.000 3173.750 3184.000 3184.667 3200.000 3215.000
##
## Group sample size: 1
## Number of groups: 36
## Center of group statistics: 3200
## Standard deviation: 19.26062
##
## Smoothing parameter: 0.1
## Control limits:
## LCL UCL
## [1,] 3194.800 3205.200
## [2,] 3193.004 3206.996
## ...
## [36,] 3188.073 3211.927e)
Discuta a escolha de λ e L na parte (d) em relação ao desempenho do gráfico MMEP.
Podemos dizer, que na escolha desse parâmetros não cai muito bem para o desempenho desse gráfico, dado que, podemos perceber que apresenta um comprimento médio da sequência baixo, pois ele aponta alarmes falsos em poucas amostras observadas, não conseguindo otimizar o processo na detecção de ‘erros’ desnecessários no processo, sendo inferior ao gráfico de Shewhart.
Isso pode ser mostrado pelos códigos abaixo.
#MMEP - Usando o pacote qcc
rep = 10000
n = 4000
resultado = NULL
for(j in 1:rep){
X = NULL
X = rnorm(n,mu1,sig0)
z = LIC = LSC = NULL
for(i in 1:length(X))
{ z[i] = lam*sum((1-lam)^seq(0,i-1)*X[i:1])+(1-lam)^i*mu0
LIC[i] = mu0 - L*sig0*sqrt(lam/(2-lam)*(1-(1-lam)^(2*i)))
LSC[i] = mu0 + L*sig0*sqrt(lam/(2-lam)*(1-(1-lam)^(2*i)))
if((z[i]>LSC[i]) || (z[i]<LIC[i]))
{ resultado[j] = i
break
}
}
}
CMS_MMEP = mean(resultado)
CMS_MMEP## [1] 84.46922)
Simule 20 valores de X~N(10,1) e outros 50 valores para X~N(9,1). Verifique qual dos dois gráficos: CUSUM(\(h = 4.774\); \(k = 0.5\)) ou MMEP(\(λ = 0.2\); \(L = 2.859\)) sinaliza com maior rapidez o deslocamento da média do processo (de 10 para 9). Avalie também o CMS de cada gráfico gerado.
#Dados
set.seed(123)
X=NULL
X[1:20] = rnorm(20,10,1)
X[21:70] = rnorm(50,9,1)
h<-4.774
k<-0.5
sig0<-sd(X)
mu0<-10
#mu1<-mu0 + k*sig0
mu1<-9Utilizando o pacote qcc do R. No gráfico de cusum tabular, podemos perceber que depois da amostra de número 20 (com processo ainda sob controle) o gráfico tinha alertado que o processo tinha saído do controle na amostra de número 23-24, ou seja, quando foi submetido uma pequena mudança na média do processo (mudança em uma unidade), o gráfico cusum foi extremamente eficiente.
#CUMSUM - Usando o pacote qcc
q <- cusum(X, decision.interval = h, se.shift = k, center=mu0, std.dev=sig0)
summary(q)##
## Call:
## cusum(data = X, center = mu0, std.dev = sig0, decision.interval = h, se.shift = k)
##
## cusum chart for X
##
## Summary of group statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.313 8.674 9.346 9.360 10.033 11.787
##
## Group sample size: 1
## Number of groups: 70
## Center of group statistics: 10
## Standard deviation: 1.034308
##
## Decision interval (std.err.): 4.774
## Shift detection (std. err.): 0.5#Comprimento médio da sequência
b = h + 1.166
delta = abs(mu1-mu0)/sig0
#Dado o processo estar sob controle, como delta=0 então CMS_mais = CMS_menos
#delta <-0
Delta0 = -delta - k
Delta1 = delta - k
CMS0_mais = (exp(-2*Delta0*b) + 2*Delta0*b - 1)/(2*Delta0^2)
CMS0_menos = (exp(-2*Delta1*b) + 2*Delta1*b - 1)/(2*Delta1^2)
CMS0_CUMSUM = 1 / ((1/CMS0_mais) + (1/CMS0_menos))
CMS0_CUMSUM## [1] 10.43874Utilizando o pacote qcc, podemos fazer o gráfico de médias móveis exponencialmente ponderadas. Nesse gráfico, podemos perceber que após a amostra de número 20 (com processo ainda sob controle), foi alertado que o processo tinha saído do controle na amostra de número 24, ou seja, em uma pequena mudança na média, em apenas 4 amostras o gráfico conseguiu detectar essa mudança. Sendo superior ao gráfico de Shewhart.
#MMEP - Usando o pacote qcc
lam<-0.2; L<-2.859;
q <- ewma(X, lambda= lam, nsigmas= L, center=mu0, std.dev=sig0)
summary(q)##
## Call:
## ewma(data = X, center = mu0, std.dev = sig0, lambda = lam, nsigmas = L)
##
## ewma chart for X
##
## Summary of group statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.313307 8.673604 9.346393 9.359541 10.033448 11.786913
##
## Group sample size: 1
## Number of groups: 70
## Center of group statistics: 10
## Standard deviation: 1.034308
##
## Smoothing parameter: 0.2
## Control limits:
## LCL UCL
## [1,] 9.408583 10.59142
## [2,] 9.242617 10.75738
## ...
## [70,] 9.014305 10.98570#Comprimento médio da sequência
rep = 10000
n = 4000
resultado = NULL
for(j in 1:rep){
X = NULL
X = rnorm(n,mu1,sig0)
z = LIC = LSC = NULL
for(i in 1:length(X))
{ z[i] = lam*sum((1-lam)^seq(0,i-1)*X[i:1])+(1-lam)^i*mu0
LIC[i] = mu0 - L*sig0*sqrt(lam/(2-lam)*(1-(1-lam)^(2*i)))
LSC[i] = mu0 + L*sig0*sqrt(lam/(2-lam)*(1-(1-lam)^(2*i)))
if((z[i]>LSC[i]) || (z[i]<LIC[i]))
{ resultado[j] = i
break
}
}
}
CMS_MMEP = mean(resultado)Comparando ambos CMS, podemos dizer que ambos os gráficos foram eficazes em detectar pequenas mudanças na média do processo, apresentando resultados bem similares, porém, na escolha de uma, o gráfico de médias móveis exponencialmente ponderadas foi levemente superior na detecção de pequenas mudanças.
cbind(CMS0_CUMSUM,CMS_MMEP)## CMS0_CUMSUM CMS_MMEP
## [1,] 10.43874 9.41863)
Simule 20 valores de X~N(10,1) e outros \(50\) valores para X~N(7,1). Verifique qual dos dois gráficos: Shewhart(\(L = 3\)) ou MMEP(\(λ = 0,2; L = 2,859\)) sinaliza com maior rapidez o deslocamento da média do processo (de 10 para 7). Avalie também o CMS de cada gráfico gerado.
#Dados
set.seed(123)
X=NULL
X[1:20] = rnorm(20,10,1)
X[21:70] = rnorm(50,7,1)
sig0<-sd(X)
mu0<-10
mu1<-7Gráfico de médias móveis exponencialmente ponderadas utilizando o pacote qcc.
#MMEP - Usando o pacote qcc
lam<-0.2; L<-2.859;
q <- ewma(X, lambda= lam, nsigmas= L, center=mu0, std.dev=sig0)
summary(q)##
## Call:
## ewma(data = X, center = mu0, std.dev = sig0, lambda = lam, nsigmas = L)
##
## ewma chart for X
##
## Summary of group statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.313307 6.696760 7.501064 7.930970 9.277099 11.786913
##
## Group sample size: 1
## Number of groups: 70
## Center of group statistics: 10
## Standard deviation: 1.674728
##
## Smoothing parameter: 0.2
## Control limits:
## LCL UCL
## [1,] 9.042391 10.95761
## [2,] 8.773662 11.22634
## ...
## [70,] 8.403985 11.59602No gráfico EWMA podemos dizer que a escolha desses parâmetros foi bastante eficaz, pois ele conseguiu identficar uma mudança na média na amostra de número 23, ou seja, após apenas 3 amostras que sairam de controle ele conseguiu identificar o afastamento da média em 3 unidades.
Gráfico de Shewart utilizando o pacote qcc.
L<- 3
q <- qcc(X, type="xbar.one", center=mu0, std.dev=sig0, nsigmas=L)
summary(q)##
## Call:
## qcc(data = X, type = "xbar.one", center = mu0, std.dev = sig0, nsigmas = L)
##
## xbar.one chart for X
##
## Summary of group statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.313307 6.696760 7.501064 7.930970 9.277099 11.786913
##
## Group sample size: 1
## Number of groups: 70
## Center of group statistics: 10
## Standard deviation: 1.674728
##
## Control limits:
## LCL UCL
## 4.975817 15.02418No gráfico de Shewart, o processo conseguiu identificar a sua primeira observação fora de controle a partir da amostra de número 26. Porém, todas as observações ficaram dentro dos limites de controle, apenas apresentandos duas flutuações, uma das 20 primeiras amostras e outra após essas 20 amostras, dando uma ideia que para esses limites 3\(\sigma\) a capacidade de identificar pequenas mudanças é ineficaz, dado que essas pequenas flutuações da média são ignoradas (sem pontos vermelhos) no gráfico.
É sensato dizer que as escolhas dos parâmetros desse gráfico EWMA foi eficaz, sendo levemente superior ao gráfico de SHEWART na detecção de afastamentos da média.
Esses resultados podem ser constatados abaixo:
# CMS SHEWHART
beta = pnorm(mu0+3*sig0, mu1,sig0) - pnorm(mu0-3*sig0, mu1,sig0)
alpha<-0.0027
CMS1_SHEWHART<-1/(1-beta)
L<- 3
lam<-0.2
#Comprimento médio da sequência - MMEP
rep = 10000
n = 4000
resultado = NULL
for(j in 1:rep){
X = NULL
X = rnorm(n,mu1,sig0)
z = LIC = LSC = NULL
for(i in 1:length(X))
{ z[i] = lam*sum((1-lam)^seq(0,i-1)*X[i:1])+(1-lam)^i*mu0
LIC[i] = mu0 - L*sig0*sqrt(lam/(2-lam)*(1-(1-lam)^(2*i)))
LSC[i] = mu0 + L*sig0*sqrt(lam/(2-lam)*(1-(1-lam)^(2*i)))
if((z[i]>LSC[i]) || (z[i]<LIC[i]))
{ resultado[j] = i
break
}
}
}
CMS_MMEP = mean(resultado)A capacidadede identificar que a média do processo saiu de controle do gráfico MMEP é superior pois precisa de apenas 3.4888 amostras em média para identificar essa mudança. Já pelo gráfico de Shewart precisa de 8.818588 amostras em média para relatar sobre essa mudança.
cbind(CMS1_SHEWHART, CMS_MMEP)## CMS1_SHEWHART CMS_MMEP
## [1,] 8.818588 3.48884)
Suponha que um produto seja embalado em lotes de tamanho \(\mathcal{N} = 500\). O procedimento de inspeção consiste em uma única amostragem de \(n = 200\) itens, adotado-se Ac = 20.
N<-500
n<-200
Ac<-20a)
Calcule o risco do produtor para \(p_0 = 0.05\) e o risco do consumidor para \(p_1 = 0.15\) com base na distribuição binomial;
p0<-0.05
p1<-0.15
risco_produto<-pbinom(20,n,prob = p0)
risco_produto## [1] 0.9988401risco_consumidor<-pbinom(20,n,prob = p1)
risco_consumidor## [1] 0.02547343b)
Calcule o risco do produtor para \(p_0 = 0.05\) e o risco do consumidor para \(p_1 = 0.15\) com base na distribuição hipergeométrica;
p0<-0.05
p1<-0.15
risco_produto<-phyper(20,p0*N, N-p0*N, n)
risco_produto## [1] 0.9999951risco_consumidor<-phyper(20,p1*N, N-p1*N, n)
risco_consumidor## [1] 0.006743276c)
Obtenha as curvas características de operação baseadas nas duas distribuições.
p<-seq(0,0.40,0.01)
pa1<-phyper(20,p*N, N-p*N, n)
pa2<-pbinom(20,n,prob = p)
plot(p,pa1, ylim=c(0,1), type="l", col="blue")
lines(pa1)
par(new=TRUE)
plot(p,pa2, ylim=c(0,1), type="l", col="red")
lines(pa2)
legend("topright",legend=list("Hipergeométrica","Binomial"),fill = c("blue","red"), col=c("blue","red"))
5)
Suponha que um plano de amostragem única com \(n = 150\) e \(Ac = 2\) tenha sido aplicado para lotes de tamanho \(N = 3000\).
n<-150
Ac<-2
N<-3000a)
Apresente a curva característica de operação para esse plano;
p<-seq(0,0.40,0.01)
pa1<-phyper(2,p*N, N-p*N, n)
plot(p,pa1, ylim=c(0,1.1), type="l", col="blue", ylab="Probabilidade de Aceitação", xlab="Fração de defeituosos do lote")
lines(pa1)
legend("topright",legend=c("Hipergeometrica"),fill = c("blue"), col=c("blue"))
b)
Considerando \(p = 0.02\), calcule a Qualidade Média de Saída (QMS) e a Inspeção Total Média (ITM) para o plano apresentado ao adotar inspeção retificadora;
p<-0.02
pa1<-phyper(Ac,p*N, N-p*N, n)- Qualidade de Saída Média (QSM): mede a qualidade no lote resultante da aplicação da inspeção de retificação, sendo definida por:
QMS<-((pa1*(N-n))/(N))*p
QMS## [1] 0.007889343Como o \(QMS < p\), então a qualidade do lote está em boas condições.
- Inspeção Total Média (ITM): mede o número médio de itens inspecionados, devido ao uso de um programa de inspeção por retificação.
ITM <- n + (1 - pa1)*(N - n)
ITM## [1] 1816.599Note que o \(n < ITM < N\), isso é bom.
c)
Considerando \(p = 0.05\), calcule a Qualidade Média de Saída (QMS) e a Inspeção Total Média (ITM) para o plano apresentado ao adotar inspeção retificadora;
p<-0.05
pa1<-phyper(Ac,p*N, N-p*N, n)
QMS<-((pa1*(N-n))/(N))*p
QMS## [1] 0.0007692877Como o \(QMS < p\), então a qualidade do lote está em boas condições.
ITM <- n + (1 - pa1)*(N - n)
ITM## [1] 2953.843Note que ainda mantém o \(n < ITM < N\), isso é bom, porém fica bem próximo ao lote total, ou seja, quase todo o lote foi revistado.
d)
Apresente as curvas para QSM e ITM versus p. Identifique o máximo de QSM.
#Repetindo os valores definidos anteriormente.
n<-150
Ac<-2
N<-3000
p<-seq(0,0.20,0.001)#Gráfico QMS
pa1<-phyper(Ac,p*N, N-p*N, n)
QMS<-((pa1*(N-n))/(N))*p
plot(p,QMS, type="l",ylim=c(0,0.010),xlab = "Qualidade do lote que entra")
max(QMS)## [1] 0.008653284abline(h=max(QMS))
No gráfico do ITM, podemos perceber que cada vez que aumenta a fração de não conformes no lote, maior o número de itens a ser especionado.
#Gráfico ITM
ITM <- n + (1 - pa1)*(N - n)
plot(p,ITM, type="l",xlab = "Qualidade do lote que entra")
6)
Considerando um plano de amostragem dupla com \(n1 = 50\), \(Ac1 = 2\), \(n2 = 100\) e \(Ac2 = 6\). A fração de itens defeituosos produzida é \(p = 0.05\).
Quantidades necessárias para o cálculo.
n1<-50
n2<-100
c1<-2
c2<-6
p<-0.05a)
- Qual a probabilidade de aceitação na primeira amostra
P1A = pbinom(c1,n1,p)
P1A## [1] 0.5405331b)
- Qual a probabilidade de aceitação na segunda amostra?
P2A=0
for(d1 in (c1+1):c2)
{
for(d2 in 0:(c2-d1))
{
aux = dbinom(d1,n1,p)*dbinom(d2,n2,p)
P2A = P2A + aux
}
}
P2A## [1] 0.07536839c)
Apresente a curva característica de operação para o plano apresentado. Para qual valor de p teríamos probabilidade de aceitação próxima a \(0.9\)?
p<-seq(0,0.40,0.01)
P2A=0
P1A = pbinom(c1,n1,p)
for(d1 in (c1+1):c2)
{
for(d2 in 0:(c2-d1))
{
aux = dbinom(d1,n1,p)*dbinom(d2,n2,p)
P2A = P2A + aux
}
}
#Probabilidade de aceitação
PA = P1A + P2A
PA## [1] 1.000000e+00 9.996175e-01 9.846868e-01 9.146147e-01 7.799286e-01
## [6] 6.159015e-01 4.607475e-01 3.332377e-01 2.359896e-01 1.645840e-01
## [11] 1.132300e-01 7.684545e-02 5.143227e-02 3.394557e-02 2.209742e-02
## [16] 1.419266e-02 8.997653e-03 5.632672e-03 3.483194e-03 2.128393e-03
## [21] 1.285419e-03 7.674297e-04 4.529963e-04 2.643953e-04 1.525944e-04
## [26] 8.708769e-05 4.914744e-05 2.742525e-05 1.513120e-05 8.253258e-06
## [31] 4.449916e-06 2.371307e-06 1.248700e-06 6.496480e-07 3.338517e-07
## [36] 1.694262e-07 8.488790e-08 4.197840e-08 2.048273e-08 9.857995e-09
## [41] 4.678152e-09plot(p,PA, type="l", ylim=c(0,1.05),ylab = "Probabilidade Aceitação", xlab="Fração de defeituosos no lote")
abline(h=0.9, col="red")
abline(v=0.032, col="red")
O valor de \(p\) tal que a probabilidade de aceitação é de \(0.9\) é de 0.032 aproximadamente
d) Qual o tamanho amostral médio (TAM) do plano apresentado?
Informações anteriores
n1<-50
n2<-100
c1<-2
c2<-6
p<-0.05Probabilidade de aceitação da primeira amostra.
P1A = pbinom(c1,n1,p)
P1A## [1] 0.5405331#Probabilidade de rejeição da primeira amostra.
P1R = 1 - pbinom(c2,n1,p)
P1R## [1] 0.01178645Tamanho amostral médio. O TAM equivale ao tamanho de \(n1\) somada ao tamanho de \(n2\), ponderada pela chance de ocorrência de uma \(2º\) amostra.
TAM <- n1 + n2*(1-P1A - P1R)
TAM## [1] 94.76804e)
Obtenha a curva TAM versus \(p\)
Informações anteriores
p<-seq(0,0.40,0.001)
n1<-50
n2<-100
c1<-2
c2<-6Probabilidade de aceitação da primeira amostra.
P1A = pbinom(c1,n1,p)Probabilidade de rejeição da primeira amostra.
P1R = 1 - pbinom(c2,n1,p)#Tamnho amostral médio
TAM <- n1 + n2*(1-P1A - P1R)
plot(p,TAM,type="l",ylab="Tamanho amostral médio",xlab="Fração de defeituosos no lote")