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

split_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , walkint , t , hyper , k , IN , out);
function [K , tp_wk , D , DtD , M , MDt , ytPy , out] = split_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , walkint , 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] ------%%

position          = ceil(k1*rand);

proposal          = tp_wk(position);

uu                = rand*hyper.merge;

wk1               = proposal - uu;

wk2               = proposal + uu;


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

Dproposal         = [D(: , 1: position - 1) D(: , position + 1 : k1)];      

Dproposal         = [Dproposal , bsinu(wk1 , x) ,  bsinu(wk2 , x)];

Dproposalt        = Dproposal';

DtDproposal       = Dproposalt*Dproposal;

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

MproposalDt       = Mproposal*Dproposalt;

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


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


acceptance        = min(1 , rsplit); 


dist2             = sqrt(sum((wk1 - wk2).^2 , 2));

temp_wk           = wk1 - tp_wk(1 : k1);

dist1             = sqrt(sum( temp_wk.*temp_wk , 2));

dist1(position)   = Inf;

if(any(dist1 < dist2))
    
    acceptance = 0;
    
end


if (rand < acceptance)
    
    previouswk                 = tp_wk;
    
    wktrunc                    = [previouswk(1 : position - 1) ; previouswk(position  + 1 : k1)];
    
    tp_wk                      = [wktrunc ; wk1 ; wk2];
    
    K                          = k2;
    
    D                          = Dproposal;
    
    DtD                        = DtDproposal;
    
    M                          = Mproposal;
    
    MDt                        = MproposalDt;
    
    ytPy                       = ytPproposaly;
    
    out.asplit                 = out.asplit + 1;
    
else
        
    K                          = k1;
        
    out.rsplit                 = out.rsplit + 1;
    
end

Contact us at files@mathworks.com