求解二维马尔科夫链

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

2011-7-28 10:32:50 上传

  

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楼提供的文献了.

epoh 发表于 2011-7-27 21:55
lambda 初始值影响不大
ARL0=200; 500;1000;
可设相同0.1,lb=0.05,ub=0.2

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

我知道你问题出在那里了

page 16/18

transition probability,a1,a2

都已不同(没了b1,b2)

且需定义inverse of the Huber score function

这部分麻烦你先列一下

我再帮你编程

epoh 发表于 2011-8-4 21:52
我知道你问题出在那里了page 16/18transition probability,a1,a2都已不同(没了b1,b2)且需定义inverse of th ...

老师，您好！我传的这篇文献是想让您看一下我想应用这篇文献里面的这个对数变换的思想，我的意思是想把对数变换的思想应用到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},

难道你要回头做一维

epoh 发表于 2011-8-5 18:28
刚要改程序仔细看才发现这是一维的P(i, j) = Pr {Wt in state j |Wt−1 in state i}        = Pr {a1  ...

老师，您好！我不是要做一维的，我给您看这篇文献是在帖子里面您说我为什么做对数变换，其实我还是针对我在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

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