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

Kara Maki (view profile)

 

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