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

merge_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] = merge_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] = merge_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] ------%%

position        = ceil(k1*rand);

wk1             = tp_wk(position);

temp_wk         = wk1 - tp_wk;

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

dist(position)  = Inf;


position2       = find(dist == min(dist));

position2       = position2(1); % case where there are several solutions for a second basis

wk2             = tp_wk(position2);

wkg             = 0.5*(wk1 + wk2);



if (position2 > position)
    
    
    Dproposal = [D(: , 1:position - 1) D(: , position + 1: position2 - 1) D(: , position2 + 1 : k1)];      
    
else
    
    Dproposal = [D(: , 1:position2 - 1) D(: , position2 + 1: position - 1) D(: , position + 1 : k1) ];      
    
end

Dproposal         = [Dproposal,  bsinu(wkg , x)];

Dproposalt        = Dproposal';

DtDproposal       = Dproposalt*Dproposal;

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

MproposalDt       = Mproposal*Dproposalt;

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


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


acceptance        = min(1 , ratio);  

if (min(dist) < 2*hyper.merge)
    
    acceptance = 0;   
    
end


if (rand < acceptance)
    
    
    if (position2 > position)
             
        wk_temp = [tp_wk(1:position - 1) ; tp_wk(position + 1: position2 - 1) ; tp_wk(position2 + 1 : k1)];      
        
    else
        
        wk_temp = [tp_wk(1:position2 - 1) ; tp_wk(position2 + 1: position - 1) ; tp_wk(position + 1 : k1) ];      
        
    end
        
    tp_wk                      = [wk_temp ; wkg];
    
    K                          = k2;
    
    D                          = Dproposal;
    
    DtD                        = DtDproposal;
    
    M                          = Mproposal;
    
    MDt                        = MproposalDt;
    
    ytPy                       = ytPproposaly;
    
    out.amerge                 = out.amerge + 1;
    
else
        
    K                          = k1;
        
    out.rmerge                 = out.rmerge + 1;
    
end

Contact us at files@mathworks.com