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

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