clc; clear; load data.mat Y = data; nx = size(Y,1); ny = size(Y,2); nt = size(Y,3); p=nx*ny; Y=reshape(Y,p,nt); Y=(Y-mean(Y,2))./std(Y,0,2);%standardization Y0=Y; Y0(:,187:202)=[]; Y0(:,200:252)=[];%delete the probable anomaly data ntt=size(Y0,2); [coeff,score,latent,tsquared,explained,mu]=pca(Y0','NumComponents',20);%use pca to decide sigma of noise and B00 fit=score*coeff'+repmat(mu,ntt,1); noise=Y0'-fit; sq=(noise-mean(noise)).^2; sigma=sqrt(sum(sq(:))/(p*ntt));%standard deviation of noise B00=coeff;% normal bases cov_theta=cov(score); sigma0=mean(diag(cov_theta)); skx = round(nx/15); B= bsplineBasis(ny,skx,3); B=kron(B,B); for i=1:292 B((232*i+1):(232*(i-1)+292),:)=[]; end B(:,[15,16,31,32,47,48,63,64,79,80,95,96,111,112,127,128,143,144,159,160,175,176,191,192,207,208,223,224,239,240,255,256])=[]; nq = size(B,2); w=0.1*ones(nq,1); sigma9=10*ones(nq,1); nrep=100; m = 400; niter = 30; controlimit=20; lambda0=0.9; %ARL0=1100; T=3000; parpool('local',4); parfor ss=1:nrep c1=1:100; myrand=randsample(c1,T,true); YY=Y(:,myrand)-repmat(mu,T,1)'; ss [index_ini,index8,theta8,theta80,Test,runlength] = thom66(YY',B00,B,m,sigma,niter,w,sigma9,lambda0,controlimit,sigma0); TTest(ss,:)=Test; RRunlength(ss)=runlength; end ARL= mean(RRunlength); sdrl=std(RRunlength); delete(gcp('nocreate'))