使用R做ARL的蒙地卡羅模擬

這是上學期上品質管理的時侯 老師出的作業之一

用任何一種程式語言模擬製程偏移之後 ARL的表現

其實我本來也不會寫 後來在網路上查到有人用C++寫好了 但我不好意思照抄 所以就把C++轉成R

因為品管的老師很喜歡出這個題目 我就放上來解救大家

但是每一個老師的假設都不一樣，所以如果有人要抄的話還是要注意一下

這是假設製程在第100點以後開始偏移，所以前100點是服從N~(0,1)，100點以後就看老師要求的偏移量是什麼，去改k的值

a<-0              #在陣列中'列'的位移

n<-20           #樣本數

m1<-0          #原始mean

sigma<-1     #原始sigma

L<-3              #管制界線的寛度

K<-1.5          #偏移程度

m2<-m1+K*sigma                                                   #計算偏移後的mean

UCL<-m1+L*sigma/sqrt(n)                                     #計算UCL

LCL<-m1-L*sigma/sqrt(n)                                       #計算LCL

D<-matrix(0, nrow=200000, ncol=3, byrow=F)   #宣告一個200000*3的矩陣以紀錄樣本點與xar

RL<-matrix(0, nrow=1000, ncol=2, byrow=F)      #宣告一個1000*2的矩陣以紀錄發生警示的樣本點與xar

for(i in seq(1,1000,by=1))     #開始模擬

{ j<-0                                         #j為記綠RL

 repeat                                     #產生前100次In control的normal(0,1)抽樣

 {

   a<-a+1

   j<-j+1           #記綠RL

   xbar1<-mean(rnorm(n,mean=0,sd=1)) #產生100次normal(0,1)的抽樣

   D[a,1]=i        #寫入記錄表的第1列，用於模擬次數

   D[a,2]=j        #寫入記錄表的第2列，用於記錄抽樣次數

   D[a,3]=xbar1 #寫入記錄表的第3列，用於記錄每次抽樣的xbar

   if (j>99) break

 }

 repeat            #假設製程於第100次抽樣後開始失控，記綠RL與xbar

 {

   j<-j+1

   a<-a+1

   xbar2<-mean(rnorm(n,mean=m2,sd=sigma)) #抽取n個normal(m2,sd)樣本，並計算xbar

   D[a,1]=i #寫入記錄表的第1列，用於模擬次數

   D[a,2]=j #寫入記錄表的第2列，用於記錄抽樣次數

   D[a,3]=xbar2 #寫入記錄表的第3列，用於記錄每次抽樣的xbar

   if (xbar2>UCL | xbar2<LCL)  #判斷xbar是否超出管制線

   {

     RL[i,1]=j                 #如果超出管制線則記錄其RL(run length)

     RL[i,2]=xbar2        # mark the xbar

     break

   }

 }

}

print(n) #顯示n

print(K) #顯示k

print(ARL<-(mean(RL[,1])-100))   #計算ARL

print(ARLSD<-sd(RL[,1]))               #計算ARL的標準差