Question

Write a program (MATLAB) which generates exponential random variables, and use it to test the Central Limit Theorem as follows. For various values of n (say, 5, 10, 20, 50, 100, 1000), generate samples of the random variable Mn = 1/n \tiny \sum Xi where Xi are iid exponential random variables. A very simple method for generating exponential random variables is given on page 196 of the textbook (read that section). Plot the discretized CDF for Mn (an approximation of the CDF using discrete bins) for each n superimposed on the theoretical Gaussian CDF according to the Central Limit Theorem (or plot them next to each other). Calculate the sample mean and variance according to the material in the book and compare these to the theoretical values. Comment briefly on your results. Show all graphs and calculations. Include a copy of your code. Remember, you will need to generate a large number of samples of Mn for each n in order to get decent histograms. I suggest at least 100+ samples. 100 samples of M1000 require 100,000 exponential samples. You may choose any mean and variance you want for your exponential random numbers, but keep these fixed for all samples so that the results are comparable. Using the program written for the previous problem, repeat the whole experiment, but now choose a random mean and variance for each exponential sample. Thus, the Mn now will not be sums of identically distributed random variables. Show your results and comment. For this case, make sure that the distributions you use to choose your mean and variance at each sample are fixed throughout. Write the code in MATLAB.

Answer 1

% Distribution of Random variables
% Initialization
clear all;clc;
% Rate (Exponential Distribution)
lambda=4;
% Probability for bernoulli distributons
prob=0.6;
% maximum number of iterations
maxItr=10e3;
% Uniform Distribution (Linear Congruential Generator, LCG)
u=zeros(maxItr,1);
for i=1:maxItr
    
    % parameters of LCG
    t = datetime('now');
    z=minute(t)+randi([1,10]); c=second(t)+randi([1,10]);
    m=hour(t)+randi([1,10]); a=second(t)+randi([1,10]);
    
    u(i)=(mod((a*z+c),m))/m;
end
figure
histogram(u)
% Exponential Distribution
% Inverse transform method
exp=-log(u)/lambda;
figure
histogram(exp)
% Erlang Distribution
 sum=0;num_dist=10;
for i=1:num_dist
    
    % Generate exponential random numbers
    exp=exprnd(lambda,maxItr,1);
    sum=sum+exp;
    
end
erlang=sum/num_dist;
figure
histogram(erlang)
% Bernoulli Distribution
ber(u <= prob)=1;
figure
histogram(ber)
% Binomial Distribution
sum=0;num_dist=10;
for i=1:num_dist
    
    % Generate bernoulli random number
    u=rand(maxItr,1);
    ber(u <= prob)=1;
    ber(u > prob)=0;
    
    sum = sum + ber;
end
bimonial = sum;
figure
histogram(bimonial)
% Geometric Distribution
% Inverse transform method
geo=log(u)/log(1-prob);
figure
histogram(geo)
% Negative Binomial Distribution
sum=0;num_dist=10;
for i=1:num_dist
    
    % Generate geometric random number
    geo=geornd(prob,maxItr,1);
    
    sum=sum+geo;
end
bimonial=sum;
figure
histogram(bimonial)
% Normal Distribution (using Box Muller Method)
u1=zeros(maxItr,1);
u2=zeros(maxItr,1);
v1=zeros(maxItr,1);
v2=zeros(maxItr,1);
s=zeros(maxItr,1);

for i=1:maxItr
    
    while (1)
    
        %generating u1 = u~(0,1) 
        
        t = datetime('now');
        z=minute(t)+randi([1,10]); c=second(t)+randi([1,10]);
        m=hour(t)+randi([1,10]); a=second(t)+randi([1,10]);
        u1(i)=(mod((a*z+c),m))/m;
        
        %generating u1 = u~(0,1) 
        
        t = datetime('now');
        z=minute(t)+randi([1,10]); c=second(t)+randi([1,10]);
        m=hour(t)+randi([1,10]); a=second(t)+randi([1,10]);
        u2(i)=(mod((a*z+c),m))/m;
        
        % calculating 's'
        
        v1(i)=2*u1(i)-1;
        v2(i)=2*u2(i)-1;
        s(i)=v1(i)^2+v2(i)^2;
        
        if s(i)<=1
            break;
        end

    end

end

% This method generates two different Normal Variates
Normal_X=sqrt(-2*log(s)./s).*v1;
Normal_Y=sqrt(-2*log(s)./s).*v2;
figure
hold on
histogram(Normal_X,30)
histogram(Normal_Y,30)
hold off