propose.m.txt

这是马尔可夫-蒙特卡罗算法的MATLAB源程序.

function q = Propose(p,m)
% q = Propose(p,m)
% propose q given state p and model m

% Colin Fox, 19 Jan 2003

ratio = [1 4 1 4]; % ratio of move probabilities
%moveprob = ratio/sum(ratio);
%movenum = disample(moveprob);
movenum = min(find(rand<cumsum(ratio/sum(ratio))));

switch movenum
case 1
    % Move # 1 is pick a resistor and random and set it to a possible value at random
    [mr,nr] = find(p.DR);               % list of resistor indices
    ind = PickOneFromM(length(mr));     % pick an index at random
    NewR = PickOne(m.Resistances);      % pick an allowable resistance value at random
    q = struct('mr',mr(ind),'nr',nr(ind),'NewR',NewR,'OldR',p.DR(mr(ind),nr(ind)),'MoveType',1,'lg',NaN);
case 2
    % Move # 2 is pick a node (**with different-valued resistors** not implemented) at random, pick a pair
    % of resistors on that node, at random, and swap resistor values
    [nn,mm] = size(p.DR);               % nn (or mm) is the number of nodes
    nind = PickOneFromM(nn);            % pick a node at random
    [mr,nr] = find(p.DR);               % list of resistor indices
    ind1 = find(mr == nind);            % index of resistors "from" node
    ind2 = find(nr == nind);            % index of resistors "to" node
    inds = ChooseN([ind1;ind2],2);      % pick a pair of resistors
    OldR = [p.DR(mr(inds(1)),nr(inds(1))) ; p.DR(mr(inds(2)),nr(inds(2)))]; % existing resistances
    NewR = [p.DR(mr(inds(2)),nr(inds(2))) ; p.DR(mr(inds(1)),nr(inds(1)))]; % swapped over
    q = struct('mr',mr(inds),'nr',nr(inds),'NewR',NewR,'OldR',OldR,'MoveType',2,'lg',NaN);
case 3
    % Pick two resistors at random and swap them
    [mr,nr] = find(p.DR);               % list of resistor indices
    inds = ChooseNFromM(length(mr),2);  % pick a pair of indices at random
    OldR = [p.DR(mr(inds(1)),nr(inds(1))) ; p.DR(mr(inds(2)),nr(inds(2)))]; % existing resistances
    NewR = [p.DR(mr(inds(2)),nr(inds(2))) ; p.DR(mr(inds(1)),nr(inds(1)))]; % swapped over
    q = struct('mr',mr(inds),'nr',nr(inds),'NewR',NewR,'OldR',OldR,'MoveType',3,'lg',NaN);
case 4
    % Pick a resistor at random, then a resistor at each end, and swap the end resistors
    [mr,nr] = find(p.DR);               % list of resistor indices
    ind = PickOneFromM(length(mr));     % pick an index at random
    mind = mr(ind); % one end
    nind = nr(ind); % the other end
    
    mres = PickOne(find((mr == mind & nr ~= nind) | (nr == mind & mr ~= nind)));
    nres = PickOne(find((mr == nind & nr ~= mind) | (nr == nind & mr ~= mind)));
    OldR = [p.DR(mr(mres),nr(mres)) ; p.DR(mr(nres),nr(nres))]; % existing resistances
    NewR = flipud(OldR); % swapped over
    q = struct('mr',mr([mres nres]),'nr',nr([mres nres]),'NewR',NewR,'OldR',OldR,'MoveType',4,'lg',NaN);
end