t1=2.4390e-02; t2=2.4994e-02; t3=1.2748e-01; t4=9.3777e-02;

% ********* hc mine ************* 5: filter 4 s=2, 6: filter 4 s=3
x=[[	2.9362e-02   3.0444e-02   1.1499e-01   9.0615e-02   4.1318e-04 ...
	4.3212e-04   1.3099e-03   1.2893e-04]
[	2.7541e-02   2.8362e-02   1.1714e-01   9.1951e-02   4.1773e-04 ...
	 4.3110e-04  1.4755e-03   1.2963e-04]
[	2.5723e-02   2.6320e-02   1.1991e-01   9.3349e-02   4.1474e-04 ...
	4.7779e-04   1.5966e-03   1.3156e-04]
[	2.5170e-02   2.5734e-02   1.2017e-01   9.3810e-02   4.3111e-04 ... 
	4.8084e-04   1.8418e-03   1.3558e-04]
[	3.0604e-02   3.2285e-02   1.1663e-01   8.8617e-02   3.3056e-04 ...
	3.4482e-04   1.3194e-03   1.1490e-04]
[	3.4323e-02   3.6710e-02   1.1391e-01   8.5190e-02   2.9875e-04 ...
	3.1011e-04   1.0840e-03   1.0661e-04]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[4 2 6];
figure(1),p1=plot(x(sel,1)-100,x(sel,5),'s'); hold on
        plot(x(sel,1)-100,x(sel,5)); hold off; plerr(x(sel,1)-100,x(sel,5));
figure(2),p2=plot(x(sel,2)-100,x(sel,6),'s'); hold on
        plot(x(sel,2)-100,x(sel,6)); hold off; plerr(x(sel,2)-100,x(sel,6));
%subplot(223),p3=plot(abs(x(sel,3)-t3),x(sel,7),'s'); hold on
%        plot(abs(x(sel,3)-t3),x(sel,7)); hold off;
%	plerr(abs(x(sel,3)-t3),x(sel,7));
%subplot(224),p4=plot(abs(x(sel,4)-t4),x(sel,8),'s'); hold on
%        plot(abs(x(sel,4)-t4),x(sel,8)); hold off; 
%	plerr(abs(x(sel,4)-t4),x(sel,8));

% ********* huber delta=0.004 *************
x=[[	4.5140e-02   4.5854e-02   1.3875e-01   1.0416e-01   1.4987e-03 ...
	1.4192e-03   3.8931e-03   4.7973e-04]
[	2.6890e-02   2.7434e-02   1.2623e-01   9.4994e-02   8.2589e-04 ...
	6.7251e-04   2.4495e-03   1.9957e-04]
[	2.5046e-02   2.5499e-02   1.2383e-01   9.4477e-02   5.9997e-04 ...
	6.6327e-04   2.4047e-03   1.9584e-04]
[	2.5666e-02   2.6197e-02   1.2065e-01   9.4169e-02   5.6715e-04 ...
	6.1229e-04   2.2241e-03   1.9653e-04]
[	2.7088e-02   2.7792e-02   1.1662e-01   9.3590e-02   5.4059e-04 ...
	5.5952e-04   1.7389e-03   1.9456e-04]
[	2.9484e-02   3.0469e-02   1.1302e-01   9.2445e-02   4.9760e-04 ...
	5.1271e-04   1.1960e-03   1.8387e-04]
[	3.7291e-02   3.9580e-02   1.0759e-01   8.7959e-02   3.4659e-04 ...
	3.7091e-04   5.5348e-04   1.3405e-04]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[2:7];
figure(1),hold on; q1=plot(x(sel,1)-100,x(sel,5),'o');
	plot(x(sel,1)-100,x(sel,5),'--'); hold off; plerr(x(sel,1)-100,x(sel,5))
figure(2),hold on; q2=plot(x(sel,2)-100,x(sel,6),'o');
	plot(x(sel,2)-100,x(sel,6),'--'); hold off; plerr(x(sel,2)-100,x(sel,6))
%subplot(223),hold on; q3=plot(abs(x(sel,3)-t3),x(sel,7),'o');
%	plot(abs(x(sel,3)-t3),x(sel,7),'--'); hold off
%	plerr(abs(x(sel,3)-t3),x(sel,7));
%subplot(224),hold on; q4=plot(abs(x(sel,4)-t4),x(sel,8),'o');
%	plot(abs(x(sel,4)-t4),x(sel,8),'--'); hold off
%	plerr(abs(x(sel,4)-t4),x(sel,8));

% ********* huber delta=0.002 *************
x=[[	3.3213e-02   3.3873e-02   1.2975e-01   9.8175e-02   1.1417e-03 ...
	9.7699e-04   2.6090e-03   3.0513e-04]
[	2.6815e-02   2.7374e-02   1.2631e-01   9.5140e-02   8.3241e-04 ...
	6.7748e-04   2.4163e-03   2.1555e-04]
[	2.5414e-02   2.5946e-02   1.2082e-01   9.4287e-02   5.6203e-04 ...
	6.2437e-04   2.2868e-03   1.9624e-04]
[	2.6526e-02   2.7201e-02   1.1591e-01   9.3859e-02   5.4315e-04 ...
	5.8040e-04   2.1112e-03   2.0165e-04]
[	3.2233e-02   3.3664e-02   1.0548e-01   9.0996e-02   4.4494e-04 ...
	4.7211e-04   6.8214e-04   1.6987e-04]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[1:5];
figure(1),hold on; r1=plot(x(sel,1)-100,x(sel,5),'+');
        plot(x(sel,1)-100,x(sel,5),'-.'); hold off; %plerr(x(sel,1)-100,x(sel,5))
figure(2),hold on; r2=plot(x(sel,2)-100,x(sel,6),'+');
        plot(x(sel,2)-100,x(sel,6),'-.'); hold off; %plerr(x(sel,2)-100,x(sel,6))
%subplot(223),hold on; r3=plot(abs(x(sel,3)-t3),x(sel,7),'+');
%        plot(abs(x(sel,3)-t3),x(sel,7),'-.'); hold off
%	plerr(abs(x(sel,3)-t3),x(sel,7))
%subplot(224),hold on; r4=plot(abs(x(sel,4)-t4),x(sel,8),'+');
%        plot(abs(x(sel,4)-t4),x(sel,8),'-.'); hold off
%	plerr(abs(x(sel,4)-t4),x(sel,8))

% ********* huber delta=0.008 *************
x=[[	2.5159e-02   2.5603e-02   1.2402e-01   9.4432e-02   6.1318e-04 ...
	6.7017e-04   2.3600e-03   1.9883e-04]
[	2.5983e-02   2.6534e-02   1.2142e-01   9.3989e-02   5.7065e-04 ...
	6.0428e-04   2.1317e-03   1.9507e-04]
[	3.3998e-02   3.5631e-02   1.1342e-01   8.9901e-02   3.9906e-04 ...
	4.0144e-04   1.1666e-03   1.6310e-04]
[	4.3865e-02   4.7084e-02   1.0556e-01   8.2987e-02   2.7271e-04 ...
	2.7605e-04   6.3454e-04   1.2614e-04]
[	5.8700e-02   6.1714e-02   8.9655e-02   6.8060e-02   1.5362e-04 ...
	1.5584e-04   1.7743e-04   6.0570e-05]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[1:5];
figure(1),hold on; s1=plot(x(sel,1)-100,x(sel,5),'x');
        plot(x(sel,1)-100,x(sel,5),':'); hold off; %plerr(x(sel,1)-100,x(sel,5))
figure(2),hold on; s2=plot(x(sel,2)-100,x(sel,6),'x');
        plot(x(sel,2)-100,x(sel,6),':'); hold off; %plerr(x(sel,2)-100,x(sel,6))
%subplot(223),hold on; s3=plot(abs(x(sel,3)-t3),x(sel,7),'x');
%        plot(abs(x(sel,3)-t3),x(sel,7),':'); hold off
%        plerr(abs(x(sel,3)-t3),x(sel,7))
%subplot(224),hold on; s4=plot(abs(x(sel,4)-t4),x(sel,8),'x');
%        plot(abs(x(sel,4)-t4),x(sel,8),':'); hold off
%        plerr(abs(x(sel,4)-t4),x(sel,8))

figure(1),legend([q1 r1 s1 p1],'Huber:  \delta=0.004',...
        'Huber:  \delta=0.002','Huber:  \delta=0.008','Proposed method');
axis([0 80 1 5])
%xlabel('normalized bias (%)'); ylabel('normalized standard deviation (%)');
figure(2),legend([q2 r2 s2 p2],'Huber:  \delta=0.004',...
        'Huber:  \delta=0.002','Huber:  \delta=0.008','Proposed method'); 
axis([0 80 1 5])
%xlabel('normalized bias (%)'); ylabel('normalized standard deviation (%)');
%subplot(223),legend([p3 q3 r3 s1],'Proposed method','Huber:  \delta=0.004',...
%        'Huber:  \delta=0.002','Huber:  \delta=0.008');
%subplot(224),legend([p4 q4 r4 s1],'Proposed method','Huber:  \delta=0.004',...
%        'Huber:  \delta=0.002','Huber:  \delta=0.008');

% ********* lc mine ************* 5: filter 1 s=1, 6: filter 1 s=2
%	7: filter 1 s=4
x=[[	2.8997e-02   3.0366e-02   1.1713e-01   9.1178e-02   6.3987e-04 ...
	7.6710e-04   2.4986e-03   2.2285e-04]
[	2.8865e-02   3.0277e-02   1.1694e-01   9.1224e-02   6.4969e-04 ...
	8.1283e-04   2.7745e-03   2.1529e-04]
[	2.8877e-02   3.0395e-02   1.1673e-01   9.1205e-02   6.9958e-04 ...
	8.8605e-04   3.0168e-03   2.1496e-04]
[	2.7545e-02   2.8878e-02   1.1785e-01   9.2329e-02   7.7841e-04 ...
	9.9055e-04   3.2918e-03   2.3216e-04]
[	3.0030e-02   3.1634e-02   1.1645e-01   9.0015e-02   6.1790e-04 ...
	7.2806e-04   2.4160e-03   2.1693e-04]
[	3.2480e-02   3.4635e-02   1.1475e-01   8.7492e-02   5.8475e-04 ...
	6.5753e-04   2.2201e-03   2.0611e-04]
[	3.8946e-02   4.2221e-02   1.0983e-01   8.1270e-02   5.1498e-04 ...
	5.3217e-04   1.7254e-03   1.8445e-04]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[1 5 6];
figure(3),p1=plot(x(sel,1)-100,x(sel,5),'s'); hold on
        plot(x(sel,1)-100,x(sel,5)); hold off; plerr(x(sel,1)-100,x(sel,5))
figure(4),p2=plot(x(sel,2)-100,x(sel,6),'s'); hold on
        plot(x(sel,2)-100,x(sel,6)); hold off; plerr(x(sel,2)-100,x(sel,6))
%subplot(223),p3=plot(abs(x(sel,3)-t3),x(sel,7)-2.9e-4,'s'); hold on
%        plot(abs(x(sel,3)-t3),x(sel,7)-2.6e-4); hold off
%	plerr(abs(x(sel,3)-t3),x(sel,7)-2.6e-4)
%subplot(224),p4=plot(abs(x(sel,4)-t4),x(sel,8),'s'); hold on
%        plot(abs(x(sel,4)-t4),x(sel,8)); hold off
%	plerr(abs(x(sel,4)-t4),x(sel,8))

% ********* huber delta=0.004 *************
x=[[	6.0159e-02   6.1436e-02   1.5352e-01   1.1846e-01   2.8677e-03 ...
	3.1587e-03   7.6374e-03   8.6796e-04]
[	2.6922e-02   2.7578e-02   1.2302e-01   9.4014e-02   9.2443e-04 ...
	1.0940e-03   3.3598e-03   3.2921e-04]
[	2.6827e-02   2.7649e-02   1.1857e-01   9.3393e-02   7.8877e-04 ...
	9.4385e-04   2.7944e-03   2.9861e-04]
[	2.8834e-02   2.9832e-02   1.1455e-01   9.2522e-02   7.5390e-04 ...
	8.5511e-04   1.9315e-03   2.7846e-04]
[	3.1864e-02   3.3290e-02   1.1086e-01   9.0958e-02   7.0709e-04 ...
	8.0397e-04   1.4905e-03   2.6513e-04]
[	3.6380e-02   3.8614e-02   1.0789e-01   8.8395e-02   6.2722e-04 ...
	7.0940e-04   1.1342e-03   2.3775e-04]
[	4.6761e-02   5.0766e-02   1.0527e-01   8.1663e-02   3.6911e-04 ...
	4.2781e-04   4.4311e-04   1.4414e-04]
[	5.2659e-02   5.6699e-02   1.0185e-01   7.4692e-02   1.1541e-04 ...
	1.3456e-04   9.9207e-05   5.9333e-05]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[2:7];
figure(3),hold on; q1=plot(x(sel,1)-100,x(sel,5),'o');
        plot(x(sel,1)-100,x(sel,5),'--'); hold off; plerr(x(sel,1)-100,x(sel,5))
figure(4),hold on; q2=plot(x(sel,2)-100,x(sel,6),'o');
        plot(x(sel,2)-100,x(sel,6),'--'); hold off; plerr(x(sel,2)-100,x(sel,6))	
%subplot(223),hold on; q3=plot(abs(x(sel,3)-t3),x(sel,7),'o');
%        plot(abs(x(sel,3)-t3),x(sel,7),'--'); hold off
%	plerr(abs(x(sel,3)-t3),x(sel,7))
%subplot(224),hold on; q4=plot(abs(x(sel,4)-t4),x(sel,8),'o');
%        plot(abs(x(sel,4)-t4),x(sel,8),'--'); hold off
%	plerr(abs(x(sel,4)-t4),x(sel,8))

% ********* huber delta=0.002 *************
x=[[	3.2663e-02   3.3415e-02   1.2677e-01   9.6873e-02   1.3443e-03 ...
	1.5066e-03   3.9583e-03   4.0885e-04]
[	2.6807e-02   2.7471e-02   1.2301e-01   9.4045e-02   9.3481e-04 ...
	1.0819e-03   3.3953e-03   3.2987e-04]
[	2.8249e-02   2.9143e-02   1.1200e-01   9.2856e-02   7.7336e-04 ...
	8.8753e-04   2.3134e-03   2.8368e-04]
[	3.1215e-02   3.2507e-02   1.0665e-01   9.1361e-02   7.2662e-04 ...
	8.4508e-04   1.2894e-03   2.6894e-04]
[	4.3792e-02   4.7506e-02   9.7810e-02   8.4315e-02   6.4300e-04 ...
	6.8716e-04   7.3977e-04   2.2103e-04]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[1:5];
figure(3),hold on; r1=plot(x(sel,1)-100,x(sel,5),'+');
        plot(x(sel,1)-100,x(sel,5),'-.'); hold off; %plerr(x(sel,1)-100,x(sel,5))
figure(4),hold on; r2=plot(x(sel,2)-100,x(sel,6),'+');
        plot(x(sel,2)-100,x(sel,6),'-.'); hold off; %plerr(x(sel,2)-100,x(sel,6))
%subplot(223),hold on; r3=plot(abs(x(sel,3)-t3),x(sel,7),'+');
%        plot(abs(x(sel,3)-t3),x(sel,7),'-.'); hold off
%	plerr(abs(x(sel,3)-t3),x(sel,7))
%subplot(224),hold on; r4=plot(abs(x(sel,4)-t4),x(sel,8),'+');
%        plot(abs(x(sel,4)-t4),x(sel,8),'-.'); hold off
%	plerr(abs(x(sel,4)-t4),x(sel,8))

% ********* huber delta=0.008 *************
x=[[	2.7053e-02   2.7865e-02   1.1985e-01   9.3241e-02   7.8497e-04 ...
	9.3605e-04   2.7914e-03   3.0057e-04]
[	2.9373e-02   3.0384e-02   1.1696e-01   9.2177e-02   7.5715e-04 ...
	8.3748e-04   2.2303e-03   2.8270e-04]
[	3.7100e-02   3.9204e-02   1.1085e-01   8.7800e-02   5.9820e-04 ...
	6.3890e-04   1.4920e-03   2.3451e-04]
[	4.2285e-02   4.5286e-02   1.0672e-01   8.4129e-02   5.2043e-04 ...
	5.5363e-04   1.1676e-03   2.0750e-04]
[	5.4811e-02   5.8575e-02   9.5393e-02   7.3729e-02   4.1500e-04 ...
	4.5143e-04   6.4406e-04   1.4590e-04]
[	5.9453e-02   6.1496e-02   8.1321e-02   6.1429e-02   1.5634e-04 ...
	1.6787e-04   1.6228e-04   6.0026e-05]];
x(:,1)=x(:,1)/t1*100; x(:,5)=x(:,5)/t1*100;
x(:,2)=x(:,2)/t2*100; x(:,6)=x(:,6)/t2*100;

sel=[1:5];
figure(3),hold on; s1=plot(x(sel,1)-100,x(sel,5),'x');
        plot(x(sel,1)-100,x(sel,5),':'); hold off; %plerr(x(sel,1)-100,x(sel,5))
figure(4),hold on; s2=plot(x(sel,2)-100,x(sel,6),'x');
        plot(x(sel,2)-100,x(sel,6),':'); hold off; %plerr(x(sel,2)-100,x(sel,6))
%subplot(223),hold on; s3=plot(abs(x(sel,3)-t3),x(sel,7),'x');
%        plot(abs(x(sel,3)-t3),x(sel,7),':'); hold off
%        plerr(abs(x(sel,3)-t3),x(sel,7))
%subplot(224),hold on; s4=plot(abs(x(sel,4)-t4),x(sel,8),'x');
%        plot(abs(x(sel,4)-t4),x(sel,8),':'); hold off
%        plerr(abs(x(sel,4)-t4),x(sel,8))

figure(3),legend([q1 r1 s1 p1],'Huber:  \delta=0.004',...
        'Huber:  \delta=0.002','Huber:  \delta=0.008','Proposed method');
axis([0 80 1 6])
%xlabel('normalized bias (%)'); ylabel('normalized standard deviation (%)');
figure(4),legend([q2 r2 s2 p2],'Huber:  \delta=0.004',...
        'Huber:  \delta=0.002','Huber:  \delta=0.008','Proposed method');
axis([0 80 1 6])
%xlabel('normalized bias (%)'); ylabel('normalized standard deviation (%)');
%subplot(223),legend([p3 q3 r3 s1],'Proposed method','Huber:  \delta=0.004',...
%        'Huber:  \delta=0.002','Huber:  \delta=0.008');
%subplot(224),legend([p4 q4 r4 s1],'Proposed method','Huber:  \delta=0.004',...
%        'Huber:  \delta=0.002','Huber:  \delta=0.008');

% **************** EM ************8
% hc mine 4
%1.1421e-02 7.8574e-02 1.8858e-02 9.6443e-02 1.3193e-01 1.4834e-01 1.6121e-01
% hc huber 12
%4.0956e-02 9.1109e-02 5.9379e-02 1.6777e-01 9.9986e-02 1.3755e-01 1.3848e-01
% hc huber 13
%3.1873e-02 6.7723e-02 4.7220e-02 1.1004e-01 8.2457e-02 1.0840e-01 1.2213e-01
% hc huber 11
%5.5391e-02 1.3513e-01 7.7144e-02 2.4624e-01 1.4602e-01 2.0037e-01 1.9050e-01
% hc huber 11.3
%4.7688e-02 1.1008e-01 6.7367e-02 2.0649e-01 1.1748e-01 1.6356e-01 1.5455e-01
% lc mine 4
%1.7243e-02 8.3898e-02 4.0734e-02 1.6743e-01 1.6209e-01 2.1551e-01 1.6730e-01
% lc mine 3
%1.7031e-02 8.7804e-02 3.9059e-02 1.7657e-01 1.5328e-01 2.0737e-01 1.6935e-01
% lc huber 12
%4.2051e-02 1.2033e-01 7.4016e-02 1.1490e-01 1.3797e-01 1.4124e-01 1.5071e-01
% lc huber 11
%5.4864e-02 1.5617e-01 9.5103e-02 1.6415e-01 1.6126e-01 1.8206e-01 1.7500e-01
% lc huber 11.3
%5.2339e-02 1.4792e-01 9.2562e-02 1.4827e-01 1.5418e-01 1.7227e-01 1.6959e-01

% **************** EM ************ good
% hc huber
%4.7788e-02   1.1040e-01   6.7470e-02   2.0675e-01   1.1799e-01   1.6418e-01
%   1.5532e-01
% hc mine
%1.1446e-02   8.0830e-02   1.8888e-02   1.0042e-01   1.4384e-01   1.7053e-01
%   1.8063e-01
% lc huber
%5.0283e-02   1.4068e-01   8.9725e-02   1.4197e-01   1.3999e-01   1.5725e-01
%   1.5309e-01
% lc mine
%1.6888e-02   8.5284e-02   3.8653e-02   1.7868e-01   1.3420e-01   1.7967e-01
%   1.4580e-01
















