% CDA6530:  Project 1 Solution  (Fall 2010)
% Question 1: generate discrete r.v.

% Question 1(a):  Binomial distribution (40, 0.7)
sampleN = 100; n= 40; p = 0.7; P = zeros(n+1, 1);
for k=0:n,   % theoretic value of P(X=k)
 P(k+1) = nchoosek(n,k)*p^k * (1-p)^(n-k);  % remember in Matlab the index starts from 1
end
Xvalue = 0:n;   % remember that the r.v. X is from 0 to n 
num =length(P); SumP = zeros(num, 1);
for i=1:num,
    SumP(i) = sum(P(1:i));
end
U = rand(sampleN,1); Sample = zeros(sampleN,1);   %simulation value
for i=1:sampleN,
    for j=1:num,
        if U(i)< SumP(j),
            Sample(i) = Xvalue(j);
            break;
        end
    end
end
% draw probability mass function (pmf)
PMF = zeros(num,1);
for i=1:num,
    for j=1:sampleN,
        if Sample(j) == Xvalue(i),
            PMF(i)= PMF(i)+1;
        end
    end
end
PMF = PMF./sampleN; figure; plot(Xvalue, PMF, '-ob');
xlabel('Random Variable X'); ylabel('P(X)'); hold on;
plot(Xvalue, P, '--^r'); legend('Simulation', 'Analytical');
title('Probability Mass Function for Binomial Distribution');

%====================================================
% Question 1(b):  Geometric distribution p= 0.5
sampleN = 1000; p=0.5; 
n = 20;  % the maximum value we will simulate and show in our experiment
P = zeros(n, 1);  
for k=1:n,   % theoretic value of P(X=k)
 P(k) = (1-p)^(k-1) * p;  % remember in Matlab the index starts from 1
end
Xvalue = 1:n;  
num =length(P); SumP = zeros(num, 1);
for i=1:num,
    SumP(i) = sum(P(1:i));
end
U = rand(sampleN,1); Sample = zeros(sampleN,1);   %simulation value
for i=1:sampleN,
    for j=1:num,
        if U(i)< SumP(j),
            Sample(i) = Xvalue(j);
            break;
        end
    end
    if Sample(i) == 0,   % the r.v. is bigger than our bound defined in n
        disp('r.v. is bigger than the defined bound n! Geometric Distr.');
    end
end
% draw probability mass function (pmf)
PMF = zeros(num,1);
for i=1:num,
    for j=1:sampleN,
        if Sample(j) == Xvalue(i),
            PMF(i)= PMF(i)+1;
        end
    end
end
PMF = PMF./sampleN; figure; plot(Xvalue, PMF, '-ob');
xlabel('Random Variable X'); ylabel('P(X)'); hold on;
plot(Xvalue, P, '--^r'); legend('Simulation', 'Analytical');
title('Probability Mass Function of Gemetric Distribution');

%====================================================
% Question 1(c):  Poisson Distribution lambda= 6
sampleN = 1000; lambda=6;  
n = 20;  % the maximum value we will simulate and show in our experiment
P = zeros(n+1, 1); 
for k=0:n,   % theoretic value of P(X=k)
 P(k+1) = exp(-lambda) * lambda^k / factorial(k);  % remember in Matlab the index starts from 1
end
Xvalue = 0:n;  
num =length(P); SumP = zeros(num, 1);
for i=1:num,
    SumP(i) = sum(P(1:i));
end
U = rand(sampleN,1); Sample = zeros(sampleN,1);   %simulation value
for i=1:sampleN,
    for j=1:num,
        if U(i)< SumP(j),
            Sample(i) = Xvalue(j);
            break;
        end
    end
    if Sample(i) == 0,   % the r.v. is bigger than our bound defined in n
        disp('r.v. is bigger than the defined bound n! Poission Dist.');
        % if your code display too many of this warning messsage, 
        % you need to increase the bound n. 
    end
end
% draw probability mass function (pmf)
PMF = zeros(num,1);
for i=1:num,
    for j=1:sampleN,
        if Sample(j) == Xvalue(i),
            PMF(i)= PMF(i)+1;
        end
    end
end
PMF = PMF./sampleN; figure; plot(Xvalue, PMF, '-ob');
xlabel('Random Variable X'); ylabel('P(X)'); hold on;
plot(Xvalue, P, '--^r'); legend('Simulation', 'Analytical');
title('Probability Mass Function of Poission Distribution');

%====================================================
%====================================================
% Question 2(a):  Exponential distribution, lambda= 1.2

sampleN = 1000; lambda=1.2;  
Xvalue = 0:0.2:5;  % define the drawing points on the figure
                   % X=5 is the upper bound for our figure to draw
n = length(Xvalue); P = zeros(n, 1); 
for k=1:n,   % theoretic value of P(X=k)
 P(k) = 1- exp(-lambda * Xvalue(k));  % remember in Matlab the index starts from 1
end
U = rand(sampleN,1);  Sample = - log(1- U)./lambda;   %simulation value
F = zeros(n,1);
for i = 1:n,
    for j=1:sampleN,
        if Sample(j)<= Xvalue(i),
            F(i) = F(i) + 1;
        end
    end
end
F = F./sampleN; figure; plot(Xvalue, F, '-ob'); xlabel('Random Variable X'); 
ylabel('F(X)'); hold on; plot(Xvalue, P, '--^r');
legend('Simulation', 'Analytical'); title('CDF of Exponential Distribution');

%====================================================
% Question 2(b):  Normal distribution, mu=3.1,  sigma^2 = 5

sampleN = 1000; mu=3.1; sigma = sqrt(5);  
Xvalue = -1:0.2:7;  % define the drawing points on the figure
n = length(Xvalue); P = zeros(n, 1); 
for k=1:n,   % theoretic value of P(X=k)
 P(k) = quad(@(x)myfun(x,mu, sigma), -10, Xvalue(k));  
end
Sample = zeros(sampleN, 1);   %simulation value
for i=1:2:sampleN-1,
    S = 2;
    while S >1,
      U = rand(2,1); V = 2*U-1;  S = V(1)*V(1) + V(2)*V(2);
    end
    Sample(i:i+1) = sqrt(-2*log(S)/S) .* V;
end
Sample = Sample .* sigma +  mu; 
F = zeros(n,1);
for i = 1:n,
    for j=1:sampleN,
        if Sample(j)<= Xvalue(i),
            F(i) = F(i) + 1;
        end
    end
end
F = F./sampleN; figure; plot(Xvalue, F, '-ob'); xlabel('Random Variable X'); 
ylabel('F(X)'); hold on; plot(Xvalue, P, '--^r');
legend('Simulation', 'Analytical'); title('CDF of Normal Distribution');

%----------------------%
%function y = myfun(x, mu, sigma)
% normal distribution function for integral computing
%    y = exp(-(x-mu).^2./(2*sigma^2)) ./ (sqrt(2*pi)*sigma);
%end
%----------------------%



