CUSUM1,EWMA2,的作者都是Sven Knoth

也就是 R package "spc"的作者.

如果你是要做ACUSUM

而且是用Markov chain,

你可能漏掉了一篇重要文献

Chang,T.C.,Gan,F.F.,1995.A cumulative sum control chart for monitoring process variance

ACUSUM:

目前完成Markov chain程序的只有

10楼那篇ACUSUM Procedure for Signaling Process Variance

要转成 ln(S^2)可能需要参考Chang,T.C.,Gan,F.F.,1995

AEWMA:

AEWMA Control Chart for Monitoring Process Variances(ln(S^2))

根据page 16/18的程序,跑出的结果,有些怪异

但程序应该无误,可能是参数设置转换问题

我把程序先传给你参考

请检查,试调整,及试跑结果.

由于程序尚未成熟,

为避免不必要误解.

我设了pw

pw就在你的短信息.

lns2.txt,Xs.txt

  

谢谢您的指导，我先找找那篇文献再跟您交流~~谢谢啦！！！

CUSUM1.pdf page 12/14

Table 4, 就是依据Table 5 of Chang and Gan (1995)

in-control ARL value 100

内有四种方法

  (1).CUSUM-S2

  (2).CUSUM-ln(S2)

  (3).EWMA-S2

  (4).EWMA-ln(S2)

%%%%%%

(1)CUSUM-S2

   可用Hawkins＇anyarl.exe算出 ARL 99.92

   N=5  (df=4)

   kh=1.285

   hh=2.922

   anyarl.exe.rar

      

(2)CUSUM-ln(S2)

   等你抓到Chang and Gan (1995)文献后

   再看参数如何转换

(3)EWMA-S2

   可用R package "spc",function sewma.arl()算出 ARL 99.99888

   lambda=0.15

   c=2.4831

   ###

   df=4

   cu <- 1 + 2.4831*sqrt(lambda/(2-lambda))*sqrt(2/df) #1.4999964

   cu的算法请参考:EWMA2.pdf page 2/12 formula(3)

   sewma.arl(lambda, 1, cu, 1, df,sided="Rupper", r=100)#99.99888

   

(4)EWMA-ln(S2)

   lambda(ln)=0.28

   c(ln)=1.4085

   ARL=99.9547

   EWMA_ln(S2).rar

   

   pls aslo see EWMA1.pdf  page 11/13

暂时保存。。

只要适当的转换

EWMA_S2,EWMA_ln(S2),

都可轻易的算出ARL.

给定lambda,L

%%%%EWMA_ln(S2)  ARL=199.7620
n=5
lambda=0.32
L=1.60861
sigmay=sqrt((2/(n-1))+(2/(n-1)^2)+(4/(3*(n-1)^3)-(16/(15*(n-1)^5))))  
h=L*sqrt(lambda/(2-lambda))*sigmay
......
temp=b2;
u=exp(temp);
R(i,j)=gamcdf(u,2.5,0.4);
......

%%%%EWMA_S2      ARL= 200.8928
lambda=0.16;
L=2.96352;     
h=1+L*sqrt(lambda/(2-lambda))*sqrt(2/n)
......
temp=b2;
u=(temp);
R(i,j)=gamcdf(u,2,0.5);
.......

EWMA_S2.rar

   

EWMA_S2.rar
大小:(419 Bytes)

 马上下载

%%%%

同理

既然10楼WACUSUM_S2程序都已完成

依你的专业应该能顺利转换WACUSUM_ln(S2)

PS:

不过还是有隐忧

因为都没看过有关

ACUSUM-ln(S2)的文献

或是CUSUM-ln(S2)的文献

老师您好！根据您的指点，我做了关于 EWMA-S2 和 EWMA-lnS2 的控制限的搜索以及链长的分析，都求解出来了，谢谢您的指点。我也搜索了好多，确实没有相关CUSUM-lnS2方面的文献，但是根据常规CUSUM控制图计算均值的程序，我就是运用了同样的对数变换，还是求解出了控制限，但是不知道为什么，应用到ACUSUM的链长计算程序中就出现了问题。我想问一下，Qt=（1-lamda）Qt-1 +lamda*Zt    其中，Zt=ln(St2/σ2)  这个二维马尔科夫链尺度范围为： 0<Wt<h, 0<Qt<L
h是其控制限，这里Qt的控制限范围我按照公式 L=E(Zt)+5*sqrt(lambda/(2-lambda))*sigmay  
我只是大概算了一下，n=5,sigmay=0.8
我想问一下，应用到这个二维马尔科夫链中，是因为这个L值可能不是很准确，还是加了对数变换之后，相应里面这个转换矩阵的计算公式有所改变，为什么控制限的搜索值老是等于初始值，计算出来的误差相当大，我自己也再研究一下，想麻烦您看一下，是不是里面这个转换概率加了对数变换之后有什么改变呢？谢谢您啦！！！

EWMA_S2,       两个参数, lambda,L   一维

EWMA_ln(S2),  两个参数, lambda,L   一维

AEWMA_S2,       ??

AEWMA_ln(S2), 三个参数, lambda,r,h   一维

ACUSUM_ln(S2),  ??                            ??

ACUSUM_ln(S2),比较麻烦的是,

没有文献可参考比对,

来验证程序的正确性.

如果方便的话看能不能找到:

Brook,D.and Evans,D.A.(1972)An Approach to the Probability Distribution of CUSUM Run Length

Howell, D.K., 1987. A CUSUM scheme for the control of process variance

Srivastava, M.S., Chow, W., 1992. Comparison of the CUSUM procedure with other procedures

          that detect an increase in the variance and a fast

          accurate approximation for the ARL of the cusum procedure

老师，您好！我把adaptive CUSUM控制图那篇研究均值的程序改为计算方差的，也就是这篇帖子里面7楼那个程序，中间Qt的取值范围我按照公式
L=5*sqrt(lambda/(2-lambda))*sigmay

参数设置分别为：lamda=0.1  dertamin=1.1 ARL0=200
控制限搜索出来为：0.2325当lamda为0.2时，控制限为0.2266，相差不大，不知道计算的对不对，
并且我还做了σt2/σ2分别偏移量为1.1,1.2，1.5,2 等情况下的链长分析，感觉好像都没有常规CUSUM控制图做出来效果好，不知道这个控制限计算的对不对？还有那篇关于均值的adaptive CUSUM控制图的程序不知道能不能像AEWMA控制图那样经过对数变换就可以应用到方差方面的研究上面了，您说的那几篇文献我正在找，等下载到了我再给您交流，谢谢您了！！

直接用差分方程不能解吗

dertamin=1.1, ARL0=200,

lamda=0.1,  控制限搜索出来为：0.2325

当lamda为0.2时，控制限为0.2266，

相差不大，不知道计算的对不对?

应该是不对.

我说过现在最大的困扰是

没有文献可对照

无法验证程序的正确性

我个人觉得ACUSUM,求ARL

既然给定了L,hw,

那么最重要的就是

h,a1,a2,b1,b2公式的正确性了

正如你由mean转为S2时,

不就是更改了h,a1,a2,b1,b2

您分析的有道理，我下载到几篇关于CUSUM-S2的文献，现在传上去您看一下，我也自己再研究一下，谢谢您了！！！

A comparison of control charting....
Table 1 ARL 200
1.CUSUM_ln(S2), k=0.068, h=2.66
2.CUSUM_Xs,      k=0.38,  h=4.28
The performance of control........
Table X ARL 500
  CUSUM_S2,      k=1.46,   h=3.72787

3个程序都OK

CUSUM_S2,Xs,ln(S2).rar

目前就只剩ACUSUM_ln(S2),

缺乏数据可以比照

既然这是你的重头戏

我看教科书,...,

都得尽量找找看罗.

毕竟你是做ACUSUM非AEWMA

哈哈!这是我最近才发现的.

谢谢老师的指点，之前我可能没说清楚，我是研究ACUSUM，并不是AEWMA，之前还麻烦您给写了关于AEWMA的程序，谢谢您啦！！我自己也会努力研究的，关于ACUSUM-lnS2的文献好像特别少，我也找找看，谢谢啦！！祝您天天快乐！！

仔细看了三篇adaptive 文献

7楼:A Markov Chain Model for the Adaptive CUSUM....

10楼:An Adaptive CUSUM Procedure for Signaling Process Variance

24楼:An adaptive EWMA control chart....

所选用的weight function,score function都不相同

就以10楼 WACUSUM_S2而言

m=1.2

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;

适用于WACUSUM_S2,然不适用于WACUSUM_ln(S2)

因为几乎exp(temp)之后都是1,

造成 T2-T1=0,所以ARL非常小.

adaptive应该不能如此转换.

建议你想办法拿到最基础的这篇文献

Brook and Evans (1972),An Approach to the Probability Distribution of the CUSUM Run Length

再来参考修改

另这篇文献也试试

Shu,L.,Jiang,W.and Wu,Z.,Adaptive CUSUM Procedures with Markovian Mean Estimation 2008

哦，原来是这个样子啊，那篇最基础的文献我再想想办法吧，可能有点太老了，那篇ACUSUm-C的文献我下载到了，它是把AEWMA的思想又加进来了，我传上去您看看，谢谢您的热心指导！！！ACUSUM-c.pdf

补充:

10楼,10年加权page 2/17,公式(3)

the increment of the CUSUM statistic at time t

推导公式就在Appendix A: Derivation of the CUSUM Dispersion Chart

这是根据S2 gamma distribution with shape parameter alpha=n/2

           and scale parameter beta=2*sigma^2/n.

公式(10) page 3/17

  (1-/delta) can be viewed as certain weight functions

             on the deviations (Zt-kt)

所以如果你是要做WACUSUM_ln(S2),

也要有类似的推导

找出正确的weight function & reference value

老师，您好！您说的有道理，我也觉得这个ACUSUM可能是变量改了之后转移概率的求法上面有了较大改动，而不是像其他控制图那样直接做转换就可以了，我再仔细研究研究，有什么问题再跟您交流，祝您天天开心！！！

      老师，您好！如果我做WACUSUM-lns2 chart，令统计量为Wt={0，Wt-1+Zt-kt} ，其中，Zt=St2/σ2，   常规CUSUM控制图的统计量可以写为Ut={0，Ut-1+Yt-k},
      也就是我把 常规控制图中的delta 做指数加权变换，也就是adaptive CUSUM的思想，我觉得这个 weigh function 应该是合理的吧，如果这样做WACUSUM-lns2 chart，这个转移概率是不是也要改呢？我算了一下，还是控制限搜索不出来，误差很大，想麻烦老师给看一下，谢谢您了！！

这个地方是不是还是像前面那样经对数变换T2-T1会很小呢，我自己也再仔细算一下，想麻烦老师先看看这个思路对不对，谢谢您了！！！

10年_加权,page 3/17,公式(10)是WACUSUM

page 3/17,公式(8)是ACUSUM,少掉了weight functions (1-1/delta).

你的程序用于 page 4/17,Table 1

ACUSUM_S2,是正确的

ARL= 495.4425

%%%%%%%%

m1=37;m2=37;   

n=5;

p=0.1

m=1.21

hw=6.231

w=(2*hw)/(2*m1-1);

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

q=L/m2;

h=max(m,(l+0.5)*q);

aa=log(h)/(1-1/h);

a1=(k-i-0.5)*w+aa;

a2=(k-i+0.5)*w+aa;

b1=(l-(1-p)*(j+0.5))*q/p;

b2=(l+1-(1-p)*(j+0.5))*q/p;

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

但用于ACUSUM_ln(S2),依旧是T2-T1很小

这个地方的转换看似简单,

但可惜没抓到窍门,

我看还是找文献参考比较有学理基础,

至少你在写论文时才有所凭据.

老师，您好！这方面的文献好像比较少，我先找找看，自己也再考虑一下，再跟您交流，谢谢您的指导！！祝您天天开心！！

相关文献很难找吧!

稍微修改原程序,

可以找到optimal lambda,hw,L;

然带入sigma1/sigma=1.1,数据不对

还得继续找出原因.

请问有人会写(2008)Adaptive CUSUM procedures with Markovian mean estimation的程序吗?

马尔科夫过程(MarKov Process)是一个典型的随机过程。设X(t)是一随机过程，当过程在时刻t0所处的状态为已知时，时刻t(t>t0)所处的状态与过程在t0时刻之前的状态无关，这个特性成为无后效性。无后效的随机过程称为马尔科夫过程。马尔科夫过程中的时同和状态既可以是连续的，又可以是离散的。我们称时间离散、状态离散的马尔科夫过程为马尔科夫链。马尔科夫链中，各个时刻的状态的转变由一个状态转移的概率矩阵控制。  
二维隐马尔科夫模型  
例如，采样窗是垂直方向的，自然想到是不是也可以在水平方向建立状态这就是二维隐马尔科夫模型，但是二维HMM的训练和识别算法的复杂性，使用不是很理想，嵌入式隐马尔科夫模型由一系列超状态组成，每个超状态又包含若干状态，称为嵌入状态，超状态反应其中的一维信息，嵌入状态反映另一维的信息，但是由于超状态内的状态之间没有状态的转移，所以不是真正的二维。只能看作是一个简化的二维。  

人脸的超状态模型就是刚才的从上到下的五个状态（前面的HMM中的五个状态），在各个超状态之间增加水平信息。转移关系还是从一个超状态到另一个超状态。嵌入的水平状态数为3，6，6，6，3用一个或多个分量的混合高斯密度函数表示。  

观测向量的提取，首先把人脸分为图像快。然后取图像块的灰度值或者变换系数组成一个观测向量，图像快采用遍历的方法进行采样，就是从上到下，从左到右，来获取图像的采样快。2D-DCT的低频分量。