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

mcmc1( X, burnin, sigma, P, binedge )
% Run MCMC for the burnin period.


function [ Accept, X ] = mcmc1( X, burnin, sigma, P, binedge )

    acceptval = rand(burnin,1);     % all acceptance tests
    N = length(X);                  % size of matrices
    Accept = 0;                     % number of matrices accept during run
   
    % Calculating growth factor of initial matrix
    [L,U] = lu(X);                    
    growth = max( abs(U(:)) )/max( abs(X(:)) );
    [junk, bin] = histc(growth,binedge);
    
    % *****MCMC starts*****
    for m = 1:burnin         
        Xpro = X + sigma*trnd(8,N,N);       % proposal matrix
        
        % calculating growth factor of proposal matrix
        [L,Upro] = lu(Xpro);      
        growthXpro = max( abs(Upro(:)))/max( abs(Xpro(:)) );
        [junk, newbin] = histc(growthXpro,binedge);
        
        % calculation 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
    end
end

Contact us