radialDeath.m
上传用户:fsbooksir
上传日期:2013-10-19
资源大小:14k
文件大小:2k
源码类别:

matlab例程

开发平台:

Matlab

  1. function [k,mu,M,aDeath,rDeath] = radialDeath(aDeath,rDeath,k,mu,M,delta,x,y,hyper,t,nabla);
  2. %
  3. if nargin < 11, error('Not enough input arguments.'); end
  4. [N,d] = size(x);      % N = number of data, d = dimension of x.
  5. [N,c] = size(y);      % c = dimension of y, i.e. number of outputs.
  6. insideDeath=1;
  7. uD=rand(1);
  8. % INITIALISE H AND P MATRICES:
  9. % ===========================
  10. invH=zeros(k(t)+1+d,k(t)+1+d,c);
  11. P=zeros(N,N,c);
  12. invHproposal=zeros(k(t)+d,k(t)+d,c);
  13. Pproposal=zeros(N,N,c);
  14. for i=1:c,
  15.   invH(:,:,i) = (M'*M + (1/delta(t,i))*eye(k(t)+1+d));
  16.   P(:,:,i) = eye(N) - M*inv(invH(:,:,i))*M';
  17. end;
  18. % CHOOSE UNIFORMLY A BASIS FUNCTION TO BE DELETED:
  19. % ===============================================
  20. proposalPos= d+1+unidrnd(length(mu{t}(:,1)),1,1);
  21. if (proposalPos==d+1+k(t)),
  22.   Mproposal = [M(:,1:proposalPos-1)];     
  23. else
  24.   Mproposal = [M(:,1:proposalPos-1) M(:,proposalPos+1:k(t)+d+1)];      
  25. end;
  26. for i=1:c,
  27.   invHproposal(:,:,i) = (Mproposal'*Mproposal + inv(delta(t,i))*eye(k(t)+d)); 
  28.   Pproposal(:,:,i) = eye(N) - Mproposal*inv(invHproposal(:,:,i))*Mproposal'; 
  29. end;
  30. % PERFORM A DEATH MOVE:
  31. % ====================
  32. ratio= k(t) * sqrt(delta(t,1)) * sqrt((det(invH(:,:,1)))/(det(invHproposal(:,:,1)))) * ((hyper.gamma+y(:,1)'*P(:,:,1)*y(:,1))/(hyper.gamma + y(:,1)'*Pproposal(:,:,1)*y(:,1)))^((hyper.v+N)/2);  
  33. for i=2:c,
  34.   ratio= ratio * sqrt(delta(t,i)) * sqrt((det(invH(:,:,i)))/(det(invHproposal(:,:,i)))) * ((hyper.gamma+y(:,i)'*P(:,:,i)*y(:,i))/(hyper.gamma + y(:,i)'*Pproposal(:,:,i)*y(:,i)))^((hyper.v+N)/2);  
  35. end;
  36. acceptance = min(1,ratio);  
  37. if (uD<acceptance),
  38.   previousMu = mu{t};
  39.   if (proposalPos==(1+d+1)),
  40.     mu{t+1} = [previousMu(2:k(t),:)]; 
  41.   elseif (proposalPos==(1+d+k(t))),
  42.     mu{t+1} = [previousMu(1:k(t)-1,:)];
  43.   else
  44.     mu{t+1} = [previousMu(1:proposalPos-1-d-1,:); previousMu(proposalPos-d-1+1:k(t),:)];
  45.   end;
  46.   k(t+1) = k(t)-1;
  47.   M=Mproposal;
  48.   aDeath=aDeath+1;
  49. else
  50.   mu{t+1} = mu{t};
  51.   k(t+1) = k(t);
  52.   rDeath=rDeath+1;
  53.   M=M;
  54. end;