function mcm3multidim % Fig. 9.11 mcm3multidim: Monte Carlo Multidimensional Integration (3/2007), % nx = 2:5 dims, normal distribution, % I = int(F(x),x=a..b), F(x) = exp(-sum(x.^2)/2)/sqrt(2*pi)^nx, % so f(x) = I_{a<= x <= b}, an indicator function using vector inequalities. % clc; clear % fprintf('Monte Carlo Multidimenstional Integration:'); fprintf('\n including Acceptance-Rejection Technique Application,'); fprintf('\n with Normal Dist. on (a,b)'); fprintf('\n and with int(F(x),x=a..b), F(x) = exp(-x.^2/2)/sqrt(2pi)^nx;\n'); % nxmax = 5; % dimension kmax = 6; % power of 10 % f = inline ('exp(-sum(x.*x)/2) / sqrt(2*pi)^length(x)','x'); muhatn = zeros(kmax,nxmax); stderrn = zeros(kmax,nxmax); for nx = 2:nxmax a = -2*ones(1,nx); % lower vector limit b = 2*ones(1,nx); % upper vector limit nmc = zeros(1,kmax); for k = 1:kmax nmc(k) = 10^k; % sample size nac = 0; for n = 1:nmc(k) x = randn(1,nx); % MATLAB 1Xnmc vector normal distribution; if (x >= a) & (x <= b) % von Neumann accept-reject technique nac = nac + 1; % counts accepted points; end end muhatn(k,nx) = nac/nmc(k); stderrn(k,nx) = sqrt(muhatn(k,nx)*(1-muhatn(k,nx))/(nmc(k)-1)); end end % fprintf('\nMonte Carlo results in mutlidimension,'); fprintf('\n by von Neumann''s Acceptance-Rejection technique:\n'); fprintf('\nMonte Carlo Mean Estmate, muhatn:'); fprintf('\n k n nx=2 nx=3 nx=4 nx=5\n'); for k = 1:kmax fprintf ('%2i %8i %8.6f %8.6f %8.6f %8.6f\n'... ,k,nmc(k),muhatn(k,2:nxmax)); end fprintf('\nMonte Carlo Std Error Estmate, stderrn:'); fprintf('\n k n nx=2 nx=3 nx=4 nx=5\n'); for k = 1:kmax fprintf ('%2i %8i %9.3e %9.3e %9.3e %9.3e\n'... ,k,nmc(k),stderrn(k,2:nxmax)); end % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% scrsize = get(0,'ScreenSize'); ss = [3.0,2.8,2.6,2.4,2.2,2.0]; % nfig = 1; figure(nfig); kv = 1:kmax; plot(kv,muhatn(:,2),'k-o'... ,kv,muhatn(:,3),'k-x'... ,kv,muhatn(:,4),'k-+'... ,kv,muhatn(:,5),'k-*'... ,'linewidth',3,'MarkerSize',14); axis([min(kv) max(kv) 0.5 1]); title('Monte Carlo Means, Normal Dist., F(x) = \phi_n(x) on [a,b]'... ,'Fontsize',44,'FontWeight','Bold'); xlabel('log(n), Log_{10} Sample Size'... ,'Fontsize',44,'FontWeight','Bold'); ylabel('Mean Estimates, \mu_n'... ,'Fontsize',44,'FontWeight','Bold'); hlegend=legend('nx = 2','nx = 3','nx = 4','nx = 5','Location','Best'); set(hlegend,'Fontsize',36,'FontWeight','Bold'); set(gca,'Fontsize',36,'FontWeight','Bold','linewidth',3); set(gcf,'Color','White','Position' ... ,[scrsize(3)/ss(nfig) 70 scrsize(3)*0.60 scrsize(4)*0.80]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % nfig = nfig+1; figure(nfig); kv = 1:kmax; plot(kv,stderrn(:,2),'k-o'... ,kv,stderrn(:,3),'k-x'... ,kv,stderrn(:,4),'k-+'... ,kv,stderrn(:,5),'k-*'... ,'linewidth',3,'MarkerSize',14); axis tight; % axis([min(kv) max(kv) 0 0.2]); title('Monte Carlo Std. Errors, Normal Dist. on [a,b]'... ,'Fontsize',44,'FontWeight','Bold'); xlabel('log(n), Log_{10} Sample Size'... ,'Fontsize',44,'FontWeight','Bold'); ylabel('Std. Errors, stderr_n'... ,'Fontsize',44,'FontWeight','Bold'); hlegend=legend('nx = 2','nx = 3','nx = 4','nx = 5','Location','Best'); set(hlegend,'Fontsize',36,'FontWeight','Bold'); set(gca,'Fontsize',36,'FontWeight','Bold','linewidth',3); set(gcf,'Color','White','Position' ... ,[scrsize(3)/ss(nfig) 70 scrsize(3)*0.60 scrsize(4)*0.80]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % end mcm3multidim.m %