% 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