image thumbnail
from RJ-MCMC algorithm for sinusoids parameter estimation by Sebastien PARIS
RJ-MCMC algorithm for sinusoids estimation

death_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , t , hyper , k , IN , out);
function [K , tp_wk , D , DtD , M , MDt , ytPy , out ] = death_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , t , hyper , k  , IN , out);

%
% [K , tp_wk , D , DtD , M , MDt , ytPy , out ] = death_move3(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , t , hyper , k  , IN , out);
%


N                 = length(y);

k1                = k(t)  + 0;

k2                = k(t)  - 1;

inv_delta         = (1/delta(: , t));



%%---- Draw an index base to remove  ------%%

proposalpos       = ceil( k1*rand );

%----------------------------------------%%


Dproposal         = [D(: , 1 : proposalpos - 1) , D(: , proposalpos + 1 : end)];      

Dproposalt        = Dproposal';

DtDproposal       = Dproposalt*Dproposal;

Mproposal         = inv(DtDproposal + inv_delta*eye(k2) );

MproposalDt       = Mproposal*Dproposalt;

ytPproposaly      = (yt*(IN - Dproposal*MproposalDt))*y;


rdeath            = ((k1*(1 + delta(t))))*((hyper.gamma + ytPy)/(hyper.gamma + ytPproposaly))^((hyper.upsilon + N)*0.5);


acceptance        = min(1 , rdeath); 


if (rand < acceptance)
    
    tp_wk                    = [tp_wk(1 : proposalpos - 1) ; tp_wk(proposalpos + 1 : k1)];
    
    K                        = k2;
    
    D                        = Dproposal;
    
    DtD                      = DtDproposal;
    
    M                        = Mproposal;
    
    MDt                      = MproposalDt;
    
    ytPy                     = ytPproposaly;
      
    out.adeath               = out.adeath + 1;
    
else
       
    K                         = k1;
        
    out.rdeath                = out.rdeath + 1;
      
end

Contact us at files@mathworks.com