你好！谢谢上次你对我的帮助，我想请问一下，那个R是说的R语言吗？可以应用到MATLAB中吗？我想再向你请教一个问题，这个ARL的值是和控制限L和参数r有关的，如果在既定ARL情况下，如何用搜索的方法得到控制限L，同学跟我说用搜索的方法，比如二分法或黄金分割法来得到，我没怎么明白，您这几次都热心帮助我，所以想问问您知道吗。或者说ARL、控制限L和参数r如何进行优化使得控制图性能最优，您了解这方面吗？想再麻烦您再给解答一下，谢谢啦！！！

给定ARL,列出objective function
要找出最佳参数(r,L),方法很多
你可参考底下网站:
http://numericalmethods.eng.usf.edu/resource.html
2楼我给的是用Genetic Algorithm
也就是基因算法.
如果你需要的话.
我可以给你OptMEWMA_FINAL.exe

R package "spc"
xewma.crit() 计算critical values of EWMA charts
是采secant rule

其实安装R很简单,只需十分钟.
使用起来也很灵活.
http://cran.r-project.org/
Windows\base\Download R 2.11.1 for Windows

谢谢指导，哈哈~~~~~

你好！非常感谢您去年对我的帮助，我想再问一下，还是那个问题，R是一个k×k阶的矩阵，其中每个元素Pij可以表示为：
Pij=Ф[((j*d)-(1-r)*(i-0.5)*d)/(r)]- Ф[((j-1)*d-(1-r)*(i-0.5)*d)/(r)] ；Ф为标准正态概率分布函数，其中，i,j=1,2,…,k;r=0.1,0.2,….0.9;
ARL=P(I-R)-11
其中，P=[0,…,1,…,0]1×k,P(1,k/2)=1;I为k阶单位阵，1为所有元素均为1的k×1的列向量.
在求控制限的时候，我对ARL的值进行矩阵计算，得到下面这个式子：
Rk×k=I-1/ARL*Ik×1*P1×k
等式右边是已知的，根据这个式子，可以建立目标函数求解L和r吗？
我不知道是不是这个问题很幼稚，可能对matlaB和算法了解太少了，我想问一下我这个思路对不对，能不能再麻烦您给指点一下？

能麻烦您再跟我说下那个已知ARL的值，怎么列出目标函数？能不能麻烦您再指点一下？谢谢啦！！！

麻烦能不能再跟我说一下啊那个已知ARL的值，如何列出目标函数，我算了很久，没有找到思路，能不能麻烦您再给指点一下？谢谢啦！！！

optimpara.m

%x(1)=r smoothing parameter
%x(2)=L control limit
%in matlab command window

ARL0=500;  %in-control ARL
k=100;
x0 = [0.13;2.5];  %initial value of r & L
lb=[0.1;2.3];      %lower bound of r & L
ub=[0.2;3.3];     %upper bound of r & L
[x,fval] = fmincon(@(x)optimpara(x,ARL0,k),x0,[],[],[],[],lb,ub)
%x =
%     0.1093  (r)
%     2.8374  (L)

谢谢您的帮助，真是十分感谢！！！！

您好！可以帮我看看这段程序吗？好像是矩阵a\b维度不同问题，我修改不好，麻烦您指导一下，谢谢了！！！
echo off all
clear
clc
m=3;mu=0;v=1;s=0;d=3;fi=0.6;
rfi(1)=0;rfi(2)=0.3;rfi(3)=-0.3;rfi(4)=0.6;rfi(5)=-0.6;rfi(6)=0.9;rfi(7)=-0.9;
deta=0:0.25:d;h=0.05;n=2*fix(m/h);nn=n+1;
one=eye(nn,nn); b=ones(nn,1)*-1;
ai=zeros(nn,1);
bi=ai;n1=ai;n2=ai;
y=zeros(7,length(deta));
for i=1:nn
    bi(i)=mu+m;
    ai(i)=mu-m;
    n1(i)=fix((ai(i)+m-mu)/h)+1;
    n2(i)=fix((bi(i)+m-mu)/h)+1;
end
for r=1:7
for k=1:length(deta)
    a=zeros(nn,nn);
    for i=1:nn
        at=zeros(1,nn);
        for j=1:nn
            if j>n1(i) && j<n2(i)
                xx=ai(i)+h*(j-n1(i))-rfi(r)*s*(-m+h*(i-1))-(1-rfi(r))*(mu+deta(k)*v);
                at(1,j)=h*normpdf(xx,0,v);
            elseif j==n1(i)
                xx=ai(i)-rfi(r)*s*(-m+h*(i-1))-(1-rfi(r))*(mu+deta(k)*v);
                at(1,j)=(h/2)*normpdf(xx,0,v);
            elseif j==n2(i)
                xx=bi(i)-rfi(r)*s*(-m+h*(i-1))-(1-rfi(r))*(mu+deta(k)*v);
                at(1,j)=(h/2)*normpdf(xx,0,v);
            end
        end
            a(i,:)=at(1,:);
    end
       a=a-one;
       if det(a)==0
           disp "ZERO  DETERMINANT FOR LIEAR SYSTEM"
       else
           x=a/b;
           arl(k)=(h/2)*(normdf(mu-m-(1-rfi(r))*(mu+deta(k)*v),0,1)*x(1)+normpdf(mu+m-(1-rfi(r))*(mu+deta(k)*v),0,1)*x(nn));
           arg=noempdf(h*(1:1:n-1)+mu-m-(1-rfi(r))*(mu+deta(k)*v),0,1);
           x=x(2:n);
           arl(k)=arl(k)+arg*h*x+1;
       end
end
y(r,:)=arl;
end

4# epoh

a*x=b;
       x=inv(a)*b;
or  x=a\b;     NOT x=a/b;

%%%%%%
>>arl
369.640058965427        165.297773402873        49.4015923573809        17.3397569610425       
7.36830241838697        3.75899855126688        2.27069067418857        1.59392055997436       
1.26888044175189        1.11259659132598        1.04176499631875        1.01322889570437       
1.00348493714105

您好！我运行R，结果说没有这个函数，像下面说的：

> xewma.ad(l=0.13,c=2.88,mu1=0,mu0=0,zr=0,sided="two",limits="fix",r=40)
错误: 没有"xewma.ad"这个函数

请问这是什么原因？  谢谢！！！
4# epoh

pls install package "spc"
in R command window
install.packages("spc")

then
library(spc)

谢谢高人指导！！！！ 18# epoh

您好！再请教一个问题，还是这个积分函数，开始是基于AR（1）过程求休哈特图的ARL，如果求解EWMA图的ARL，怎样修改积分限，谢谢！
echo off all
clear
clc
m=3;mu=0;v=1;s=0;d=3;fi=0.6;
rfi(1)=0;rfi(2)=0.3;rfi(3)=-0.3;rfi(4)=0.6;rfi(5)=-0.6;rfi(6)=0.9;rfi(7)=-0.9;
deta=0:0.25:d;h=0.05;n=2*fix(m/h);nn=n+1;
one=eye(nn,nn); b=ones(nn,1)*-1;
ai=zeros(nn,1);
bi=ai;n1=ai;n2=ai;
y=zeros(7,length(deta));
for i=1:nn
    bi(i)=mu+m;
    ai(i)=mu-m;
    n1(i)=fix((ai(i)+m-mu)/h)+1;
    n2(i)=fix((bi(i)+m-mu)/h)+1;
end
for r=1:7
for k=1:length(deta)
    a=zeros(nn,nn);
    for i=1:nn
        at=zeros(1,nn);
        for j=1:nn
            if j>n1(i) && j<n2(i)
                xx=ai(i)+h*(j-n1(i))-rfi(r)*s*(-m+h*(i-1))-(1-rfi(r))*(mu+deta(k)*v);
                at(1,j)=h*normpdf(xx,0,v);
            elseif j==n1(i)
                xx=ai(i)-rfi(r)*s*(-m+h*(i-1))-(1-rfi(r))*(mu+deta(k)*v);
                at(1,j)=(h/2)*normpdf(xx,0,v);
            elseif j==n2(i)
                xx=bi(i)-rfi(r)*s*(-m+h*(i-1))-(1-rfi(r))*(mu+deta(k)*v);
                at(1,j)=(h/2)*normpdf(xx,0,v);
            end
        end
            a(i,:)=at(1,:);
    end
       a=a-one;
       if det(a)==0
           disp "ZERO  DETERMINANT FOR LIEAR SYSTEM"
       else
           x=a\b;
           arl(k)=(h/2)*(normpdf(mu-m-(1-rfi(r))*(mu+deta(k)*v),0,1)*x(1)+normpdf(mu+m-(1-rfi(r))*(mu+deta(k)*v),0,1)*x(nn));
           arg=normpdf(h*(1:1:n-1)+mu-m-(1-rfi(r))*(mu+deta(k)*v),0,1);
           x=x(2:n);
           arl(k)=arl(k)+arg*h*x+1;
       end
end
y(r,:)=arl;
end

18# epoh

4楼code是求解ARL.

13楼optimpara.m
是求解smoothing parameter
           &control limit

两个你可试过?

您好！！！非常感谢你的热心帮助。
4楼的编程我试过可以得到ARL值，若过程的均值偏移为delta=0.5,方差为1，自相关系数fi为0，上下限从（0.759352，-0.759352）变为（1.259352，-0.259352）可得到的ARL值都是503.44，请问要考虑偏移和自相关如何修改程序？



21# epoh

您好，请指点一下，好吗？这是基于AR(1)模型求解残差EWMA图的ARL,但运行结果是0，与实际结果370不符，请高人指点，万分感激！

function y=arl_ewma(lmda,K)
rn=randn(1,10000+1);
beta=0;  %AR(1)自回归系数
delta=0; %均值偏移系数
mu=(1-beta)*delta*sqrt(1/(1-beta^2))*var(rn);%残差均值偏移量
UCL=mu+K*sqrt(lmda/(2-lmda))*var(rn);%上限
LCL=mu-K*sqrt(lmda/(2-lmda))*var(rn);%下限
sum=0;
rep=10000;  Z=1;
for i=1:rep
    r1=0;n=1;
    while n<10000   
       if Z>UCL && Z<LCL  
            Z=lmda*rn(1,i)+(1-lmda)*Z;
        r1=r1+1;
        break
        end
       n=n+1;
    end
    sum=sum+r1;
end
y=sum/rep;
21# epoh

in R
xewma.arl(l=0.13,c=2.88,mu=0,zr=0,hs=0,sided="two",limits="fix",q=1,r=40)  
#504.8908  
                           
xewma.arl(l=0.13,c=2.88,mu=0.5,zr=0,hs=0,sided="two",limits="fix",q=1,r=40)  
# 34.08857

%%in Matlab
T1=normcdf(temp,0,1);  => T1=normcdf(temp,0.5,1);
T2=normcdf(temp,0,1);  => T2=normcdf(temp,0.5,1);

ans:34.1734

看不太懂,请解释一下
你设置beta=0;  delta=0;
当然 mu=0

你设置 Z=1;
主要是让 if Z>UCL && Z<LCL  条件达不到??

r1一直是0,r1代表啥么?

lmda=?
K=?

谢谢您的指导！deltao=0;beta=0;是在过程没有自相关的影响，没有均值偏移的情况下的ARL值，对于Z为ewma残差的统计量初值应该设为0，上面错了。 if Z>UCL && Z<LCL  是EWMA统计量满足上下控制限，若超出则跳出循环，r1是记录循环次数，即运行长度RL。lmda是EWMA的平滑系数，值为0.2,；K为参数值为3.



25# epoh

lmda=0.2;
K=3;
beta=0;  %AR(1)自回归系数
delta=0; %均值偏移系数
n=1;
sum=0;
rep=10001;
while(n<rep)
rn=randn(1,rep);
mu=(1-beta)*delta*sqrt(1/(1-beta^2))*var(rn);%残差均值偏移量
UCL=mu+K*sqrt(lmda/(2-lmda))*var(rn);%上限   
LCL=mu-K*sqrt(lmda/(2-lmda))*var(rn);%下限  

Z(1)=0;
for i=2:rep
   Z(i)=lmda*rn(1,i)+(1-lmda)*Z(i-1) ;
end
Z=Z(1,2:rep);
idx=find(Z>UCL | Z<LCL,1);
if(isempty(idx)==1) idx=0;
end
sum=sum+idx;
n=n+1;
end      
y=sum/(rep-1)     %558.6505

%%%%%%%%%%%%%%%%%%%%%%%%%%%
in R
Computation of the (zero-state) ARL
l=0.2
c=3
mu=0
xewma.arl(l,c,mu,zr=0,hs=0,sided="two")
#     arl
#559.8741

非常感谢您的帮助！！！
27# epoh

好东西啊 顶一个

epoh老师好，我想请问您还有OptMEWMA_FINAL.exe这个吗？您能发我一份吗？邮箱 s18237106122@gmail.com   万分感谢！我是质量管理专业的学生，现做本科毕业论文，题目是多元统计过程控制。我想请教个问题：我求得的mcusum控制图在受控状态下的状态转移概率矩阵，每行所有元素之和不等于1，这是不是定错了？