function [risk,eps1,eps2,inter1]=bayeserr(p1,m1,m2,c1,c2)
% BAYESERR computes the Bayesian risk for optimal classifier.
% [risk,eps1,eps2,inter1]=bayeserr(p1,m1,m2,c1,c2)
%
% BAYESERR computes a value of the Bayesian risk for the 0/1-cost
%  function and an optimal Bayesian classifier minimizing the 
%  risk. The binary clasification problem is supposed, where 
%  the conditional probabilities are one-dimensional Gaussians 
%  p(x|k=1)=N(m1,c1) and p(x|k=2)=N(m2,c2). 
%
% Input:
%  p1 [0<=p1<=1]      Apriori probability of the first class.
%  m1 [real]          Mean value of the Gaussian p(x|k=1)=N(m1,c1).
%  m2 [real]          Mean value of the Gaussian p(x|k=2)=N(m2,c2).
%  c1 [positive real] Variance of the Gaussian p(x|k=1)=N(m1,c1).
%  c2 [positive real] Variance of the Gaussian p(x|k=2)=N(m2,c2).
%
% Output:
%  risk [positive real] Bayesian risk for an optimal classifier.
%  eps1 [positive real] Integral of p(x|k=1) over x in L2, where
%    L2 is an arrea where x is classified to the 2nd class.
%  eps2 [positive real] Integral of p(x|k=1) over x in L1, where
%    L1 is an arrea where x is classified to the 1nd class.
%  inter1 [1x2] or [1x4] One or two intervals describing L1.
%
% See also BAYES.

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Modifications:
% 27-Oct-2001, V.Franc


% p2 is uniquely determined since p1+p2=1.
p2=1-p1;

% finds out decision function which is generaly quadratic
[a,b,c]=bayesnd(p1,p2,m1,m2,c1,c2);


% Splits X into X1 and X2 according to the computed quadratic 
% discriminat function ax^2 + bx + c = 0 and computes 
% eps1 and eps2.

if a==0,
   % The decision function is linear, i.e. in 1D it is a
   % single threshold.
   th=-c/b;
   inter1=[th,inf];
   
   % gets label for the interval (th,inf)
   class=classify(th+1,p1,p2,m1,m2,c1,c2);
   
   if abs(c)==inf,
     risk=0;
     if class==1,,
       eps1=0;
       eps2=1;
       inter1=[-inf,inf];
       return;
     else
      eps1=1;
      eps2=0;
      inter1=[];
     end
   end   
      
   eps1=1-erfc2(th,m1,sqrt(c1));
   eps2=erfc2(th,m2,sqrt(c2));  
         
   if class==2,
      % swap eps1 and eps2
      tmp=eps2;
      eps2=eps1;
      eps1=tmp;
   end
   
else
   % The decision function is quadratic, i.e. in 2d
   % there exis two thresholds which determine three intervals.
   
   D=b^2-4*a*c;
   if D > 0,
      th1=(-b-sqrt(D))/(2*a);
      th2=(-b+sqrt(D))/(2*a);
      
      if th1 > th2,
         tmp=th1;
         th1=th2;
         th2=tmp;         
      end;
      
      % finds out label for the interval [th1,th2].
      class=classify((th2+th1)/2,p1,p2,m1,m2,c1,c2);
      
      if class==2
         % integral eps2 = int_th2^inf + int_{-inf}^th1
         eps2 = 1 + erfc2(th2,m2,sqrt(c2)) - erfc2(th1,m2,sqrt(c2));          
         % integral eps1= int_th1^th2
         eps1 = erfc2(th1,m1,sqrt(c1)) - erfc2(th2,m1,sqrt(c1));         
         
         inter1=[-inf,th1,th2,inf];
      else
         % integral eps1 = int_th2^inf + int_{-inf}^th1
         eps1 = 1 + erfc2(th2,m1,sqrt(c1)) - erfc2(th1,m1,sqrt(c1));  
         % integral eps2= int_th1^th2
         eps2=erfc2(th1,m2,sqrt(c2))-erfc2(th2,m2,sqrt(c2));
         
         inter1=[th1,th2];
      end
      
   else
      % finds out label for the interval [-inf,inf].
      class=classify(0,p1,p2,m1,m2,c1,c2);
      
      if class == 1,
         eps1=0;
         eps2=1;
         inter1=[-inf,inf];
      else
         eps1=1;
         eps2=0;
         inter1=[];
      end
      risk=0;
      return;
   end
end

% computes the Bayesian risk 
risk = p1*( eps1 - eps2 ) + eps2;

return;

%-----------------------------------------------
function class = classify(x,p1,p2,m1,m2,c1,c2)
% finds out to which class the given x belongs

if p1==1,     % only the 1st class is possible
   class=1;
elseif p2==1, % only the 2nd class is possible
   class=2;
elseif normald(x,m1,c1)*p1 > normald(x,m2,c2)*p2,
   class=1;
else
   class=2;
end

return;