sc.m

score检验统计

clear;
clc;
tic;
randn('state',0);

for n=25:25:100;
    a=1
   for gamma=-0.9:0.1:0.9;
% delta=1.0e-003;

      for j=1:1000
          [x,y,m,epsilon,u,G]=data(n,gamma);
   %求score检验统计量
          C=eye(n)-(1/n)*ones(n);
          G1=C*G';
          H=inv(G1'*G1);
          SC(j)=u*G1*H*G1'*u'/2;
        if (abs(SC(j))>3.841)
            sc(j)=1;
        else
           sc(j)=0;
        end
      
      end
    
 psc(n/25,a)=sum(sc)/1000; 
 a=a+1;

%   psc=sum(sc)/1000;
    end
end