close all;clear all
% make "fundamental" probability distribution function (PDF)
NP = 1000000;
X = (randn(NP,1)-2)+(randn(NP,1)-4).^2;

binNumP = 100;
figure;hist(X,binNumP);axis([-10, 80 -inf inf])

mu = mean(X); % mean of the fundamental pdf
sig2 = var(X,1); % variance of the fundamental pdf
sig = sqrt(sig2); % standard deviation of the fundamental pdf
[mu sig]
ll = 4*NP/binNumP; % not important
line([mu mu],[0 ll],'Color','c')
line([mu mu]+sig,[0 ll],'Color','m')
line([mu mu]-sig,[0 ll],'Color','m')

% myPercentileTest is an "anonymous" function which performs the following:
%
% function lookForZero = myPercentileTest(x,X,percentileSought)
% lookForZero = sum(X<=x)/length(X) - percentileSought;

myPercentileTest = @(x,X,percentileSought) (sum(X<=x)/length(X) - percentileSought);

%
p50 = fzero(myPercentileTest,mu,[],X,0.50); %  = median
% fzero finds the value of x which causes myPercentileTest to vanish
% i.e. it finds the value of x for which the CDF is equal to percentileSought.
% mu is the value to start searching from (as good a guess as any).
p05= fzero(myPercentileTest,mu,[],X,0.05); %  = 5th percentile
p95 = fzero(myPercentileTest,mu,[],X,0.95); %  = 95th percentils
[p50 p05 p95]
line([p50 p50],[0 ll],'Color','r')
line([p05 p05],[0 ll],'Color','r')
line([p95 p95],[0 ll],'Color','r')

% now draw 20 samples from fundamental distribution
NS = 20;
binNumS = 5;

figure
IndRand = ceil(NP*rand(NS,1));
x = X(IndRand); % 20 variables drawn at random
xBar = mean(x); %  estimate of mean
sx2 = var(x); % estimate of variance
sx = sqrt(sx2); % estimate of standard deviation
[xBar sx]
stdErrMean = sx/sqrt(NS); % estimate of error in estimate of mean
xp05 = fzero(myPercentileTest,xBar,[],x,0.05); % estimate of 5th percentile
xp50 = fzero(myPercentileTest,xBar,[],x,0.50); % estimate of median
xp95 = fzero(myPercentileTest,xBar,[],x,0.95); % estimate of 95th percentile
[xp50 xp05 xp95]
hist(x,binNumS);axis([-10, 80 -inf inf])
ll = 5*NS/binNumS; % not important
line([mu mu],[0 ll],'Color','c','LineWidth', 3)
line([xBar xBar],[0 ll],'Color','b','LineWidth', 2)
line([xBar xBar]+sx,[0 ll],'Color','m')
line([xBar xBar]-sx,[0 ll],'Color','m')
line([xBar xBar]+stdErrMean,[0 ll],'Color','b')
line([xBar xBar]-stdErrMean,[0 ll],'Color','b')
line([xp50 xp50],[0 ll],'Color','k','LineWidth', 2)

figure
for mm = 1:10 % and draw those 20 samples 10 times (e.g. 10 different experiments).
    IndRand = ceil(NP*rand(NS,1));
    x = X(IndRand); % 20 variables drawn at random
    xBar = mean(x); %  estimate of mean
    sx2 = var(x); % estimate of variance
    sx = sqrt(sx2); % estimate of standard deviation
    [xBar sx]
    stdErrMean = sx/sqrt(NS); % estimate of error in estimate of mean
 	xp05 = fzero(myPercentileTest,xBar,[],x,0.05); % estimate of median
	xp50 = fzero(myPercentileTest,xBar,[],x,0.50); % estimate of 5th percentile
	xp95 = fzero(myPercentileTest,xBar,[],x,0.95); % estimate of 95th percentile
    [xp50 xp05 xp95]
    subplot(10,1,mm);hist(x,binNumS);axis([-10, 80 -inf inf])
    ll = 5*NS/binNumS; % not important
	line([mu mu],[0 ll],'Color','c','LineWidth', 3)
    line([xBar xBar],[0 ll],'Color','b','LineWidth', 2)
    line([xBar xBar]+sx,[0 ll],'Color','m')
    line([xBar xBar]-sx,[0 ll],'Color','m')
    line([xBar xBar]+stdErrMean,[0 ll],'Color','b')
    line([xBar xBar]-stdErrMean,[0 ll],'Color','b')
	line([xp50 xp50],[0 ll],'Color','k','LineWidth', 2)
end
