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

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

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


N                 = length(y);

k1                = k(t)  + 0;

k2                = k(t)  + 0;


inv_delta         = 1/delta(t);

std_sRW           = sqrt(hyper.sRW);


cte1              = pi  + hyper.walk;

cte2              = 0   - hyper.walk;

cte_hyper         = (hyper.upsilon + N)*0.5;

Ik                = inv_delta*eye(k2);


u                 = rand(k1 , 1);

v                 = rand(k1 , 1);

%delta_w           = w(2) - w(1);


for j = 1 : k1
    
    Dt            = D';
    
    DtD           = Dt*D;
    
    invM          = DtD + Ik;
    
    detinvM       = det(invM);
    
    M             = inv(invM);
    
    MDt           = M*Dt;
    
    ytPy          = (yt*(IN - D*MDt))*y;
    
    if (u(j) < hyper.lambda) 
        
        proposal      = w(sum(rand > cdfperio_y)); % + delta_w*rand ;      %M-H step 1 : q1(w'_{j,k}|w_{j,k})
        
        %<-- bin index -------->%
        %wk' drawn uniformely inside the frequency
        %frontiers associated with the bin%
        
        %         proposal      = min([proposal ; cte1]);
        %         
        %         proposal      = max([proposal ; cte2]);
        
    else
        
        proposal      = tp_wk(j) + std_sRW*randn;                           %M-H step 2 : q2(w'_{j,k}|w_{j,k})
        
        proposal      = min([proposal ; cte1]);
        
        proposal      = max([proposal ; cte2]);
        
    end
    
%     if( ~all(tp_wk - proposal) )
%         
%         tp_wk(j) = tp_wk(j);
%         
%     else
        
        Dproposal                 = D;
        
        Dproposal(: , j)          = bsinu(proposal , x);
        
        Dproposalt                = Dproposal';
        
        DtDproposal               = Dproposalt*Dproposal;
        
        invMproposal              = DtDproposal + Ik;
        
        detinvMproposal           = det(invMproposal);
        
        Mproposal                 = inv(invMproposal);
        
        MproposalDt               = Mproposal*Dproposalt;
        
        ytPproposaly              = (yt*(IN - Dproposal*MproposalDt))*y;
        
        
        
        ratio                     = sqrt(abs(detinvM/detinvMproposal))*((hyper.gamma + ytPy)/(hyper.gamma + ytPproposaly))^(cte_hyper);
        
        
        acceptance                = min(1 , ratio); 
        
        if (v(j) < acceptance)
            
            tp_wk(j)              = proposal;
            
            D                     = Dproposal;
            
            DtD                   = DtDproposal;
            
            M                     = Mproposal;
            
            MDt                   = MproposalDt;
            
            ytPy                  = ytPproposaly;
            
            out.arwupdate         = out.arwupdate + 1;
            
        else
                        
            out.rrwupdate         = out.rrwupdate + 1;
            
        end
        
        %tp_wk = sort(tp_wk);
        
        %  end
end

K                            = k1; 


Contact us at files@mathworks.com