disp 'Monte Carlo Demo - very easy to run - just type the file name'
disp 'Uploaded January 1, 2012'
disp 'Author: Amitava Biswas, PhD, FIEI, FIETE'
disp 'The University of Texas at El Paso'
disp 'Email:abiswas@utep.edu Phone:(915)747-8307'
disp 'this demo shows an interesting fact from advanced statistics'
disp 'sample variance is expected to be population variance'
disp 'but its square root is expected below standard deviation'
disp 'to get unbiased estimate of standard deviation mutiply by'
disp 'gamma(n/2-1/2)*sqrt(n/2-1/2)/gamma(n/2)'
close all; clear all; tic % how fast is your computer?
n=4; % sample size? change if necessary
scrsz = get(0,'ScreenSize'); % full screen looks better
figure('Position', scrsz, 'Units', 'normalized');
sumvar=0; % initialize sum of variances observed
sumstd=0; % initialize sum of standard deviations observed
k=10000; % how many times to try your luck?
for m=1:k
clf % erase all previous cases
axes('Position',[0 0 1 1], 'Units','normalized');
axis([-2 2 -2 2]); % broad enough, change if necessary
axis on; grid on; % looks good
sample=zeros(1,n); % reset
for p=1:n % draw n samples Monte Carlo style
sample(p)=randn; % just a sample from standard normal distribution
switch fix(4*rand) % select a nice color
case 0; text(sample(p), randn,'\clubsuit','fontsize',200,'color','green')
case 1; text(sample(p), randn,'\diamondsuit','fontsize',200,'color','magenta')
case 2; text(sample(p), randn,'\heartsuit','fontsize',200,'color','red')
otherwise; text(sample(p), randn,'\spadesuit','fontsize',200,'color','black')
end
end
sumvar=sumvar+var(sample); % accumulate sum to check for mean bias
sumstd=sumstd+sqrt(var(sample)); % accumulate sum to check for mean bias
text(-1.8,1.5,['Monte Carlo Demo by Amitava Biswas'], 'fontsize',24);
text(-1.8,0.5,['sample size = n = ' num2str(n)], 'fontsize',24);
text(-1.8,0,['cycle counter ' num2str(k-m) ' + ' num2str(m)], 'fontsize',24);
text(-1.8,-0.5,['average VARIANCE = ' num2str(sumvar/m,4) ...
', and bias = ' num2str(100*sumvar/m-100,2) '%'],'fontsize',24)
text(-1.8,-1,['average STANDARD DEVIATION = ' num2str(sumstd/m,4) ...
', and bias = ' num2str(100*sumstd/m-100,2) '%'],'fontsize',24)
text(-1.8,-1.5, ['better multiply by gamma(n/2-1/2)*sqrt(n/2-1/2)/gamma(n/2) = ' ...
num2str(gamma(n/2-1/2)*sqrt(n/2-1/2)/gamma(n/2),4)],'fontsize',24)
drawnow
end
toc % how fast is your computer?
disp('pretty simple and fun, no?')
