求解二维马尔科夫链

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我应用二维马尔科夫链[0,c]*[a,b]来构建一个模型，其中：0<=Z<=c;a<=T<=b.把Z的受控区间划分为m1个状态，每个状态长度为f=2c/2m1-1;同样的，T的受控区间也划分为m2个状态，每个状态的长度为v=2（b-a）/2m2-1矩阵R中的每一个元素转移概率p(i,j)(k,l)可表示为：p(i,j)(k,l)=pr{Z(t)处于状态k,T(t)处于状态l/ Z(t-1)处于状态i,T(t-1)处于状态j }=pr{min[b2,max(a1,b1)]<X<max[b1,min(a2,b2)]}k不等于0，l不等于0；Pr{a1<X<max[a1,min(a2,b2)]}k不等于0,l等于0；Pr{b1<X<max[b1,min(a2,b2)]}k等于0, l不等于0；Pr{X<min(a2,b2)}k等于0, l等于0；其中，a1,a2,b1,b2是变量，涉及到状态空间长度等，公式比较麻烦，就不写了。这个最后的R阵是不是N*N的方阵，其中N=m1*m2，这个矩阵怎么编写出来的，我不知道那个p(i,j)(k,l）怎么编译，最后生成的到底是个什么样的矩阵？望诸位热心人帮忙给看一下，指点一下！！十分感谢！！！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

关键词：二维马尔科夫链 马尔科夫链 马尔科夫 马尔科 矩阵 二维马尔科夫链

Adaptive CUSUM Control Chart  ??

有现成的软件？

是的？您也研究这个吗？其实一直都是看的这个控制图，因为没有看到相关方面的资料，所以才从EWMA看起的。

假设m1=17, m2=25,
  N = 17 x 25
2011.05.09 revised
  matrix R = 425 *425

如果是Adaptive CUSUM Control Chart
初看结构应该是:
4 for loop + if Conditionally execute statements
statements 就由你写上
如果不会放进p(i,j)(k,l）
届时我再帮你写进.

%%%%%%%%%
for i=0:m1-1
for j=0:m2-1
  for k=0:m1-1
    for l=0:m2-1
     a1=...
     a2=...
     b1=..
     b2=...
     if k==0 && l==0
         statements 1
     elseif k==0 && l~=0
         statements 2
     elseif k~=0 && l==0
         statements 3
     elseif k~=0 && l~=0
          statements 4
     end  %end if   
    end
  end
end
end

function y=ARL()
m1=17;m2=25;
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*(ln(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*(ln(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);
     else if k==0&&l~=0
          temp=b1;
          T1=normcdf(temp,0,1);
          temp=max(b1,min(a2,b2));
          T2=normcdf(temp,0,1);
         else if k~=0&&l==0
               temp=a1;
               T1=normcdf(temp,0,1);
               temp=max(a1,min(a2,b2));
               T2=normcdf(temp,0,1);
             else if k~=0&&l~=0
                temp=min(b2,max(a1,b1));
                 T1=normcdf(temp,0,1);
                temp=max(b1,max(a2,b2));
                T2=normcdf(temp,0,1);


    i=i+1;
    j=j+1;
     k=k+1;
     l=l+1;
                 end;
             end;
         end;
  end;
end;

这四个循环我已经放进去了，那个P(i,j)(k,l)我还真不知道怎么放进去，还请您帮我指点一下，非常感谢您的热心，真是麻烦您了~~~

先确认几件事
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);