| 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
|
|