求经验欧式似然统计量的功效-源程序

1、m=8;N=20 400;%;20 25 30 40 50 100 150 200 250 300 400cishu=1000;GX1=;GX2=;GX3=;GX4=;GX5=;%用于检验重新分割的%for k=1:length(N)gx1=0;TJ1=0;gx2=0;TJ2=0;gx3=0;TJ3=0;PP2=;PPerson=;T1=;T2=;PLE=;PQ=;gx1xin=0;TJ1xin=0;gx2xin=0;TJ2xin=0;gx3xin=0;TJ3xin=0;PP2xin=;PPersonxin=;PP2xin=0;T1xin=;T2xin=;PLExin=;PQxin=;for

2、i=1:cishuX=normrnd(0,1,N(k),1);%不同分布用于对比%Y=trnd(4,N(k),1);%Y=trnd(1,N(k),1);%Y=normrnd(0,2,N(k),1);%Y=normrnd(0,1.5,N(k),1);%Y=poissrnd(1,N(k),1);%Y=exprnd(1,N(k),1);%Y=betarnd(1,2,N(k),1);a1=min(X);b1=max(X);a2=min(Y);b2=max(Y);aa1=min(X);aa2=min(Y);p1=normcdf(aa1+(b1-a1)/m);p=p1; %P为备择假设的各个区间的理论概率

3、for j=2:mif j<ms=normcdf(aa1+j*(b1-a1)/m)-normcdf(aa1+(j-1)*(b1-a1)/m);elses=1-normcdf(aa1+(m-1)*(b1-a1)/m);endp=p;s;endpp=;for i=1:N(k)pp=pp;p;endP=reshape(pp,m,N(k); %p为理论概率生成的矩阵% S1=(b1-aa1)/m;S2=(b2-aa2)/m;xduandian=linspace(a1,b1,m+1);yduandian=linspace(a2,b2,m+1);O12=zeros(m,N(k);O22=zeros(

4、m,N(k);for j=1:N(k)for i=1:mif (xduandian(i)<=X(j)&(xduandian(i+1)>X(j) O12(i,j)=1;endendif X(j)=b1O12(m,j)=1;endendO11=sum(O12')'for j=1:N(k)for i=1:mif (yduandian(i)<=Y(j)&(yduandian(i+1)>Y(j) O22(i,j)=1;endendif Y(j)=b2O22(m,j)=1;endendO21=sum(O22')'Person=sum(

5、O11-N(k).*p).2)./(N(k).*p); %person统计量的值 P2=sum(O21-N(k).*p).2)./(N(k).*p);PPerson=PPerson;Person;PP2=PP2;P2;g1=;q1=O12-P;for i=1:mt=1/N(k)*sum(q1(i,:);g1=g1;t;ends1=zeros(m,m);for i=1:N(k)s1=s1+q1(:,i)*q1(:,i)'endg2=;q2=O22-P;for i=1:mt=1/N(k)*sum(q2(i,:);g2=g2;t;ends2=zeros(m,m);for i=1:N(k)s2

6、=s2+q2(:,i)*q2(:,i)'ends1=1/N(k)*s1-g1*g1's2=1/N(k)*s2-g2*g2'LE=N(k)*g1'*(pinv(s1)*g1;PLE=PLE;LE;L=N(k)*g2'*(pinv(s2)*g2;PQ=PQ;L;%重新分割区间m=4,k1=2,b1=2 分割数为fennum=m!/(b!)%根据频数合% 降序排列得到的组pailie,ind1=sort(O11'); %O11就是频数，对其升序排列I=eye(m,m);I=rot90(I);pailie=pailie*I; %得到降序排列ind1=in

7、d1*I; %ind即为索引，就是对应原来的位置indexfyueshu=zeros(1,m-2); %fyueshu为可能的约数j=1;for i=2:m-1if m/i=fix(m/i)fyueshu(j)=m/i;j=j+1;endendyueshunum=j-1; %约数的数目K=zeros(1,yueshunum); %K为实际得到的约数B=zeros(1,yueshunum); %B为对应的m/KK=fyueshu(1:yueshunum);for i=1:yueshunumB(i)=m/K(i);endfengeX=cell(2,yueshunum); %创建一个元胞数组，用来存

8、放K所对应的aX(i)就是前面的O11，和O12Pxin=cell(1,yueshunum);Personxin=zeros(1,yueshunum);g1xin=cell(1,yueshunum);q1xin=cell(1,yueshunum);s1xin=cell(1,yueshunum);LExin=zeros(1,yueshunum);P2xin=zeros(1,yueshunum);g2xin=cell(1,yueshunum);q2xin=cell(1,yueshunum);s2xin=cell(1,yueshunum);Lxin=zeros(1,yueshunum);for ys

9、n=1:yueshunum %ysn就是第ysn个约数 这是一个大循环，求每一种组合的统计量 fengeX1,ysn=zeros(K(ysn),N(k); %按降序顺序重新组合区间，for j=1:N(k)for i=1:B(ysn):mfor ii=1:B(ysn)if(O12(ind1(i+ii-1),j)=1)fengeX1,ysn(i+B(ysn)-1)/B(ysn),j)=1;endendendendfengeX2,ysn=sum(fengeX1,ysn')'%重新计算区间概率生成的矩阵% %2个区间组合%Pxinysn=zeros(K(ysn),N(k); %Pxi

10、n1为k1下备择假设的各个区间的理论概率%Pxin2为改变分布后理论概率生成的矩阵for i=1:K(ysn)for j=1:B(ysn)Pxinysn(i,1)=Pxinysn(i,1)+P(ind1(i-1)*B(ysn)+j),1);endendfor i=2:N(k)Pxinysn(:,i)=Pxinysn(:,1);endPersonxin(ysn)=sum(fengeX2,ysn-N(k).*Pxinysn(:,1).2)./(N(k).*Pxinysn(:,1);%计算LEg1xinysn=; %q1xinysn=fengeX1,ysn-Pxinysn;%for i=1:K(ys

11、n)txin=1/N(k)*sum(q1xinysn(i,:);g1xinysn=g1xinysn;txin;end%s1xinysn=zeros(K(ysn),K(ysn);for i=1:N(k)s1xinysn=s1xinysn+q1xinysn(:,i)*q1xinysn(:,i)'ends1xinysn=1/N(k)*s1xinysn-g1xinysn*g1xinysn'%LExin(ysn)=N(k)*g1xinysn'*(pinv(s1xinysn)*g1xinysn;end%对O22进行重新组合%pailie1,ind2=sort(O21');

12、%O11就是频数，对其升序排列pailie1=pailie1*I; %得到降序排列ind2=ind2*I; %ind即为索引，就是对应原来的位置indexfengeY=cell(2,yueshunum);for ysn=1:yueshunum %ysn就是第ysn个约数 这是一个大循环，求每一种组合的统计量 fengeY1,ysn=zeros(K(ysn),N(k); %按降序顺序重新组合区间，for j=1:N(k)for i=1:B(ysn):mfor ii=1:B(ysn)if(O22(ind2(i+ii-1),j)=1)fengeY1,ysn(i+B(ysn)-1)/B(ysn),j)

13、=1;endendendendfengeY2,ysn=sum(fengeY1,ysn')'P2xin=sum(fengeY2,ysn-N(k).*Pxinysn(:,1).2)./(N(k).*Pxinysn(:,1);%计算LEg2xinysn=; %q2xinysn=fengeY1,ysn-Pxinysn;%for i=1:K(ysn)txin=1/N(k)*sum(q2xinysn(i,:);g2xinysn=g2xinysn;txin;ends2xinysn=zeros(K(ysn),K(ysn);for i=1:N(k)s2xinysn=s2xinysn+q2xiny

14、sn(:,i)*q2xinysn(:,i)'ends2xinysn=1/N(k)*s2xinysn-g2xinysn*g2xinysn'%Lxin=N(k)*g2xinysn'*(pinv(s2xinysn)*g2xinysn;endPersonX=max(Personxin);P2xinY=max(P2xin);PPersonxin=PPersonxin;PersonX;PP2xin=PP2xin;P2xinY;LExinX=max(LExin);PLExin=PLExin;LExinX;LxinY=max(Lxin);PQxin=PQxin;LxinY;h1,p,k

15、sstat1,cv=kstest(X,X tcdf(X,4),0.05); %计算KSh2,p,ksstat2,cv=kstest(X,X normcdf(X,0,1),0.05); T1=T1;ksstat1;T2=T2;ksstat2;t=t+h1;endxz=0.05;pos1,c=sort(PPerson);pos1xin,cxin=sort(PPersonxin);TJ1=pos1(1-xz)*cishu,1);TJ1xin=pos1xin(1-xz)*cishu,1);gx1=length(find(PP2>TJ1)/cishu;gx1xin=length(find(PP2xin>TJ1xin)/cishu;GX1=GX1;gx1; %Person的功效GX4=GX4;gx1xin; %pos2,c=sort(PLE);pos2xin,cxin=sort(PLExin);TJ2=pos2(ceil(1-xz)*cishu),1);TJ2xin=pos2xin(ceil(1-xz)*cishu),1);g