% Start by looking at some of the functions in the library for 
%    generating PDFs and PMFs

help dfuni
  dfuni(rangelow, rangehigh, samples)  Generates a uniform pdf over the 
     given range.  samples sets the size of the numerical representation, 
     and is optional.  The default sample size is 1000.

f = dfuni(0,1);
plotf(f)

help dfnorm
  dfnorm(u, sigsq, samples)  Generates a Gaussian pdf with mean u and 
     variance sigma squared (sigsq).  samples sets the size of the numerical
     representation, and is optional.  The default sample size is 1000.

plotf( dfnorm(0,1) )

help mfbino
  mfbino(n, p)  generates a binomial pmf with parameters n, p

mb = mfbino(10, .3);
plotf(mb)

% Now look at the matrices that are produced by the above functions
%    Note the first row in each matrix contains identifiers that are
%    described in the helpfile.txt 

mb 

dfuni(0,1,10)

% The following example shows how to create one of these matrices for the 
%     density f(x) = 3x^2 between 0 and 1, and o elsewhere

x = [0: .001: 1];
fx = 3 * x .^2;
fnew = [ 0 1; x' fx' ];
plotf(fnew)

% Generating CDFs

Funi = getcf( f );
plotf(Funi)
Fnew = getcf(fnew);
plotf(Fnew)

% Calculating mean and variance 

mmean(f)

mmean(fnew)

var(f)

% Convolving PDFs

% Find PDF of sum of 2 independent uniform(0,1) random variables 
f2 = convpf(f, f);
plotf(f2)

% resize cuts out 0 part of matrix resulting from convolution
f2 = resize(f2,0,2);
plotf(f2)

% Now sum 3 

f3 = convpf( f2, f);
f3 = resize( f3, 0,3);
plotf(f3)


% Now add the sum of 4 independent U(0,1) random variables 

f4 = resize( convpf( f3, f), 0,4);
plotf(f4, 'g:')

f8 = resize( convpf(f4,f4), 0,8);
plotf(f8)

F8 = getcf(f8);
pconv = pval(F8,5) - pval(F8,4)


------------------------------------------------------------
%Now let's see what Mathematica can do
------------------------------------------------------------

%MATLAB is less accurate, but can go on and on and on
f16 = resize( convpf(f8,f8), 0,16);
f32 = resize( convpf(f16,f16), 0,32);
f48 = resize( convpf(f32,f16), 0,48);
plotf(f48,'m')

% Check the CLT result of #6.48

F48 = getcf(f48);
pconv = pval(F48,25) - pval(F48,22)

% See what CLT answer is for sum of 8 uniforms

% Simulation is another way to look at random variables

s1 = rand(8,1000);
colsum = sum(s1);
psim8 = sum(colsum>=4 & colsum<=5) / length(colsum)

s2 = rand(8,1000);
s3 = rand(8,1000);
s4 = rand(8,1000);
s5 = rand(8,1000);
s6 = rand(8,1000);
stot = s1+s2+s3+s4+s5+s6;
colsum = sum(stot);
psim = sum( colsum>=22 & colsum<=25 ) / length(colsum)

psim =

    0.5310


% rvgen is a function in the probability libarary that "draws" or 
%     samples from any PDF you have defined.  Actually it uses a CDF, 
%     but it will take either as input.

help rvgen
  rvgen(f, M, N, seed)  Generates a matrix of N samples of M experiments
     from a given probability or cumulative function.  Each row of the 
     matrix represent an experiment, and each column represents a sample.  
     The seed sets the random number generator seed.  So the resulting
     matrix is M by N.  Trials may be duplicated by using the same seed.
 
     rvgen would rather have a cumulative function.

suni = rvgen(f, 4, 1000, 123);

% See how uniform random numbers fill the unit square when you plot 
%     one row as the x-coordinate and another row as the y-coordinate.

plot( suni(1,:), suni(2,:) , '*' )

% Create a new simulation matrix which is the sum of two U(0,1) simulations

sim2 = [ suni(1,:) + suni(2,:); suni(3,:) + suni(4,:)];

% These simulations can all be plotted with the MATLAB function hist, or
%   if you know the PDF that the simulation is from (or should be from)
%   you can use the library function histo to plot the histogram and the
%   PDF simultaneously.

hist(sim2(1,:))
histo(sim2, f2, 20, 'g', 'b:')  
% we can use f2 as the PDF, since we figured that out above with convpf

% Finally, a magic trick.  Take the cube root of each number in the U(0,1)
%    simulated matrix.  Then plot its histogram and see that it has the
%    same PDF as our "new" random variable, i.e. f(x) = 3x^2 for x in [0,1].
%    You should learn in class soon why this works.  
% Hint: what is Fnew, and how does it relate to the cube root? 

plotf(fnew)
simnew = suni .^(1/3);
histo(simnew, fnew, 20, 'g', 'm:')

 8436197 flops.