maximize_params.m

Bayesian网络工具箱.

function CPD = maximize_params(CPD)
% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
% CPD = maximize_params(CPD)


% For details, see "Fitting a Conditional Gaussian Distribution", Kevin Murphy, tech. report,
% 1998, available at www.cs.berkeley.edu/~murphyk/publ.html
% Refering to table 2, we use equations 1/2 to estimate the covariance matrix in the untied/tied case,
% and equation 9 to estimate the weight matrix and mean.
% We do not implement spherical Gaussians - the code is already pretty complicated!

if ~adjustable_CPD(CPD), return; end

%assert(approxeq(CPD.nsamples, sum(CPD.Wsum)));

[self_size cpsize dpsize] = size(CPD.weights);

check_real = 1;
if check_real
  % Sometimes, with limited data, we can get imaginary components!
  % These are basically harmless, but ugly.
  if ~isreal(CPD.WXYsum) | ~isreal(CPD.WXXsum) | ~isreal(CPD.WYYsum)
    disp(['removing imaginary components for CPD ' num2str(CPD.self)]);
    CPD.WXYsum = real(CPD.WXYsum);
    CPD.WXXsum = real(CPD.WXXsum);
    CPD.WYYsum = real(CPD.WYYsum);
  end
end

% Estimate mean and weights simultaneously
% Let x2 = [x 1]' and B2 = [B mu]
XY = zeros(cpsize+1, self_size, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)' 
XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)' 
YY = zeros(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' 
for i=1:dpsize
  XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y
	       CPD.WYsum(:,i)']; % 1*Y
  % [x  * [x' 1]  = [xx' x
  %  1]              x'  1]
  XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i);
	       CPD.WXsum(:,i)'   CPD.Wsum(i)];
  YY(:,:,i) = CPD.WYYsum(:,:,i);
end

clamped_zero_mean =  CPD.clamped_mean & isequal(CPD.mean, zeros(self_size, 1, dpsize));
if ~clamped_zero_mean
  B2 = zeros(self_size, cpsize+1, dpsize);
  for i=1:dpsize
    if det(XX(:,:,i))==0  % fix by U. Sondhauss 6/27/99
      B2(:,:,i)=0;          
    else                    
      % Eqn 9 in table 2 of TR
      %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
      B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))';
    end                   
    if ~CPD.clamped_mean
      CPD.mean(:,i) = B2(:,cpsize+1,i);
    end
    if ~CPD.clamped_weights
      CPD.weights(:,:,i) = B2(:,1:cpsize,i);
    end
  end
else
  % Estimating B2 and then setting the last column (the mean) to 0s is *not* equivalent
  % to estimating B and then adding 0s to the last column.
  if ~CPD.clamped_weights
    B = zeros(self_size, cpsize, dpsize);
    for i=1:dpsize
      if det(CPD.WXXsum(:,:,i))==0
	B(:,:,i) = 0;
      else
	% Eqn 9 in table 2 of TR
	%B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i));
	B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))';
      end
    end
    CPD.weights = reshape(B, [self_size cpsize dpsize]);
  end
end

% Let B2 = [W mu]
if dpsize == 1  % bug fix due to Rainer Deventer
  if cpsize>0 
    B2(:,1:cpsize) = reshape(CPD.weights, [self_size cpsize]);
  end
  B2(:,cpsize+1) = CPD.mean;
else
  if cpsize>0
    B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]);
  end
  B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]);
end

if ~CPD.clamped_cov
  if CPD.tied_cov
    S = zeros(self_size, self_size);
    % Eqn 2 from table 2 in TR
    for i=1:dpsize
      S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
    end
    S = S / CPD.nsamples;
    if strcmp(CPD.cov_type, 'diag')
      S = diag(diag(S));
    end
    CPD.cov = repmat(S, [1 1 dpsize]);
  else 
    w = CPD.Wsum(:);
    % Set any zeros to one before dividing
    % This is valid because w(i)=0 iff WYsum(:,i)=0, etc.
    w = w + (w==0);
    for i=1:dpsize      
      % Eqn 1 from table 2 in TR
      CPD.cov(:,:,i) = (YY(:,:,i) - B2(:,:,i)*XY(:,:,i))/w(i);
    end
    if strcmp(CPD.cov_type, 'diag')
      for i=1:dpsize      
	CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i)));
      end
    end
  end
end

check_covars = 1;
min_covar = 1e-5;
if check_covars % prevent collapsing to a point
  for i=1:dpsize
    if min(svd(CPD.cov(:,:,i))) < min_covar
      disp(['resetting singular covariance for node ' num2str(CPD.self)]);
      CPD.cov(:,:,i) = CPD.init_cov(:,:,i);
    end
  end
end