clear
echo on
%This script illustrates how error propagates
%through a calculation when there is error in
%each of the variables.  Suppose we have two
%variables x and y which are normally distributed
%and independent.  Let's set these up so that x
%has a mean of 2 and y has a mean of 3 and both
%have standard deviations of unity:
n=1000;
x=randn(n,1)+2;
y=randn(n,1)+3;
pause
%We can plot the cumulative probability distribution
%of these two variables:
cumprob=([1:n]-.5)/n;
figure(1)
plot(sort(x),cumprob)
set(gca,'FontSize',18)
hold on
plot(sort(y),cumprob,'r')
hold off
xlabel('x , y')
ylabel('cumulative probability')
legend('x','y','Location','SouthEast')
pause
%Note that the slopes of the cumulative probabilities
%are the same, indicating that each variable has the
%same standard deviation.  

%Let's define a new variable z such that z = x + y.  
%What is its cumulative probability?
z=x+y;
figure(2)
plot((sort(z)-5),cumprob,'g')
set(gca,'FontSize',18)
hold on
plot((sort(x)-2),cumprob,'b')
plot((sort(y)-3),cumprob,'r')
hold off
xlabel('x,y,z')
legend('z - ave(z)','x - ave(x)','y - ave(y)','Location','SouthEast')
ylabel('cumulative probability')
title('cumulative distribution of x + y')
pause
%Here we have subtracted the mean off of x, y and z so
%that their slopes can be compared more readily.  Note
%that the slope of the cumulative probability of z is less
%than that of x and y, indicating a larger standard deviation.

%We can look at the distribution of z by calculating the mean
%and the standard deviation of z directly:
mean=sum(z)/n
%and the variance:
stdev = (sum((z-mean).^2)/(n-1))^.5
%Note that this standard deviation is very close to the 
%expected value of sqrt(2).  Subtraction would yield exactly
%the same standard deviation!
pause
%We can compare the actual cumulative probability distribution
%with a normal distribution with the same mean and variance:
figure(2)
hold on
plot(sort(z-5),(0.5 + erf((sort(z)-mean)/stdev/sqrt(2))/2),'k')
legend('z - ave(z)','x - ave(x)','y - ave(y)','fit','Location','SouthEast')
hold off
%which fits the data quite well.  
pause
%Now how about multiplication?  This is a little different.
%Let's define z = x * y and plot its cumulative probability:
z = x .* y;
figure(3)
plot((sort(z)-6),cumprob,'g')
set(gca,'FontSize',18)
hold on
plot((sort(x)-2),cumprob,'b')
plot((sort(y)-3),cumprob,'r')
hold off
xlabel('x,y,z')
legend('z - ave(z)','x - ave(x)','y - ave(y)','Location','SouthEast')
ylabel('cumulative probability')
title('cumulative distribution of x * y')
pause
%where again we have subtracted off the values  and .
%The first thing you should note is that z is -not- normally
%distributed!  This is because it is only normally distributed
%if the -fractional error- is small, that is sigma(x)/x <<1, and
%the same for y.  This is violated in this case, and it distorts
%the distribution of z.  We can calculate the mean and the
%standard deviation of z:
mean=sum(z)/n
%and the variance:
stdev = (sum((z-mean).^2)/(n-1))^.5


%We can compare the distribution of z to a normal distribution:
figure(3)
hold on
plot(sort(z-6),(0.5 + erf((sort(z)-mean)/stdev/sqrt(2))/2),'k')
legend('z - ave(z)','x - ave(x)','y - ave(y)','fit','Location','SouthEast')
hold off
%Note that this time the normal distribution is not as good a
%fit.  This is because we have violated the constraint that the
%error in x and y be small relative to the magnitude of x and y.
pause
%Suppose x and y were centered about 20 and 30 instead,
%but still had standard deviations of unity.  Now we get:
x=x+18;
y=y+27;
z=x .* y;
figure(4)
plot((sort(z)-600),cumprob,'g')
set(gca,'FontSize',18)
hold on
plot((sort(x)-20),cumprob,'b')
plot((sort(y)-30),cumprob,'r')
hold off
xlabel('x,y,z')
legend('z - ave(z)','x - ave(x)','y - ave(y)','Location','SouthEast')
ylabel('cumulative probability')
title('cumulative distribution of (x+18)*(y+27)')
pause
%z is now almost a normally distributed variable.  We can calculate
%the standard deviation:
mean=sum(z)/n
%and the variance:
stdev = (sum((z-mean).^2)/(n-1))^.5
pause
%We can graphically compare this to the observed probability
%distribution:
figure(4)
hold on
plot(sort(z-600),(0.5 + erf((sort(z)-mean)/stdev/sqrt(2))/2),'k')
legend('z - ave(z)','x - ave(x)','y - ave(y)','fit','Location','SouthEast')
hold off
pause
%This standard deviation could have been calculated from
%the formula:
sigmax=1;
sigmay=1;
stdev = ((sigmax/20)^2+(sigmay/30)^2)^.5*600
%which yields approximately the same value.
echo off