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
%