bayesboundary.m

內涵模糊理論與類神經網路的程式碼...提供初學者做研究參考

% ==========================================================
% 
%           Neural Networks A Classroom Approach
%                     Satish Kumar
%             Copyright Tata McGraw Hill, 2004
%
%     MATLAB code to plot Bayesian decision boundaries
%              Reference:Table 7.3;Page 235
%
% ==========================================================

n = 200;
dim = 2;
ncentres = 2;

priors = [0.5 0.5];
centres = [1 1; 2.5 2.5];
covars(:,:,1) = [0.6 -0.2; -0.2 0.6];
covars(:,:,2) = [0.2 -0.1; -0.1 0.2];

% Determine number to sample from each component
samplepriors = rand(1, n);

data = zeros(n, dim);	% Pre-allocate data array
label = zeros(n, 1);
cum_prior = 0;				% Cumulative sum of priors
total_samples = 0;		% Cumulative sum of number of sampled points

for i = 1:ncentres
  num_samples = sum(samplepriors >= cum_prior &  samplepriors < cum_prior + priors(i));
  covar = covars(:,:,i);
  d = size(covar, 1);
  [evec, eval] = eig(covar);
  coeffs = randn(num_samples, d)*sqrt(eval);
  samples(total_samples+1:total_samples+num_samples, :)=ones(num_samples,1)*centres(i,:)+coeffs*evec';
  label(total_samples+1:total_samples+num_samples) = i;
  cum_prior = cum_prior + priors(i);
  total_samples = total_samples + num_samples;
end

x0 = min(samples(:,1));
x1 = max(samples(:,1));
y0 = min(samples(:,2));
y1 = max(samples(:,2));

dx = x1-x0;
dy = y1-y0;
expand = 5/100;			% Add on 5 percent each way

x0 = x0 - dx*expand;
x1 = x1 + dx*expand;
y0 = y0 - dy*expand;
y1 = y1 + dy*expand;

resolution = 30;

step = dx/resolution;

xrange = [x0:step:x1];
yrange = [y0:step:y1];

% Grid points for probability evaluation
[X Y]=meshgrid([x0:step:x1],[y0:step:y1]);

figure(1);
hold on
x=[X(:) Y(:)];
nsamples = size(x, 1);

a = zeros(nsamples, ncentres);  % Preallocate matrix
normal = (2*pi)^(dim/2);

for i = 1:ncentres
   diffs = x - (ones(nsamples, 1) * centres(i,:));
% Use Cholesky decomposition of covariance matrix to speed computation
   c = chol(covars(:,:,i));
   temp = diffs/c;
   a(:,i) = exp(-0.5*sum(temp.*temp, 2))./(normal*prod(diag(c)));
end

px_j=a;

%Code inserted for contour plotting
P1x = reshape(a(:,1),size(X));
P2x = reshape(a(:,2),size(X));
contour(X,Y,P1x);
contour(X,Y,P2x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

px = reshape(px_j*(priors)',size(X));
post = (ones(nsamples, 1)*priors).*a;
s = sum(post, 2);
% Set any zeros to one before dividing
s = s + (s==0);
post = post./(s*ones(1, ncentres));

p1_x = reshape(post(:, 1), size(X));
p2_x = reshape(post(:, 2), size(X));

plot(samples((label==1),1),samples(label==1,2),'r^');
plot(samples((label==2),1),samples(label==2,2),'ko');
contour(xrange,yrange, p1_x,[0.5 0.5],'k-.');


axis([x0 x1 y0 y1])