function mcm1test % Figs. 9.8a, 9.8b mcm1test: Monte Carlo Method (3/3007), % nx = 1 dim, uniform dist, % I = int(F(x),x=a..b), F(x) = sqrt(1-x^2), -1 <= a < b <= +1; % technically, f(x) = (b-a)F(x) = (b-a)*sqrt(1-x^2) to account for % uniform density phi(x) = 1/(b-a) on [a,b], so I = meanf. % clc; clear % fprintf('Monte Carlo Test of 1-dim Uniform Dist. on (a,b)'); fprintf('\n with F(x)=sqrt(1-x^2) and f(x) = (b-a)F(x):\n'); a = 0; b = +1; % -1 <= a < b <= +1; % integral of f(x) = sqrt(1-x^2); on [a,b]: IntExact = 0.5*(asin(b)-asin(a))+0.5*(b*sqrt(1-b^2)-a*sqrt(1-a^2)); MufExact = IntExact; Sigf = sqrt((b-a)^2*(1-(b^2+a*b+a^2)/3)-MufExact^2); fprintf('\nk n muhatn mufExact sighatn Sigf stderrn AbsErrorf\n'); kmax = 7; n = zeros(1,kmax); meanf = zeros(1,kmax); sigf = zeros(1,kmax); sigdrn = zeros(1,kmax); error = zeros(1,kmax); for k = 1:kmax rand('state',0); % set state or seed n(k) = 10^k; % sample size, k = log10(n(k)) ; x = a+(b-a)*rand(n(k),1); % get n(k) X 1 random sample on (a,b); f = (b-a)*sqrt(1-x.^2); % vectorized f; meanf(k) = mean(f); % E[f(X)]; sigf(k) = std(f); % sqrt(sigmaf^2), sigmaf^2 = unbiased variance of f; sigdrn(k) = sigf(k)/sqrt(n(k)); error(k) = abs(meanf(k)-MufExact); fprintf('%1i %8i %6.4f %6.4f %9.3e %9.3e %9.3e %9.3e\n'... ,k,n(k),meanf(k),MufExact,sigf(k),Sigf,sigdrn(k),error(k)) 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,meanf,'k-o','linewidth',3,'MarkerSize',12); hold on plot(kv,sigf,'k-x','linewidth',3,'MarkerSize',12); axis([min(kv) max(kv) 0 1]); title('Monte Carlo Results, Uniform Dist., F(x) = sqrt(1-x^2)'... ,'Fontsize',36,'FontWeight','Bold'); xlabel('log(n), Log_{10} Sample Size'... ,'Fontsize',36,'FontWeight','Bold'); ylabel('f-Moments \mu_n, \sigma_n'... ,'Fontsize',36,'FontWeight','Bold'); legend('\mu_n, Mean-est.','\sigma_n, StdDev-est.','Location','Best'); set(gca,'Fontsize',32,'FontWeight','Bold','linewidth',3); set(gcf,'Color','White','Position' ... ,[scrsize(3)/ss(nfig) 70 scrsize(3)*0.60 scrsize(4)*0.80]); hold off % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % nfig = nfig+1; figure(nfig); kv = 1:kmax; plot(kv,log10(sigdrn),'k-o','linewidth',3,'MarkerSize',12); hold on plot(kv,log10(error),'k-x','linewidth',3,'MarkerSize',12); % ymin = min(min(log10(sigdrn)),min(log10(error))); % ymax = max(max(log10(sigdrn)),max(log10(error))); axis tight; %axis([min(kv) max(kv) ymin ymax]); title('Monte Carlo Errors, Uniform Dist., F(x) = sqrt(1-x^2)'... ,'Fontsize',36,'FontWeight','Bold'); xlabel('log(n), Log_{10} Sample Size'... ,'Fontsize',36,'FontWeight','Bold'); ylabel('f-Errors log(StdError_n), log(AbsError_n)'... ,'Fontsize',36,'FontWeight','Bold'); legend('log_{10}(StdError_n)','log_{10}(AbsError_n)','Location','Best'); set(gca,'Fontsize',32,'FontWeight','Bold','linewidth',3); set(gcf,'Color','White','Position' ... ,[scrsize(3)/ss(nfig) 70 scrsize(3)*0.60 scrsize(4)*0.80]); hold off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % end mcm1test.m %