main.m.txt

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

% Main line for resistor network sampler 

% Colin Fox 10 Jan 2003, modified 22 June 2003
% Simplified version for conversion to python, 29 Nov 2005

% flags
plotLP = 1;            % 1=plot log Likelihood and log Prior each SweepsPerSave
ShowState = 1;         % 1=display the resistor network each SweepsPerSave

nsamp = 4000;           % generate this many samples
SweepsPerSave = 200;   % save state & statistics after this many sweeps

% initialize vectors to store results
countmax = ceil(nsamp/SweepsPerSave);
llvec = zeros(1,countmax);
lpvec = zeros(1,countmax);
tvec = zeros(1,countmax); % to store CPU time

rand('state',0); % reset random number generator

% initialize model, data, permutations for solver
m = SetModel([1 1.5],0.5);                      % set model with resistance values and lumping constant
d = LoadData('Data24By24');
%pperm=1:624;rperm=1:624; % remove permutations as test   
[pperm,rperm] = precompute_symmmd(d); % optimal permutation to reduce bandwidth

p = SetInitialState(d,m,pperm,rperm);         % initial state
tlP = logPrior(d.trueDR,m,d.lic);               % log prior of true image

% initialize statistics
lLmax = -inf; MLEState = d.trueDR; % initialise maximum likelihood state
lPmax = -inf; MAPState = d.trueDR; % initialise maximum a posteriori state
[ni,mi,ri] = find(d.trueDR); [nn,mm] = size(d.trueDR);
MEANState = sparse(ni,mi,0,nn,mm); MSState = sparse(ni,mi,0,nn,mm); meancount = 0;%mean and mean-square

t0 = clock; % measure time elapsed, not quite CPU time

disp('entering main loop')
for t = 2:nsamp                                     % main loop
    q = Propose(p,m);
    if ~all(q.NewR == q.OldR)                       % reject null proposal
        [lalpha,q] = Alpha(p,q,d,m,pperm,rperm);    % MH correction using accurate Likelihood
        if rand < exp(lalpha)                       % accept
            p = Accept(p,q);                        % update state structure
        end
    end
    disp(['completed step # ',num2str(t)])
    
    % collect statistics, and display state
    if mod(t,SweepsPerSave)==0
        count = t/SweepsPerSave;
        llvec(count) = p.logLikelihood;
        lpvec(count) = p.logPrior;
        tvec(count) = etime(clock,t0); % elapsed time
        if plotLP PlotStats(llvec,lpvec,tlP,d.truelogLikelihood,count,countmax), end %else, disp(t), end
        if ShowState, figure(1),DisplayRes(p.DR,m,p.G), title(['State ',int2str(t)]), drawnow, end
        % collect statistics
        % mean and mean square
        if count > 200
            MEANState = MEANState + p.DR;
            MSState = MSState + (p.DR).^2;
            meancount = meancount+1;
        end
        % MLE and MAP
        if p.logLikelihood > lLmax, lLmax = p.logLikelihood; MLEState = p.DR; end
        if p.logLikelihood+p.logPrior > lPmax, lPmax = p.logLikelihood+p.logPrior; MAPState = p.DR; end
    end
end