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

birth_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] = birth_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] = birth_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 a proposal on [0 pi] ------%%

proposal          = pi*rand;

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

Dproposal         = [D , bsinu(proposal , x)];

Dproposalt        = Dproposal';

DtDproposal       = Dproposalt*Dproposal;

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

MproposalDt       = Mproposal*Dproposalt;

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


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


acceptance        = min(1 , rbirth); 

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


Contact us at files@mathworks.com