% LMEx040304_4th.m
% Example 4.3.4 Effects of rounding analyzed with normal approximation
% Effect of rounding on the expected value and variance of rounded 
% estimates. P. 304 in
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Using the central limit theorem, what is the probability that the error
% on 100 transactions, each rounded to the nearest $100, would exceed $500?
% What is the variance of a uniformly distributed random number
% over the interval -50 to 50?
% From definition of Variance:
% Var(Y)=E[(Y-mu)^2]= integral from -50 to 50 {(y-mu)^2 f_y(y)} dy.
% mu=0;
% f_y(y)=1; This is the height of the uniform pdf at each value of y.
% for max(y)-min(y)=1, the value of f_y(y) is 1.
% for uniform random numbers spread over different ranges, 
% for max(y)-min(y)*f_y(y) must equal 1.
% thus, for a uniform random number distributed on the interval -100 to 100,
%                          f_y(y)=1/200;
% Written by Eugene.Gallagher@umb.edu for EEOS601
% http://www.es.umb.edu/edgwebp.htm
% Written in 2001, revised 1/9/11
% The m.file goes through some applications of the symbolic math toolbox at
% the beginning to verify that the variance of 100 rounded transactions is
% 2500/3 
format rat
syms y mu
int(y^2/100,-50,50)%Use the symbolic math toolbox to find the 
                          % variance
                   % of a uniformly distributed variable on the interval
                   % -50 to 50; it is 2500/3 
VarY=eval(int(sym('y^2/100'),-50,50))
% Variance of a uniform randomly distributed number
% on the interval -50 to 50 is 2500/3
% Note how the variance remains 2500/3 with a shift in mean:
% Var(Y)=E(Y^2)-mu^2;
mu=50;
int(y^2/100,0,100)
EYsq=eval(int(sym('y^2/100'),0,100))
% The previous statement produces 10,000/3 which produces an identical
% variance after subtracting the square of the expected value, mu^2, as is
% appropriate from the definition of variance.
VARY=EYsq-mu^2
% Note that Var(Y) is also defined as
% integral -inf to inf (y-mu).^2 f_y(y) dy.
% We can create the definite integral of this equation and evaluate it
% to produce the variance
syms y, mu; mu=50;
VARY2=eval(int(sym('(y-mu)^2/100'),y,0,100))
% VARY2 should also be 2500/3
format
% using Theorem 3.13.12, in Larsen & Marx (2001, p. 223)
% Calculate the variance of all 100 transactions (multiply the variance by
% 100). Note, this is for 100 indepedently distributed random variates,
% where the positive and negative errors among transactions are not
% correlated.
Varall=100*VARY;
% If you wanted to convert the variance, calculated here in dollars to the
% variance in pennies, you'd have to use Theorem 3.13.1, p. 222.  It would
% the variance of Varall in pennies would be a=100;Var(a*Varall)=a^2 * Varall
%                                                            =10000 * Varall
% calculate the standard deviation;
sigmaall=sqrt(Varall);
% Calculate the Z statistic
Z=zeros(2,1);
Z(1)=(-500-0)/sigmaall;
Z(2)=(500-0)/sigmaall

% Use the cumulative normal probability function to find the probability of
% being more than $5 dollars off after 100 transactions:
P=1-(normcdf(Z(2))-normcdf(Z(1)))
fprintf('The probability of being more than $500 off after 100\n')
fprintf('transactions in which each transaction is rounded to the \n')
fprintf('nearest $100 is %6.4f\n',P)