急求MATLAB指点！！！

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

function y=QQ(L,r,k,UCL,LCL)
d=(UCL-LCL)/k
m1=0;
for i=-k/2+1:1:k/2
m1=m1+1;
k1=1;
m=k/2;
for j=-m+1:1:m
temp=((j*d)-(1-r)*(i-0.5)*d)/(r);
T1=normcdf(temp,0,1);
temp=((j-1)*d-(1-r)*(i-0.5)*d)/(r);
T2=normcdf(temp,0,1);
R(m1,k1)=T1-T2;
k1=k1+1;
end;
end;
PM=0*ones(1,2*m);
PM(k/2)=1;
SS=PM*inv(eye(2*m)-R)*ones(2*m,1);
y=SS;

我定义了function y=QQ(L,r,k,UCL,LCL)这个函数，但是运行时最终结果和L的值无关，这是用MATLAB解决马尔科夫链转移概率问题，由于对MATLAB了解较少，所以一直没搞清楚，希望有这方面的高手给指点一下，我这个程序问题出在哪了？万分感谢！！！！

扫码加我 拉你入群

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

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

分享0 收藏2 回帖

关键词：MATLAB matla atlab Atl Lab MATLAB 指点 马尔科夫链；MATLAB

已经给了UCL,LCL,可以不用给L.

假设ARL=500,则由Genetic Algorithms
可得optimal parameters:(r=0.13,L=2.88).
以上数据,带入你的 QQ function
得 ARL = 499.1488

%%%%%%%%%%%%
r=0.13;
L=2.88;
UCL= L*sqrt(r/(2-r)); %  0.7593524
LCL=-L*sqrt(r/(2-r)); %- 0.7593524
k=50;
ARL=QQ(r,k,UCL,LCL)

%%%%%%%%%%%%%%
function ARL=QQ(r,k,UCL,LCL)
d=(UCL-LCL)/k;
m1=0;
for i=-k/2+1:1:k/2
  m1=m1+1;
  k1=1;
  m=k/2;
   for j=-m+1:1:m
     temp=((j*d)-(1-r)*(i-0.5)*d)/(r);
     T1=normcdf(temp,0,1);
     temp=((j-1)*d-(1-r)*(i-0.5)*d)/(r);
     T2=normcdf(temp,0,1);
     R(m1,k1)=T1-T2;
     k1=k1+1;
  end;
end;
PM=0*ones(1,2*m);
PM(k/2)=1;
ARL=PM*inv(eye(2*m)-R)*ones(2*m,1);

谢谢您的指点，您后面写的是相关程序吗？我运行ARL=QQ(0.13,50,0.759352,-0.753524)，怎么显示LCl这个变量不存在呢？因为对ＭＡＴＬＡＢ了解比较少，可能问题有点幼稚，另外，我还想请问一下，如果由已知链长ＡＲＬ求解控制限Ｌ，也是基于马尔科夫链求解吗？使用这个程序可以计算吗？我要解决的问题是用马尔科夫链求解统计过程控制中的控制图控制界限问题，只是看了一些相关方面的论文，只说了结果，现在有点糊涂，这个程序是用来计算ＡＲｌ还是我可以假定一个ＡＲＬ的值求解控制限，您可以再给指点一下嘛？谢谢啦！！！！

Q1:
ARL=QQ(0.13,50,0.759352,-0.753524)，怎么显示LCl?

QQ.m 请换成2楼我改的这个:
%%%%%%%%%%%%%%
function ARL=QQ(r,k,UCL,LCL)
d=(UCL-LCL)/k;
m1=0;
for i=-k/2+1:1:k/2
  m1=m1+1;
  k1=1;
  m=k/2;
   for j=-m+1:1:m
     temp=((j*d)-(1-r)*(i-0.5)*d)/(r);
     T1=normcdf(temp,0,1);
     temp=((j-1)*d-(1-r)*(i-0.5)*d)/(r);
     T2=normcdf(temp,0,1);
     R(m1,k1)=T1-T2;
     k1=k1+1;
  end;
end;
PM=0*ones(1,2*m);
PM(k/2)=1;
ARL=PM*inv(eye(2*m)-R)*ones(2*m,1);
%%%%%%%%%%%%%
则ARL=QQ(0.13,50,0.759352,-0.753524)
可得 ARL = 483.9766
QQ.zip (326 Bytes) 本附件包括：

• QQ.m

2010-9-23 20:15:02 上传

Q2:
你应该很清楚你定义的function QQ
是应用markov chain来求解ARL
当然不能求 L

Q3:
ARL的计算方法大致有三种方法：

马氏链法,积分方程法,和随机模拟法.

你是采马氏链法

而 R package "spc"

steady-state ARL是采integral equation.

R package "spc" 有很多function可用

也值得你参考看看
    spc.pdf (240.6 KB)

2010-9-23 20:15:02 上传

R package "spc" provides
xewma.ad   steady-state ARLs of EWMA charts
xewma.arl (zero-state) ARLs of EWMA charts
xewma.crit critical values of EWMA charts
底下是在R运行结果:
###steady state
xewma.ad(l=0.13,c=2.88,mu1=0,mu0=0,zr=0,sided="two",limits="fix",r=40)
#498.8069
xewma.crit(l=0.13,L0=500,mu0=0,zr=0,hs=0,sided="two",limits="fix",r=40)
#2.876541  

呜呜，matlab好难啊

凑个热闹，你说的函数声明是：
function y=QQ(L,r,k,UCL,LCL)
而在程序运行时，竟然现实变量LCI没有定义，换成小写；
lcl在程序内部有部分被定义为lci。
解决的方法是在edit qq 打开程序文本，用查找和替换，将程序里的lci 替换为lcl。

语法问题。

你好！谢谢上次你对我的帮助，我想请问一下，那个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和算法了解太少了，我想问一下我这个思路对不对，能不能再麻烦您给指点一下？