% 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