如何用r语言采用蒙特卡罗模拟5000次计算出EWMA控制图的平均运行链长ARL-经管之家官网！

附上我参考的文献及我自己写的代码，根据参考文献写的，但是运行结果反过来了 ，求大神援助帮我改对代码就可以，事关毕业论文，谢谢啦，重赏！！！1
#指数分布，参数为R
#正态分布0.9975对应 的分位点为2.81,即ARL0取值为200时对应的L值为2.81
#lambda=1,样本数为200，每个样本容量为7，即可控下的数据位200，第201个数据开始发生偏移。
lambda=0.25
L=2.81
shift=c(1,1.1,1.2,1.3,1.4,1.5,1.6,1.8,2,2.5,3,3.5,4)#标准差偏移度i
l=length(shift) #13个不同的偏移度
n=100#j
newsigma=matrix(,n,l)
rl=matrix(1300,n,l)
for(i in 1:l){
for(j in 1:n){
m=matrix(rexp(1400,1),200,7)#生成200个样本容量为7，参数lambda为1的指数分布
#参数R
R=c()
Rsigma=c()
Rz=c()
sigma=sd(m)
for(k in 1:nrow(m)){
R[k]=max(m[k,])-min(m[k,])#求每个样本极差
Rz[k]=R[k]/sigma #求每个样本参数Rz
}
Rbar=mean(R)
Rt2=mean(Rz)
Rt3=sd(Rz)
LCL=Rbar-(L*Rt3*Rbar*(lambda/(2-lambda))^0.5)/Rt2#控制下限
UCL=Rbar+(L*Rt3*Rbar*(lambda/(2-lambda))^0.5)/Rt2#控制上限
sdw=(Rt3*Rbar*(lambda/(2-lambda))^0.5)/Rt2 #计算过程方差
newsigma[j,i]=shift*sdw #发生偏移后的标准差
ewma=c()
num=1
ewma[1]=Rbar
for(k in 2:nrow(m)){
ewma[k]=lambda*newsigma[j,i]+(1-lambda)*ewma[k-1] #发生偏移后的统计量
if (ewma[k]>LCL&&ewma[k]<UCL){ #判断是否在区间范围内
num=num+1 #计算当出现第一个不在控制区间时已经落在控制上下限的总数
}
else{break}
}
rl[j,i]=num #求每一次的运行链长
}
}
rl
arl=c()
for(i in 1:l){
arl=mean(rl[,i])
}
arl
arl