% Simulation for Bayes Decision
% m1=[1; 1];  m2=-m1;  S1=S2=I;
% p1=p2=0.5; n1=n2=50 for training and testing
% ideal error = \Phi(-\sqrt{2});
%
m1=[1;1]; m2=-m1; S1=eye(2); S2=eye(2); n1=50; n2=50;

V1=random('Normal',1,1,n1,1);  V2=random('Normal',1,1,n1,1);  Tx=[V1,V2];
W1=random('Normal',-1,1,n2,1); W2=random('Normal',-1,1,n2,1); Ty=[W1,W2];
X1=random('Normal',1,1,n1,1);  X2=random('Normal',1,1,n1,1);  X=[X1, X2];
Y1=random('Normal',-1,1,n2,1); Y2=random('Normal',-1,1,n2,1); Y=[Y1, Y2];
p1=n1/(n1+n2);  p2=n2/(n1+n2);
u1=mean(Tx)'; u2=mean(Ty)'; C1=cov(Tx); C2=cov(Ty)';
s1=sqrt(det(C1));s2=sqrt(det(C2));invC1=inv(C1);invC2=inv(C2);
errx=0;  erry=0;
for i=1:n1
  t=X(i,:);  x=t';
  t=Y(i,:);  y=t';
  q1=1/(2*pi*s1)*exp(-0.5*(x-u1)'*invC1*(x-u1))*p1;
  q2=1/(2*pi*s2)*exp(-0.5*(x-u2)'*invC2*(x-u2))*p2;
  if (q1<q2)
    errx=errx+1;
  end
  q1=1/(2*pi*s1)*exp(-0.5*(y-u1)'*invC1*(y-u1))*p1;
  q2=1/(2*pi*s2)*exp(-0.5*(y-u2)'*invC2*(y-u2))*p2;
  if (q1>q2)
    erry=erry+1;
  end
end
[errx, erry, errx+erry, n1+n2]