amde.m

Adaptive Mixtures Density Estimation的matlab实现

% Adaptive Mixtures Density Estimation
%
%   [WTS,MUS,VARS,NTERMS] = AMDE(DATA,MT,TC)
%
%   This implements the multivariate adaptive mixtures density
%   estimation procedure. It returns the model based on the
%   DATA (nxd) with MT maximum number of terms. TC gives the
%   create threshhold based on the Mahalnaobis distance of the point
%   with all of the current terms.

function [pies,mus,vars,nterms] = amde(data,maxterms,tc)

[n,d] = size(data);

if nargin ==2
    switch d
    case 2
        %tc = 2.28;     % for 1-sigma
        tc = 2;
    case 3
        tc = 3;
        %tc = 3.54;     % for 1-sigma
    case 4
        tc = 4;
        %tc = 5.4;      % for 1-sigma
    otherwise
        tc = 5;
        %tc = 6;
    end
end

% get constants, set up vectors
mus = zeros(d,maxterms);	%each col is a term
vars = zeros(d,d,maxterms);	%each dxd page (3rd dimensions) is term
pies = zeros(1,maxterms);
posterior = zeros(1,maxterms);
minpie = .00001;		% lower bound on new pies
sievebd = 1000;		% bounding the parameter space

% initialize density to first data point
  nterms = 1;
  mus(:,1) = data(1,:)';
  vars(:,:,1) = cov(data)/2.5;	% this is 'rule of thumb'
%  vars(:,:,1) = eye(d,d);
  pies(1) = 1; 
  
  for i = 2:n
      posterior(1:nterms) = rpostupm(data(i,:)',pies,mus,vars,nterms);
      
      %if any(posterior(1:nterms) > tc)    % update terms, since at least one term 'fits' the point.
          
          
    if ~createm(data(i),pies,mus,vars,tc,nterms,maxterms)	% update terms
               
         vars = rvarupm(data(i,:),pies,mus,vars,posterior,nterms,i);
          mus = rmuupm(data(i,:),pies,mus,posterior,nterms,i);
          pies(1:nterms) = rpieupm(posterior,pies,i,nterms);
      
%	vars(:,1:d*nterms)=rvarupm(data(i,:),pies,mus,vars,posterior,nterms,i);
%	mus(:,1:nterms)=rmuupm(data(i,:),pies,mus,posterior,nterms,i);

    else	% create a new term
        
        nterms = nterms+1;
        mus(:,nterms) = data(i,:)';
        pies(nterms) = max([1/(i),minpie]);	% update pies
        pies(1:nterms-1) = pies(1:nterms-1)*(1-pies(nterms));
        vars(:,:,nterms) = setvarm(mus,pies,vars,data(i,:),nterms-1);
        
    end 	% end if create statement

  % to prevent spiking of variances
  index = find(vars(1:nterms)<1/(sievebd*nterms));
  vars(index) = ones(size(index))/(sievebd*nterms);

  end  % for i loop

% clean things up.
pies(nterms+1:end) = [];
mus(:,nterms+1:end) = [];
vars(:,:,nterms+1:end) = [];

%%%%%%%%%%%%FUNCTION - UPDATE POSTERIORS %%%%%%%%%%%%%%%%%%%%%%
function posterior=rpostupm(x,pies,mus,vars,nterms)
totprob=0;
posterior=zeros(1,nterms);
for i=1:nterms	%loop to find total prob in denominator (hand, pg 37)
  posterior(i)=pies(i)*evalnorm(x',mus(:,i)',vars(:,:,i));
  totprob=totprob+posterior(i);
end
posterior=posterior/totprob;

%%%%%%%%%%%%%%%  FUNCTION EVALNORM %%%%%%%%%%%%%
function prob = evalnorm(x,mu,cov_mat);
[n,d]=size(x);
prob = zeros(n,1);
a=(2*pi)^(d/2)*sqrt(det(cov_mat));
covi = inv(cov_mat);
for i = 1:n
	xc = x(i,:)-mu;
	arg=xc*covi*xc';
	prob(i)=exp((-.5)*arg);
end
prob=prob/a;

%%%%%%%%%%%%%  FUNCTION RVARUPM %%%%%%%%%%%%%%%%%%%%%%%%%%%
function vars = rvarupm(x,pies,mus,vars,posterior,nterms,n)
[n,d] = size(x);
inertvar=10;
betan=1/(n);
sievebd=1000;		% bounding the parameter space
for i=1:nterms
  denom=(1/betan)*pies(i)+inertvar;
  vars(:,:,i)=vars(:,:,i)+posterior(i)*((x'-mus(:,i))*(x'-mus(:,i))'...
  	-vars(:,:,i))/denom;
  if det(vars(:,:,i))<1/(sievebd*nterms);
  	vars(:,:,i)=eye(d,d)*sqrt(1/(sievebd*nterms));
  end
end

%%%%%%%%%%%%%%%%%%%  FUNCTION RMUUPM  %%%%%%%%%%%%%%%%%%%%
function mus=rmuupm(x,pies,mus,posterior,nterms,n)

betan=1/(n);
for i=1:nterms
  mus(:,i)=mus(:,i)+betan*posterior(i)*(x'-mus(:,i))/pies(i);
end

%%%%%%%%%%%%%%%   FUNCTION RPIEUPM  %%%%%%%%%%%%%%%%%%
function piess=rpieupm(posterior,pies,n,nterms)
betan = 1/(n);
post = posterior(1:nterms);
piess = pies(1:nterms);
piess = piess+betan*(post-piess);

%%%%%%%%%%%%%  FUNCTION SETVARM %%%%%%%%%%%%%%%%%%%%%%

% this function will update the variances
% in the recursive amde
% call with nterms-1, since new term is based only on
% previous terms stuff
function newvar = setvarm(mus,pies,vars,x,nterms)
totprob = 0;
sievebd = 1000;		% bounding the parameter space
posterior = zeros(1,nterms);
[d,c] = size(mus);
newvar = zeros(d,d);
for i = 1:nterms	%loop to find total prob in denominator (hand, pg 37)
  posterior(i) = pies(i)*evalnorm(x,mus(:,i)',vars(:,:,i));
  totprob = totprob+posterior(i);
end
posterior = posterior/totprob;
for i = 1:nterms
  newvar = newvar+posterior(i)*vars(:,:,i);
  if det(vars(:,:,i))<1/(sievebd*nterms);
  	vars(:,:,i) = eye(d,d)*sqrt(1/(sievebd*nterms));
  end
end

%%%%%%%%%%%%%%  CREATEM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  function true = createm(datapt,pies,mus,vars,tc,nterms,maxterms)
% different create rules can be used in here as 
% desired.  just change the code appropriately
function true=createm(datapt,pies,mus,vars,tc,nterms,maxterms)

% compute functional value at data point
true=0;
% this rule says that if the point is greater than
% tc sigmas away from each term, then we create
distance=zeros(1,nterms);
for i=1:nterms
  distance(i)=(datapt'-mus(:,i))'*inv(vars(:,:,i))*(datapt'-mus(:,i));
end
if min(distance)>tc & nterms < maxterms
	true=1;
end