image thumbnail

Matlab Monte Carlo Demo on Advanced Statistical Estimation of Variance and Standard Deviation



Matlab Monte Carlo Demo on Advanced Statistical Estimation of Variance and Standard Deviation

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 ' 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')
    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) = ' ...
toc % how fast is your computer?
disp('pretty simple and fun, no?')

Contact us