monrepertoire <- "Mon repertoire de travail"
setwd(monrepertoire)
diams <- read.csv("TD_CCmes_ex1.csv")
str(diams)
## 'data.frame': 100 obs. of 2 variables:
## $ Echantillon: int 1 1 1 1 1 2 2 2 2 2 ...
## $ Mesure : num 61.9 62.1 62.4 62.1 62.2 61.8 62.1 61.9 62 62 ...
head(diams)
## Echantillon Mesure
## 1 1 61.9
## 2 1 62.1
## 3 1 62.4
## 4 1 62.1
## 5 1 62.2
## 6 2 61.8
tail(diams)
## Echantillon Mesure
## 95 19 61.8
## 96 20 61.9
## 97 20 61.8
## 98 20 62.5
## 99 20 62.3
## 100 20 62.1
#Nombre d'échantillons
(nech=length(unique(diams$Echantillon)))
## [1] 20
#Nombre de mesures total
(nmestot=nrow(diams))
## [1] 100
#Nombre de mesures par echantillon
(nmes=nrow(diams)/length(unique(diams$Echantillon)))
## [1] 5
#Calcul des ordonnées des points de la carte xbar (moyennes des échantillons)
(xbari=tapply(diams$Mesure,diams$Echantillon,mean))
## 1 2 3 4 5 6 7 8 9 10 11 12
## 62.14 61.96 62.12 62.20 62.06 62.16 61.90 62.28 62.18 62.30 61.98 62.22
## 13 14 15 16 17 18 19 20
## 62.20 62.04 62.12 62.06 62.04 62.18 62.06 62.12
unlist(lapply(split(diams$Mesure,diams$Echantillon),mean))
## 1 2 3 4 5 6 7 8 9 10 11 12
## 62.14 61.96 62.12 62.20 62.06 62.16 61.90 62.28 62.18 62.30 61.98 62.22
## 13 14 15 16 17 18 19 20
## 62.20 62.04 62.12 62.06 62.04 62.18 62.06 62.12
#Construction de la carte : ajout des échantillons
plot(xbari)
#Calcul de m0
(m0=mean(tapply(diams$Mesure,diams$Echantillon,mean)))
## [1] 62.116
#Ajout de la ligne centrale à la main
abline(h=m0,lwd=2)
#Calcul des écart-types ("bruts") de chaque échantillon
tapply(diams$Mesure,diams$Echantillon,sd)
## 1 2 3 4 5 6
## 0.18165902 0.11401754 0.28635642 0.21213203 0.15165751 0.21908902
## 7 8 9 10 11 12
## 0.20000000 0.20493902 0.14832397 0.38078866 0.28635642 0.22803509
## 13 14 15 16 17 18
## 0.07071068 0.11401754 0.33466401 0.19493589 0.31304952 0.08366600
## 19 20
## 0.23021729 0.28635642
#Calcul de s0
(s0=sqrt(mean(tapply(diams$Mesure,diams$Echantillon,var))))
## [1] 0.2274863
#Ajout des limites de controle de la carte xbar à la main
abline(h=m0-3*s0/sqrt(nmes),col="red",lty=2,lwd=2)
abline(h=m0+3*s0/sqrt(nmes),col="red",lty=2,lwd=2)
#Les limites de controle sont en dehors de la zone tracée
(ylimcarte=range(xbari,m0-3*s0/sqrt(nmes),m0+3*s0/sqrt(nmes)))
## [1] 61.8108 62.4212
plot(xbari,ylim=ylimcarte)
abline(h=m0,lwd=2)
abline(h=m0-3*s0/sqrt(nmes),col="red",lty=2,lwd=2)
abline(h=m0+3*s0/sqrt(nmes),col="red",lty=2,lwd=2)
#normalité ?
qqnorm(diams$Mesure)
qqline(diams$Mesure)
#Calcul de la ligne centrale de la carte S
(s0cor=gamma(nmes/2)/gamma((nmes-1)/2)*sqrt(2/(nmes-1))*s0)
## [1] 0.2138338
library(qcc)
## Package 'qcc', version 2.6
## Type 'citation("qcc")' for citing this R package in publications.
?qcc
diams.row <- matrix(diams$Mesure,nrow=20,ncol=5,byrow=TRUE)
#Ou bien la fonction suivante fournie par le package qcc
diams.row <- qcc.groups(diams$Mesure,diams$Echantillon)
qcc(diams.row,type="xbar",std.dev="UWAVE-R")
## List of 11
## $ call : language qcc(data = diams.row, type = "xbar", std.dev = "UWAVE-R")
## $ type : chr "xbar"
## $ data.name : chr "diams.row"
## $ data : num [1:20, 1:5] 61.9 61.8 62.1 62.4 62.3 61.9 61.8 62.5 62.4 62.2 ...
## ..- attr(*, "dimnames")=List of 2
## $ statistics: Named num [1:20] 62.1 62 62.1 62.2 62.1 ...
## ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
## $ sizes : Named int [1:20] 5 5 5 5 5 5 5 5 5 5 ...
## ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
## $ center : num 62.1
## $ std.dev : num 0.23
## $ nsigmas : num 3
## $ limits : num [1, 1:2] 61.8 62.4
## ..- attr(*, "dimnames")=List of 2
## $ violations:List of 2
## - attr(*, "class")= chr "qcc"
#Différentes estimations de l'écart-type
qcc(diams.row,type="xbar",std.dev="UWAVE-R")$std.dev
## [1] 0.2300086
qcc(diams.row,type="xbar",std.dev="UWAVE-SD")$std.dev
## [1] 0.2255871
qcc(diams.row,type="xbar",std.dev="MVLUE-R")$std.dev
## [1] 0.2300086
qcc(diams.row,type="xbar",std.dev="MVLUE-SD")$std.dev
## [1] 0.2255871
qcc(diams.row,type="xbar",std.dev="RMSDF")$std.dev
## [1] 0.2281983
qcc(diams.row,type="S")
## List of 11
## $ call : language qcc(data = diams.row, type = "S")
## $ type : chr "S"
## $ data.name : chr "diams.row"
## $ data : num [1:20, 1:5] 61.9 61.8 62.1 62.4 62.3 61.9 61.8 62.5 62.4 62.2 ...
## ..- attr(*, "dimnames")=List of 2
## $ statistics: Named num [1:20] 0.182 0.114 0.286 0.212 0.152 ...
## ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
## $ sizes : Named int [1:20] 5 5 5 5 5 5 5 5 5 5 ...
## ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
## $ center : num 0.212
## $ std.dev : num 0.226
## $ nsigmas : num 3
## $ limits : num [1, 1:2] 0 0.443
## ..- attr(*, "dimnames")=List of 2
## $ violations:List of 2
## - attr(*, "class")= chr "qcc"
#Exercice 3
viscosite=c(2838,2785,3058,3064,2996,2782,2878,2920,3050,2870,3174,3102,2762,2975,2719,2861,2797,3078,2974,2805,3163,3199,3054,3147,3156)
qcc(data = viscosite,type="xbar.one")
## List of 11
## $ call : language qcc(data = viscosite, type = "xbar.one")
## $ type : chr "xbar.one"
## $ data.name : chr "viscosite"
## $ data : num [1:25, 1] 2838 2785 3058 3064 2996 ...
## ..- attr(*, "dimnames")=List of 2
## $ statistics: Named num [1:25] 2838 2785 3058 3064 2996 ...
## ..- attr(*, "names")= chr [1:25] "1" "2" "3" "4" ...
## $ sizes : int [1:25] 1 1 1 1 1 1 1 1 1 1 ...
## $ center : num 2968
## $ std.dev : num 135
## $ nsigmas : num 3
## $ limits : num [1, 1:2] 2564 3373
## ..- attr(*, "dimnames")=List of 2
## $ violations:List of 2
## - attr(*, "class")= chr "qcc"
(mRi=abs(diff(viscosite)))
## [1] 53 273 6 68 214 96 42 130 180 304 72 340 213 256 142 64 281
## [18] 104 169 358 36 145 93 9
(m0=mean(viscosite))
## [1] 2968.28
(s0=2.66*mean(mRi)/3)
## [1] 134.7733
(LCL=m0-2.66*mean(mRi))
## [1] 2563.96
(UCL=m0+2.66*mean(mRi))
## [1] 3372.6
#Exercice 4
#cible 80 -0.65 +.8
Ti=80-0.65
Ts=80+.8
m0=80
s0=.15
alpha=.27/100
n=9
#1)
m0+c(-1,1)*qnorm(1-alpha/2)*s0/sqrt(n)
## [1] 79.85 80.15
LCI_Xbar=m0-qnorm(1-alpha/2)*s0/sqrt(n)
LCS_Xbar=m0+qnorm(1-alpha/2)*s0/sqrt(n)
#2)
LCI_S=s0/sqrt(n-1)*sqrt(qchisq(alpha/2,n-1))
LCS_S=s0/sqrt(n-1)*sqrt(qchisq(1-alpha/2,n-1))
#3)
(m1=Ts-3*s0)
## [1] 80.35
(m1prime=Ti+3*s0)
## [1] 79.8
(delta1=min((m1-m0)/s0,(m0-m1prime)/s0))
## [1] 1.333333
#4)
#Prob_delta1(LCI<=Xbar<=LCS)
#Calcul approché beta=Prob_delta1(LCI<=Xbar)
1-pnorm(LCI_Xbar,m1prime,s0/sqrt(n))
## [1] 0.1586497
#Calcul exact beta=Prob_delta1(LCI<=Xbar)-Prob_delta1(LCS<=Xbar)
1-pnorm(LCI_Xbar,m1prime,s0/sqrt(n))-(1-pnorm(LCS_Xbar,m1prime,s0/sqrt(n)))
## [1] 0.1586497
#Erreur commise par l'approximation
#Toujours majoration du risque
(1-pnorm(LCI_Xbar,m1prime,s0/sqrt(n))-(1-pnorm(LCI_Xbar,m1prime,s0/sqrt(n))-(1-pnorm(LCS_Xbar,m1prime,s0/sqrt(n)))))/min(1-pnorm(LCI_Xbar,m1prime,s0/sqrt(n)),1-pnorm(LCI_Xbar,m1prime,s0/sqrt(n))-(1-pnorm(LCS_Xbar,m1prime,s0/sqrt(n))))*100
## [1] 8.06794e-10
#u_beta=u1prime-(m0-m1prime)/s0*sqrt(n)
#Ici carte standard avec alpha=0,27%
#symétrique -> u1prime=3
u1prime=3
beta=0.05
u_beta=qnorm(beta)
((u_beta-u1prime)*s0/(m0-m1prime))^2
## [1] 12.13575
#Si risque max 5% alors n=13
LCI_Xbar_13=m0-qnorm(1-alpha/2)*s0/sqrt(13)
1-pnorm(LCI_Xbar_13,m1prime,s0/sqrt(13))
## [1] 0.03534804
#Si risque de l'ordre de 5% alors n=12
LCI_Xbar_12=m0-qnorm(1-alpha/2)*s0/sqrt(12)
1-pnorm(LCI_Xbar_12,m1prime,s0/sqrt(12))
## [1] 0.05274244
#6)
#u_beta=u1prime-(m0-m1prime)/s0*sqrt(n)
u1prime=u_beta+(m0-m1prime)/s0*sqrt(n)
#u1prime=qnorm(1-alpha/2)
u_beta=qnorm(.05)
(alpha=2-2*pnorm(u_beta+(m0-m1prime)/s0*sqrt(n)))
## [1] 0.01851541
#alpha passe de 0,27% à alpha=1,8% lorsque beta est fixé à 5% avec un effectif constant n=9.