function nma_Math504_HW1_part2()
%function nma_Math504_HW1_part2()
%
%This function simulates the pdf of the observed value from
%the following expeirment:
%
%pick a random number x from uniform[0,1], put this
%number in a box, and put twice the number in a second
%box. Next, pick one of these boxes by random, and look
%at the number inside. Call this y. If the number is smaller than
%one, then switch the box. Find the pdf of final y observed y, and
%find the estimate of the mean of y.
%Notice that  0<= y <= 2.

%by Nasser Abbasi, HW1, part 2, Math 504

clear all; close all;

seed  = 01010101;
rand('twister',seed);

% First do an initial estimate using simulation to estimate
% population mean and standard deviation, and then use these to
% obtain the needed sample size for the error level required

[s,xBar] = initialEstimate();
err = 0.01 * xBar;
sampleSize = round(((1.96*s)/err)^2);
fprintf('s=%f, xBar=%f\n',s,xBar);

%
%  C O N S T A N T S     and   P A R A M E T E R S
%
nBins = 50;
skip  = round(0.01*sampleSize); %for simulation, skip frames

%
%  I N I T I A L I Z A T I O N
%
y   = zeros(sampleSize,1);
figure(1);
set(0,'DefaultTextinterpreter','none');
h=title({'','',''});
axpos = get(gca,'pos');
extent = get(h,'extent');
set(gca,'pos',[axpos(1) axpos(2) axpos(3) axpos(4)-.20*extent(4)]);
set(h,'VerticalAlignment','Middle');

%
%  L O G I C
%
for i=1:sampleSize
    y(i) = makeAnObservation();
    if mod(i,skip)==0
        generateOneFrame(y(1:i),nBins);
    end
end

generateFinalResult(y,nBins);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [s,xbar]=initialEstimate()

sampleSize = 20000;
y = zeros(sampleSize,1);

for i = 1:sampleSize
    y(i) = makeAnObservation();
end

s    = std(y);
xbar = mean(y);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = makeAnObservation()

box = zeros(1,2);
box(1) = rand;
box(2) = 2*box(1);
boxSelected = (rand<.5)+1;       %pick a box by random

y = box(boxSelected);

if y <=1     %switch box if needed
    if boxSelected  == 1
        boxSelected = 2;
    else
        boxSelected = 1;
    end

    y = box(boxSelected);
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function generateOneFrame(y,nBins)

barWidth   = 2/nBins;
sampleSize = length(y);

firstLineTitle  = 'HW1 Mathematics 504, part(2) CSUF spring 2008';
secondLineTitle = 'simulation[$%d$],  $s=%6.5f$, $\\bar{X}=%9.8f$,  $\\mu=%5.4f$';

[n,x] = hist(y,nBins);
currentArea = barWidth*sum(n);
bar(x,n/currentArea,'y');  %relative frequency
ls=sprintf(secondLineTitle,sampleSize,std(y),mean(y),15/16);

title(char(firstLineTitle,ls),'fontsize',12,'interpreter','latex');
xlabel('x=final observed value');
ylabel('P(X=x)');
xtrue = [0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5];
line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--');
ylim([0,1]);
legend('simulation PDF','exact PDF');
drawnow;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function generateFinalResult(y,nBins)

barWidth = 2/nBins;
sampleSize = length(y);

firstLineTitle  = 'HW1 Mathematics 504, CSUF spring 2008';
secondLineTitle = 'simulation[$%d$],  $s=%6.5f$, $\\bar{X}=%9.8f$,  $\\mu=%5.4f$';

%
% Compare estimated mean with theortical mean
%
estimatedMean = mean(y);
standardError = 1.96*std(y)/sqrt(sampleSize);

lss=sprintf('95 perecent confidence interval for mean is (%6.5f ... %6.5f)',...
    estimatedMean-standardError,estimatedMean+standardError);

ls=sprintf(secondLineTitle,sampleSize,std(y),mean(y),15/16);

[n,x]=hist(y,nBins);
currentArea = barWidth*sum(n);
bar(x,n/currentArea,'y');  %relative frequency
title({firstLineTitle,ls,lss},'fontsize',12,'interpreter','latex');
xlabel('x=final observed value');
ylabel('P(X=x)');
xtrue=[0 .5 .5 1 1 2]; ytrue=[.75 .75 .25 .25 .5 .5];
line(xtrue,ytrue,'Color','r','LineWidth',2,'LineStyle','--');
ylim([0,1]);
legend('simulation PDF','exact PDF');
drawnow;

end