先确认几件事
1. elseif k~=0&&l~=0
          temp=min(b2,max(a1,b1));
          T1=normcdf(temp,0,1);
          temp=max(b1,max(a2,b2));  OR temp=max(b1,min(a2,b2));

2.a1=(k-i-0.5)*w*(log(1+250*(m+l*q)^2+2.332.....
  a2=(k-i+0.5)*w*(log(1+250*(m+l*q)^2+2.332.....
  250如何定的
  这会影响结果
3.我先把probability matrix
  放进,你可试试

m1=17;
m2=25;

R=zeros(m1*m2,m1*m2);
for i=0:1:m1-1
  for j=0:1:m2-1
    for k=0:1:m1-1
      for l=0:1:m2-1
      a1=(k-i-0.5)*w*(log(1+250*(m+l*q)^2+2.332*(m+l*q)/2)/(m+l*q)-1.166)+(m+l*q)/2;
      a2=(k-i+0.5)*w*(log(1+250*(m+l*q)^2+2.332*(m+l*q)/2)/(m+l*q)-1.166)+(m+l*q)/2;

      b1=(m+(l-0.5)*q-(1-p)*(m+j*q))/p;
      b2=(m+(l+0.5)*q-(1-p)*(m+j*q))/p;
       if k==0&&l==0
          temp=min(a2,b2);
          T=normcdf(temp,0,1);
          R(i*m2+(j+1),k*m2+(l+1))=
       elseif k==0&&l~=0
          temp=b1;
          T1=normcdf(temp,0,1);
          temp=max(b1,min(a2,b2));
          T2=normcdf(temp,0,1);
          R(i*m2+(j+1),k*m2+(l+1))=
       elseif k~=0&&l==0
           temp=a1;
           T1=normcdf(temp,0,1);
           temp=max(a1,min(a2,b2));
           T2=normcdf(temp,0,1);
           R(i*m2+(j+1),k*m2+(l+1))=
        elseif k~=0&&l~=0
           temp=min(b2,max(a1,b1));
           T1=normcdf(temp,0,1);
           temp=max(b1,min(a2,b2));
           T2=normcdf(temp,0,1);
           R(i*m2+(j+1),k*m2+(l+1))=
         end;  %end if
       end;
     end;
  end;
end;

1. elseif k~=0&&l~=0
          temp=min(b2,max(a1,b1));
          T1=normcdf(temp,0,1);
          temp=max(b1,max(a2,b2));  OR temp=max(b1,min(a2,b2));是后面这个，就是您写在程序里面的那个，我当时大意了，嘿嘿！！

2.a1=(k-i-0.5)*w*(log(1+250*(m+l*q)^2+2.332.....
  a2=(k-i+0.5)*w*(log(1+250*(m+l*q)^2+2.332.....
  这里面涉及到adaptive CUSUM控制图改进过程中的一个公式，就是h(k)=log(1+2*k^2*ARL0+2.332*k)/(2*k)-1.166,程序里面的k=(m+l*q)/2,我令ARL0=400，我把这个2*k^2*ARL0公式中的平方提出来，稍微化简了一下，就出现了250，其实可以表示为1000*((m+l*q)/2)^2.下面我又改了一下，

我放进去的程序是
function y=ARL(w,m,q,p)
m1=17;m2=25;

R=zeros(m1*m2,m1*m2);
for i=0:1:m1-1
for j=0:1:m2-1
     for k=0:1:m1-1
     for l=0:1:m2-1

         h=(m+l*q)/2;
         aa=log10(1+2*h*h*400+2.332*h)/(2*h)-1.166;
     a1=(k-i-0.5)*w*aa+h;
     a2=(k-i+0.5)*w*aa+h;

     b1=(m+(l-0.5)*q-(1-p)*(m+j*q))/p;
     b2=(m+(l+0.5)*q-(1-p)*(m+j*q))/p;
     if k==0&&l==0
         temp=min(a2,b2);
         T=normcdf(temp,0,1);
         R(i*m2+(j+1),k*m2+(l+1))=T;
     end

      if k==0&&l~=0
          temp=b1;
          T1=normcdf(temp,0,1);
          temp=max(b1,min(a2,b2));
          T2=normcdf(temp,0,1);
          R(i*m2+(j+1),k*m2+(l+1))=T2-T1;
         end

          if k~=0&&l==0
               temp=a1;
               T1=normcdf(temp,0,1);
               temp=max(a1,min(a2,b2));
               T2=normcdf(temp,0,1);
                R(i*m2+(j+1),k*m2+(l+1))=T2-T1;
          end

              if k~=0&&l~=0
                temp=min(b2,max(a1,b1));
                 T1=normcdf(temp,0,1);
                temp=max(b1,min(a2,b2));
                T2=normcdf(temp,0,1);
                R(i*m2+(j+1),k*m2+(l+1))=T2-T1;            
                 end;
             end;
         end;
     end;
     end;

PM=0*ones(1,m1*m2);
PM(212)=1;
y=PM*inv(eye(m1*m2,m1*m2)-R)*ones(m1*m2,1);


程序里面这个PM初始状态矩阵我不知道选择的对不对，我设置的参数为y=ARL(0.06455,1,0.12245,0.1)。
也就是控制限c为1.065，δmin最小值设为1，Qt的控制限L设为4，l设为0.1,
然后计算出来的 w=2*c/(2*m1-1=)0.06455;
                          δmin=1;
           q=2*(L-δmin)/(2m2-1)=0.122245

l=0.1

算出来的结果很小，是不是漏掉了什么啊，是不是还和那个Q0有关啊，我看了关于adaptive CUSUM的那篇英文文献，也没看出哪里的问题，想麻烦您再给看一下，一直麻烦您，我都不知道说什么好了，谢谢您啦！！！

为更精准,我改设m1=37,m2=37
aa=log10(1+2*h*h*400...)    应该是aa=log(1+2*h*h*400...)
One-sided 是不是应该取PM(1)=1;   ??

when delta=0,  ARL=399.8709
when delta=1,  ARL=    8.7482

不好意思，我又要麻烦您了~~附件里面我上传了一篇文献，是关于adaptive CUSUM 控制图用来监控方差变化的，我参考Appendix B计算链长，在无方差变动的情况下，分别取控制限H=1.0274，Qt的控制限参考里面的公式取L=1.0274，n=5, δmin=1.2,l=0.1.  


下面是我计算链长的程序，我把上面参数值带进去计算，应该算出来的链长值为初始值200啊，但是算出来不对，我不知道是哪里出现问题，所以能不能麻烦您给看看？又打扰您了~~~


function y=ARLacusum(w,m,q,p)
m1=15;m2=15;
R=zeros(m1*m2,m1*m2);
for i=0:1:m1-1
for j=0:1:m2-1
     for k=0:1:m1-1
     for l=0:1:m2-1

h=max(m,(l+0.5)*q);
         aa=log(h)/(1-1/h);

     a1=(k-i-0.5)*w/(1-1/h)+aa;
     a2=(k-i+0.5)*w/(1-1/h)+aa;
     b1=(l-(1-p)*(j+0.5))*q/p;
     b2=(l+1-(1-p)*(j+0.5))*q/p;

     if k==0&&l==0
         temp=min(a2,b2);
         T=gamcdf(temp,2,0.5);
         R(i*m2+(j+1),k*m2+(l+1))=T;
     end

if k==0&&l~=0
          temp=b1;
          T1=gamcdf(temp,2,0.5);
          temp=max(b1,min(a2,b2));
          T2=gamcdf(temp,2,0.5);
          R(i*m2+(j+1),k*m2+(l+1))=T2-T1;
         end

if k~=0&&l==0
               temp=a1;
               T1=gamcdf(temp,2,0.5);
               temp=max(a1,min(a2,b2));
               T2=gamcdf(temp,2,0.5);
                R(i*m2+(j+1),k*m2+(l+1))=T2-T1;
          end

if k~=0&&l~=0
                temp=min(b2,max(a1,b1));
                 T1=gamcdf(temp,2,0.5);
                temp=max(b1,min(a2,b2));
                T2=gamcdf(temp,2,0.5);
                R(i*m2+(j+1),k*m2+(l+1))=T2-T1;            
                 end;
             end;
         end;
     end;
     end;
PM=0*ones(1,m1*m2);
PM(1)=1;
y=PM*inv(eye(m1*m2,m1*m2)-R)*ones(m1*m2,1);

能否告知底下这些参数你设置多少
w,m,q,p,L
尤其是L

if Conditionally execute statements
似乎跟上次不同了
%%%%%%%%%
for i=0:m1-1
for j=0:m2-1
  for k=0:m1-1
    for l=0:m2-1
     if k==0
         statements 1   
     elseif k~=0
         statements 2     
     end  %end if
     
    end
  end
end
end

%%%%%%%%%%%%%%%%

gamcdf(temp,2,0.5)
为何shape parameter取 2
scale parameter取 0.5

我的L值是参考文献里面Appendix B里面的那个公式：L=1+5*squre（l/(2-l)*(2/n)）参数设置为：w=2*h/2*m1-

其中，n=5  那个gamcdf(x,a,b)
Zt 服从gamma 分布，程序里面我的参数取为：


文献里面St2=∑(X-μ0)2/n； 那么这个a和b的值是不是取为：
这个地方我也不是很明白~~麻烦您再给看一下，谢谢啦！！！

则m1=15;m2=15; ARL= 199.5398
   m1=39;m2=39; ARL= 209.8057
与文献相差10左右,还需再CHECK一下
         

我想再问一下，在这个程序里面，还有那篇Adaptive CUSUM关于均值的文献中，在计算链长时都提到了 Q0 的确定，我想问一下，计算链长的程序和 Q0 的值有关系吗？从程序里面看是没有关系的，难道和 Qt 的取值范围有关？另外，不同的改进 CUSUM 控制图，关于 Qt  的取值范围是不是有计算公式啊？您知道怎么计算的吗？

您好！。。。

T1=normcdf(temp,0,1);  => T1=normcdf(temp,1,1);

T1=normcdf(temp,1,1)这个是表示偏移量为1 的时候程序里面该修改的地方吗？如果是计算方差的变动量，比如当方差的变动量为2时，是不是又改为
T1=normcdf(temp,0,2）


我这个计算CUSUM控制图链长的程序是对的吧？谢谢您的指导！！！

您好！一直麻烦您好久，我想您也是做这方面的研究吧。我现在写论文想写关于方差方面的文章。对于上次向您请教的关于方差的adaptive CUSUM control chart 中，文献里面采用的是St2=∑(X-μ0)2/n ，在这里我采用的对数变换，即
Yt=ln(St2/σ2)
St2=∑(X-μ0)2/n-1
这样计算控制限的程序我把它改写为：
function y=optif(x,ARL0)
m1=15;m2=15;
n=5;
p=x(1);
m=1.1;
w=(2*x(2))/(2*m1-1);


%L=1+5*sqrt(p/(2-p)*(2/n));


L=1;
q=L/m2;
R=zeros(m1*m2,m1*m2);
for i=0:1:m1-1
  for j=0:1:m2-1
    for k=0:1:m1-1
     for l=0:1:m2-1

             end;
         end;
     end;
   end;
PM=0*ones(1,m1*m2);
PM(1)=1;
y=abs(ARL0-PM*inv(eye(m1*m2,m1*m2)-R)*ones(m1*m2,1));


我首先计算控制线，
ARL0=500; x0=[0.1,1.2];  lb=[0.1,1];   ub=[0.1,2];  [x,fval]=fmincon(@(x)optif(x,ARL0),x0,[],[],[],[],lb,ub)
但是不知道为什么区间取为多少总是初始值，为什么得不到优化参数呢？不知道您对这方面了解吗？可以麻烦您帮我看看这个程序吗？又要麻烦您了~~谢谢啦！！

lambda 初始值影响不大
ARL0=200; 500;1000;
可设相同0.1,lb=0.05,ub=0.2
但hw初始值影响很大大
请注意我的设置
以13 楼程序而言
求得的值相当理想
%ARL0=200;
x0=[0.1;1.1];  
lb=[0.05;1.0];   
ub=[0.2;1.3];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)
x =
    0.1336
    1.1587

fval =

  3.9378e-005
%%%%%%%%%%%%%%%%%
ARL0=1000;

x0=[0.1;1.8];  
lb=[0.05;1.7];   
ub=[0.2;1.9];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)
x =
    0.0800
    1.8846

fval =

  7.6930e-004

%%%%%%%%%%%%%%%%%
ARL0=500;

x0=[0.1;1.5];  
lb=[0.05;1.4];   
ub=[0.2;1.65];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)

x =
    0.0539
    1.4748

fval =

  3.1012e-007

%%%%%%%%%%%
%%%%%%%%%%%
以18 楼程序而言
求得的值相当离谱
必须重新检查

ARL0=500;
%x0=[0.1;1.2];  
%lb=[0.1;1];   
%ub=[0.1;2];  
x0=[0.1;1.5];  
lb=[0.05;1.4];   
ub=[0.2;1.65];  
[x,fval]=fmincon(@(x)optif(x,ARL0),x0,[],[],[],[],lb,ub)

x =

    0.0500
    1.5000


fval =

  493.8829

%%%%%%%%
%%%%%%%%
为何要采对数变换??

底下这样有对数变换吗?
temp=b1;
u=exp(temp);
T1=gamcdf(u,2,0.5);

%%
Appendix A: Derivation of the CUSUM Dispersion Chart
            S(t)^2 has a gamma distribution with shape parameter
            alpha= n/2 and scale parameter beta=2*sigma^2/n.

你采对数变换，
Yt=ln(S(t)^2/(sigma)^2)
那Yt服从甚么分布??

19# epoh

epoh老师，您好！
      我用您的方法进行计算，为什么会提示以下错误：

>> %ARL0=200;
x0=[0.1;1.1];  
lb=[0.05;1.0];   
ub=[0.2;1.3];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)
??? Error: "opt" is not the name of a function nor a variable,
but is used in an anonymous function either at the prompt or in the argument of EVAL.
>> ARL0=1000;
x0=[0.1;1.8];  
lb=[0.05;1.7];   
ub=[0.2;1.9];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)
??? Error: "opt" is not the name of a function nor a variable,
but is used in an anonymous function either at the prompt or in the argument of EVAL.
>> ARL0=1000;
x0=[0.1;1.8];  
lb=[0.05;1.7];   
ub=[0.2;1.9];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)
??? Error: "opt" is not the name of a function nor a variable,
but is used in an anonymous function either at the prompt or in the argument of EVAL.
>> %ARL0=200;
x0=[0.1;1.1];  
lb=[0.05;1.0];   
ub=[0.2;1.3];  
[x,fval]=fmincon(@(x)optif(x,ARL0),x0,[],[],[],[],lb,ub)
Warning: Large-scale (trust region) method does not currently solve this type of problem,
switching to medium-scale (line search).
> In fmincon at 260

非常感谢！

fmincon:
  Find minimum of constrained nonlinear multivariable function
(类似R:optim(),gauss:library cml,但功能强大)

求解garch model 系数的时候就用到
[parameters, LLF,..] =
fmincon('garchlikelihood', [omega ; alpha ; beta] ,...);

'garchlikelihood'是我们必须提供的 Objective function

在这里[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)

'opt'就是13楼所定义的Objective function


opt.m


  

opt.rar
大小:(534 Bytes)

 马上下载

本附件包括：

• opt.m

  

19# epoh

function y=opt(x,ARL0)
m1=15;m2=15;   
n=5;
p=x(1);   %p=0.1
m=1.2;
hw=x(2);  %hw=1.0274
w=(2*hw)/(2*m1-1);
L=1+5*sqrt(p/(2-p)*(2/n));
%L=1+5*(sqrt(2*p/((2-p)*n)))
q=L/m2;
R=zeros(m1*m2,m1*m2);
for i=0:1:m1-1
for j=0:1:m2-1
     for k=0:1:m1-1
     for l=0:1:m2-1
h=max(m,(l+0.5)*q);
         aa=log(h)/(1-1/h);
     a1=(k-i-0.5)*w/(1-1/h)+aa;
     a2=(k-i+0.5)*w/(1-1/h)+aa;
     b1=(l-(1-p)*(j+0.5))*q/p;
     b2=(l+1-(1-p)*(j+0.5))*q/p;

if k==0
          temp=b1;
          T1=gamcdf(temp,2.5,0.4);
          temp=max(b1,min(a2,b2));
          T2=gamcdf(temp,2.5,0.4);
          R(i*m2+(j+1),k*m2+(l+1))=T2-T1;
function y=opt(x,ARL0)
m1=15;m2=15;   
n=5;
p=x(1);   %p=0.1
m=1.2;
hw=x(2);  %hw=1.0274
w=(2*hw)/(2*m1-1);
L=1+5*sqrt(p/(2-p)*(2/n));
%L=1+5*(sqrt(2*p/((2-p)*n)))
q=L/m2;
R=zeros(m1*m2,m1*m2);
for i=0:1:m1-1
for j=0:1:m2-1
     for k=0:1:m1-1
     for l=0:1:m2-1
h=max(m,(l+0.5)*q);
         aa=log(h)/(1-1/h);
     a1=(k-i-0.5)*w/(1-1/h)+aa;
     a2=(k-i+0.5)*w/(1-1/h)+aa;
     b1=(l-(1-p)*(j+0.5))*q/p;
     b2=(l+1-(1-p)*(j+0.5))*q/p;

if k==0
          temp=b1;
          T1=gamcdf(temp,2.5,0.4);
          temp=max(b1,min(a2,b2));
          T2=gamcdf(temp,2.5,0.4);
          R(i*m2+(j+1),k*m2+(l+1))=T2-T1;
  elseif k~=0
          temp=min(b2,max(a1,b1));
          T1=gamcdf(temp,2.5,0.4);
          temp=max(b1,min(a2,b2));
          T2=gamcdf(temp,2.5,0.4);
          R(i*m2+(j+1),k*m2+(l+1))=T2-T1;            
          end;  %end if
             end;
         end;
     end;
   end;
PM=0*ones(1,m1*m2);
PM(1)=1;
x(1)
x(2)
PM*inv(eye(m1*m2,m1*m2)-R)*ones(m1*m2,1)
y=abs(ARL0-PM*inv(eye(m1*m2,m1*m2)-R)*ones(m1*m2,1))

epoh老师，您好！
p=0.1是怎么来的？有什么用？
hw=1.0274是怎么来的？有什么用？

>> ARL0=1000;
x0=[0.1;1.8];  
lb=[0.05;1.7];   
ub=[0.2;1.9];  
[x,fval]=fmincon(@(x)opt(x,ARL0),x0,[],[],[],[],lb,ub)
Warning: Large-scale (trust region) method does not currently solve this type of problem,
switching to medium-scale (line search).
> In fmincon at 260
ans =
    0.1000

ans =
    1.8000

ans =
  753.0236

y =
  246.9764

ans =
    0.1000

ans =
    1.8000

ans =
  753.0235

y =
  246.9765


ans =
    0.1000

ans =
    1.8000

ans =
  753.0236

y =
  246.9764

ans =
    0.0500

ans =
    1.9000

ans =
  791.7859

y =
  208.2141

ans =
    0.0500

ans =
    1.9000

ans =
  791.7865

y =
  208.2135

ans =
    0.0500

ans =
    1.9000

ans =
  791.7859

y =
  208.2141

ans =
    0.0941

ans =
    1.9000

ans =
  954.6719

y =
   45.3281

哈哈!zhangtao兄,
p=0.1是怎么来的？有什么用？
hw=1.0274是怎么来的？有什么用？
这可能有劳你抽空看一下
楼主10楼提供的文献了.

您好！我是想以这篇文献再稍微修改一下，因为要写论文，不能全部照搬过来，嘿嘿~~关于对数变换，我也是参考其他的文献，也是该作者，分别研究改进CUSUM和改进EWMA控制图的均值和方差，其中方差方面他提出了几种变换，我在帖子里面传了一篇文献，我是按照这篇文献里面介绍的改写的，我想把这种思想应用到基于方差的改进CUSUM控制图中，但是不知道哪里出了问题，出现了上面那种控制限计算不出来的现象，想麻烦您给您看一下~~谢谢您了！！

我知道你问题出在那里了

page 16/18

transition probability,a1,a2

都已不同(没了b1,b2)

且需定义inverse of the Huber score function

这部分麻烦你先列一下

我再帮你编程

老师，您好！我传的这篇文献是想让您看一下我想应用这篇文献里面的这个对数变换的思想，我的意思是想把对数变换的思想应用到10楼我上传的那篇 adaptive CUSUM控制图中，那一篇文献中采用的统计量是 Zt=St2/σ2  此时Zt服从gamma分布；现在我想把这个统计量换为对数变换，即统计量为：
Yt=ln(St2/σ2) ，其中，St2=∑(X-μ0)2/n-1
意思就是采用这个统计量的前提下，怎么计算10楼我上传的那篇 adaptive CUSUM控制图中的控制限以及对其链长的分析，
按照24楼我上传的那篇文献，我自己改编的程序就是之前18楼里面那个，但是我不知道哪里出现问题，找不到控制限，链长更没法分析，所以想麻烦老师给看一下，非常感谢您每次热心给我解答问题，又打扰您了！！谢谢您！！

刚要改程序仔细看

才发现这是一维的

P(i, j) = Pr {Wt in state j |Wt−1 in state i}

        = Pr {a1 < Yt ≤ a2},

难道你要回头做一维

老师，您好！我不是要做一维的，我给您看这篇文献是在帖子里面您说我为什么做对数变换，其实我还是针对我在10楼里面上传的那篇 加权 CUSUM控制图，那一篇文献里面采用的统计量为：Zt=St2/σ2 计算程序您也帮我写了就是下面这个：




PM=0*ones(1,m1*m2);
PM(1)=1;
y=PM*inv(eye(m1*m2,m1*m2)-R)*ones(m1*m2,1)


这个程序里面这个转移概率其中一部分为
Pr={b1< Zt <max[b1,min(a2,b2)]}


里面这个Zt=St2/σ2，服从gamma分布，现在我想把Zt换为
Zt=ln(St2/σ2)
就是转换一下统计量，在Zt=ln(St2/σ2)的情况下怎么再计算链长和控制限？我自己在上面这个程序里面改了一下但是无法得到控制限，所以想麻烦您给看一下，这种情况下怎么改写程序~~谢谢您了！！！

你有没发现table4,table5

h的值差异那么大

但Zero-state ARL却接近,

这篇文献我做了ln(S2),及Xs(other normalizing transformation)

程序就是依Appendix A.Appendix B.

但结果差异颇大

你有没类似的第二篇文献可供比较.

%%%%%%%%%%%

page 10/18,Table 4.Zero-state ARL

n=5
lambda=0.1
r=1.5
h=0.2881
ARL=199.8

%%%%%%%%%%%%

page 14/18,Table 5.Zero-state ARL
n=5
lambda=0.1
r=1.5
h=0.6652
ARL =199.99

老师，您好！我仔细看了一下那篇文献，是存在这个问题。