Code:
u1= load('UNITED.txt');
u2= load('STATES.txt');
u3= load('HONG.txt');
u4= load('kONG.txt');
m1= sum(abs(u1.^2)); u1= u1./sqrt(m1);
m2= sum(abs(u2.^2)); u2= u2./sqrt(m2);
m3= sum(abs(u3.^2)); u3= u3./sqrt(m3);
m4= sum(abs(u4.^2)); u4= u4./sqrt(m4);
k= [20, 50, 100, 200, 300, 400];
a= sum(u1.*u2);
aHatUS0= zeros(100,length(k)); aHatUS= zeros(100,length(k));
aHatHK0= zeros(100,length(k)); aHatHK= zeros(100,length(k));
for t= 1:100
for i= 1:length(k)
R= randn(length(u1),k(i));
v1= R'*u1; v2= R'*u2; v3= R'*u3; v4= R'*u4;
aHatUS0(t,i)= (1/k(i)).*sum(v1.*v2);
aHatHK0(t,i)= (1/k(i)).*sum(v3.*v4);
xn1= aHatUS0; xn2= aHatHK0; n= length(v1);
for j=1:10
a1= xn1; a2= xn2;
f1= n*a1.^3 - a1.^2 .* sum (v1.*v2) + a1.*(-n+sum(v1.^2+v2.^2)) - sum (v1.*v2);
fp1= 3*n*a1.^2 - 2.*a1.* sum (v1.*v2) -n + sum(v1.^2+v2.^2);
f2= n*a2.^3 - a2.^2 .* sum (v3.*v4) + a2.*(-n+sum(v3.^2+v4.^2)) - sum (v3.*v4);
fp2= 3*n*a2.^2 - 2.*a2.* sum (v3.*v4) -n + sum(v3.^2+v4.^2);
xn1= a1 - f1./fp1; xn2= a2 - f2./fp2;
end
aHatUS(t,i)=xn1(1);
aHatHK(t,i)=xn2(1);
end
end
aHatUS0Var= (m1*m2+a^2)*k.^(-1); aHatUSVar= (1-a^2)^2./((1+a^2)*k);
aHatHK0Var= (m3*m4+a^2)*k.^(-1); aHatHKVar= (1-a^2)^2./((1+a^2)*k);
mSEaHatUS0= aHatUS0Var + (mean(aHatUS0-a)).^2;
mSEaHatHK0= aHatHK0Var + (mean(aHatHK0-a)).^2;
mSEaHatUS= aHatUSVar + (mean(aHatUS-a)).^2;
mSEaHatHK= aHatHKVar + (mean(aHatHK-a)).^2;
figure; hold on; grid on;
plot(k,mSEaHatUS0,'g');
plot(k,aHatUS0Var,'b');
plot(k,mSEaHatHK0,'r');
plot(k,aHatHK0Var,'y');
title('Figure 1: Empirical MSEs and Theoretical Variances')
xlabel('k'); ylabel('MSE & Variance');
h=legend('Empirical MSE of aUS0_{hat}','Theoretical Variance of aUS0_{hat}','Empirical MSE of aHK0_{hat}','Theoretical Variance of aHK0_{hat}',4);
set(h,'Location','NorthEast')
hold off;
figure; hold on; grid on;
plot(k,mSEaHatUS,'g');
plot(k,aHatUSVar,'b');
plot(k,mSEaHatHK,'r');
plot(k,aHatHKVar,'y');
title('Figure 1: Empirical MSEs and Theoretical Variances')
xlabel('k'); ylabel('MSE & Variance');
h=legend('Empirical MSE of aUS_{hat}','Theoretical Variance of aUS_{hat}','Empirical MSE of aHK_{hat}','Theoretical Variance of aHK_{hat}',4);
set(h,'Location','NorthEast')
hold off;