No BSD License  

Highlights from
Multicanonical Monte Carlo scheme for finding rare growth factors

image thumbnail

Multicanonical Monte Carlo scheme for finding rare growth factors

by

 

05 Oct 2006 (Updated )

A multicanonical Monte Carlo (MMC) scheme developed to approximate the tails of growth factor probab

mcmc2( X, M, sigma, P, binedge )
% Run of MCMC --- statistics recorded

function [ Accept, Hits, X ] = mcmc2( X, M, sigma, P, binedge )

    acceptval = rand(M,1);                  % all acceptance test
    N = length(X);                          % size of matrix
    Accept = 0;                             % number of acceptances
    Hits = zeros( length(binedge)-1,1 );    % hits/bin for current iteration

    % Calculating initial growth factor
    [L,U] = lu(X);
    growth = max( abs(U(:)) )/max( abs(X(:)) );
    [junk, bin] = histc(growth,binedge);

    % *****MCMC Starts*****
    for m=1:M        
        
        Xpro = X + sigma*trnd(8,N,N);       % proposal matrix
        
        % Calculating growth factor of proposal
        [L,Upro] = lu(Xpro);
        growthXpro = max( abs(Upro(:)))/max( abs(Xpro(:)) );
        [junk, newbin] = histc(growthXpro,binedge);
        
        % Calculating pi(X) and pi(X_pro)
        pi_x = exp(-.5*sum(sum(reshape( X.^2, [N^2 1]))));
        pi_y = exp(-.5*sum(sum(reshape( Xpro.^2, [N^2 1]))));
        
        % Accept proposal?
        alpha = min(1, ((pi_y*P(bin))/(pi_x*P(newbin))));
        if acceptval(m) < alpha
            X = Xpro;
            bin = newbin;
            growth = growthXpro;
            Accept = Accept + 1;
        end
        
        % Record statistics
        Hits(bin) = Hits(bin)+1;
    end
end

Contact us