clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%simulation settings (model 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=200;p2=1000;
p=p2+4;
p1=p2;
nu=50;
rep=500;
m=2000;
for i=1:n
e(i)=1;
end
for i=1:4
ee(i)=1;
end
for i=1:p
ID(i,i)=1;
end
for i=1:p
sigma(i,i)=1;
end
for i=1:p
for j=1:p
sigma2(i,j)=0;
end
end
for i=5:50
for j=51:p
sigma2(i,j)=(0.5)^(j-50);
sigma2(j,i)=sigma2(i,j);
end
end

for i=1:(nu-1)
for j=(i+1):nu
sigma(i,j)=(4*log(p1)/n)^(1/2);
sigma(j,i)=(4*log(p1)/n)^(1/2);
end
end
for k=1:floor((p-nu)/10)
for i=(nu+(k-1)*10+1):(nu+k*10-1)
   for j=(i+1):(nu+k*10)
   sigma(i,j)=0.5;
   sigma(j,i)=sigma(i,j);
end
end
end
sigma1=sigma+sigma2;
ei=abs(min(eig(sigma1)));
sigma3=(sigma1+(ei+0.05)*ID)/(1+ei+0.05);
sig=sigma3^(1/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%code for the procedure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:m
t(j)=(j/m)*(2*log(p2)+1*log(log(p2)));
end
ttt=2*log(p2)+3*log(log(p2));
for j=1:m
nor(j)=1-chi2cdf(t(j),4);
end
for l=1:rep
x=randn(n,p)*sig;
bx=sum(x)/n;
x1=x-e'*bx;
y1=x1(:,1:4);
y2=x1(:,5:p);
z1=y2'*y1/n;
for k=1:(p-4)
A=y1.*(y2(:,k)*ee)-e'*z1(k,:);
TT(k)=n*z1(k,:)*inv(cov(A))*z1(k,:)';
end
for j=1:m
ren(j)=max(sum(TT>t(j)),1);
end

f0=nor-0.05*(ren)/(p-4);
if min(f0)<=0
at0(l)=min(find(f0<=0));
total0(l)=ren(at0(l));
ern0(l)=sum(TT((nu-4):(p-4))>t(at0(l)));
fdr0(l)=ern0(l)/max(1,total0(l));
power0(l)=(total0(l)-ern0(l))/(nu-4);
else
total0(l)=max(sum(TT>ttt),1);
ern0(l)=sum(TT((nu-4):(p-4))>ttt);
fdr0(l)=ern0(l)/max(1,total0(l));
power0(l)=(total0(l)-ern0(l))/(nu-4);
end

f1=nor-0.1*(ren)/(p-4);
if min(f1)<=0
at1(l)=min(find(f1<=0));
total1(l)=ren(at1(l));
ern1(l)=sum(TT((nu-4):(p-4))>t(at1(l)));
fdr1(l)=ern1(l)/max(1,total1(l));
power1(l)=(total1(l)-ern1(l))/(nu-4);
else
total1(l)=max(sum(TT>ttt),1);
ern1(l)=sum(TT((nu-4):(p-4))>ttt);
fdr1(l)=ern1(l)/max(1,total1(l));
power1(l)=(total1(l)-ern1(l))/(nu-4);
end

f2=nor-0.2*(ren)/(p-4);
if min(f2)<=0
at2(l)=min(find(f2<=0));
total2(l)=ren(at2(l));
ern2(l)=sum(TT((nu-4):(p-4))>t(at2(l)));
fdr2(l)=ern2(l)/max(1,total2(l));
power2(l)=(total2(l)-ern2(l))/(nu-4);
else
total2(l)=max(sum(TT>ttt),1);
ern2(l)=sum(TT((nu-4):(p-4))>ttt);
fdr2(l)=ern2(l)/max(1,total2(l));
power2(l)=(total2(l)-ern2(l))/(nu-4);
end

f3=nor-0.3*(ren)/(p-4);
if min(f3)<=0
at3(l)=min(find(f3<=0));
total3(l)=ren(at3(l));
ern3(l)=sum(TT((nu-4):(p-4))>t(at3(l)));
fdr3(l)=ern3(l)/max(1,total3(l));
power3(l)=(total3(l)-ern3(l))/(nu-4);
else
total3(l)=max(sum(TT>ttt),1);
ern3(l)=sum(TT((nu-4):(p-4))>ttt);
fdr3(l)=ern3(l)/max(1,total3(l));
power3(l)=(total3(l)-ern3(l))/(nu-4);
end
end
%%%%%%%%%%%%%%%%%%%
% FDR and Power
%%%%%%%%%%%%%%%%%%%
sum(fdr0)/rep
sum(power0)/rep
sum(fdr1)/rep
po=sum(power1)/rep
sum(fdr2)/rep
sum(power2)/rep
sum(fdr3)/rep
sum(power3)/rep