第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
蒙特卡洛法matlab实现.docx
运行顺序
/gf_mmc_driver.m
/mcmc1.m
/mcmc2.m
gf_mmc_driver.m
s
% "This is the driver script to be run from the command line. Also
% requires functions mcmc1 and mcmc2. This code is referred to by an
% upcoming article by T. A. Driscoll and K. Maki in the Education
% section of SIAM Review, to appear 2006-2007."
%
% MMC iteration applied to the growth factors of random matrices.
% Runs 'numberchain' instances of MMC simultaneously using the Distributed
% Computing Toolbox.
%
N = 8;
% size of nxn matrix
binedge = -.5:1:(2^(N-1)+.5);
% edges of target bins
Nbin = length(binedge) - 1;
% number of bins
M = 5e5;
% number of realizations per iteration
Niter = 70;
% number of iterations
burnin = 5e5;
% burn-in period for MCMC
numberchain = 12;
% number of independent MMC runs
% ***** Initialization *****
randn('state',sum(100*clock)), rand('state', sum(100*clock));
Accept = zeros(Niter,numberchain);
% acceptance rate per iteration per run
sigma = (1/N)*ones(numberchain,1); % sigma parameter per run
H = zeros(Nbin,Niter,numberchain); % hits/bin per iteration per run
P = ones(Nbin,numberchain)/Nbin;
% computed histogram per run
g = zeros(Nbin-1,Niter,numberchain);% parameter needed for Berg's iteration
Pnew = zeros(Nbin,1);
% updated inverse weights
GFinitial = zeros(numberchain,1);
% computed initial growth fac for MC
X = zeros(N,N,numberchain);
% initial matrices for MC runs
chaingrowth = 10*(1:numberchain-2); % desired initial growth fac for MC
% Establishing communication for distributed computing
jm = findResource('jobmanager','Name','woprdce');
% Calculating one initial matrix for Markov Chain
X(:,:,1) = randn(N,N);
% ***** Start MMC iteration *****
for i=1:Niter
% Calculating remaining initial matrices for Markov Chain
X(:,:,2)=randn(N);
% Implementing Higham and Higham's Theorem
for z=3:numberchain
while( GFiter(z)>chaingrowth(z-2)+.5 || GFiter(z) < chaingrowth(z-2)-.5 )
Mstart = toeplitz([1 -.99999*(-1+chaingrowth(z-2)^(1/(N-1)))*ones(1,N-1)], [1
zeros(1,N-1)]);
T = triu( (2*rand(N-1)-1)/N );
T = [T;zeros(1,N-1)];
Tlc = (chaingrowth(z-2).^((0:N-1)/(N-1)));
D = diag( sign(randn(N,1)) );
X(:,:,z) = D*Mstart*[T Tlc'];
[L,U]=lu(X(:,:,z));
GFiter(z) = (max(abs(U(:))))/(max(max(abs(X(:,:,z)))));
end
end
% ***** Running MCMC for burn-in *****
% Telling the Distributed Computing Toolbox to complete one job with
% 'numberchain' tasks. Each task is comprised of running a MCMC
% for the burnin time (mcmc1.m) with a different initial matrix.
j = createJob(jm);
set(j,'Timeout',1200);
set(j,'FileDependencies',{ 'mcmc1.m'});
for z = 1:numberchain
createTask( j, @mcmc1, 2, {X(:,:,z), burnin, sigma(z), P(:,z), binedge});
end
submit(j);
waitForState(j,'finished');
data = getAllOutputArguments(j);
% Received computed data and recording results
for z = 1:numberchain
Accept(i,z) = data{z,1};
X(:,:,z) = data{z,2};
end
clear j
% Adjusting the parameter sigma based on burnin results
for z = 1:numberchain
if( Accept(i,z)/burnin < .1 )
sigma(z) = max( sigma(z)/2, eps);
elseif( Accept(i,z)/burnin > .4 )
sigma(z) = 2*sigma(z);
else
sigma(z) = sigma(z);
end
end
% ***** Running MCMC for M *****
% Again telling the Distributed Computing Toolbox to complete one job
% composed of 'numberchain' tasks. Each task continues running the
% independent MCMC iterations (script mcmc2.m) from above. Now
% statistics are calculated.
j = createJob(jm);
set(j,'Timeout',1200);
set(j,'FileDependencies',{'mcmc2.m'});
for z=1:numberchain
createTask(j,@mcmc2, 3, {X(:,:,z), M, sigma(z), P(:,z), binedge});
end
submit(j);
waitForState(j,'finished');
data = getAllOutputArguments(j);
% Received computed data and recording results
for z = 1:numberchain
Accept(i,z) = data{z,1};
H(:,i,z) = data{z,2};
X(:,:,z) = data{z,3};
end
clear j
% Tuning the parameter sigma again
for z = 1:numberchain
if( Accept(i,z)/(M) < .1 )
sigma(z) = max( sigma(z)/2, eps);
elseif( Accept(i,z)/(M) > .4 )
sigma(z) = 2*sigma(z);
else
sigma(z) = sigma(z);
end
end
% *****Berg update*****
for z =1:numberchain
Pnew = P(:,z);
for k = 1:(Nbin-1)
if (H(k+1,i,z)*H(k,i,z) > 0)
g(k,i,z) = H(k+1,i,z)*H(k,i,z)/(H(k+1,i,z)+H(k,i,z));
ghat(k) = g(k,i,z)/(sum(g(k,1:i,z)));
Pnew(k+1) = Pnew(k)*(P(k+1,z)/P(k,z))*((H(k+1,i,z)/H(k,i,z))^ghat(k));
else
Pnew(k+1) = Pnew(k)*(P(k+1,z)/P(k,z));
end
end
P(:,z) = Pnew/sum(Pnew);
end
end
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
mcmc2( X, M, sigma, P, binedge )
% Run of MCMC --- statistics recorded
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
