%% ECS455 HW3
%% Prapun Suksompong, Ph.D.
%% ===============================================
close all
%% ===============================================
%%Q1
lambda = 30; T = 1000; n = 1000000; delta = T/n;
%Q1a
p1 = lambda*delta; fprintf('p1 = %5.3f \n',p1)

%Q1b
pp = binornd(1,p1,1,n); %Bernoulli(p) is binomial (1,p)
%pp = rand(1,n) < p1; %Generate Bernoulli via uniform

%Q1c
m = 1000;
expectedValue = lambda*T/m; K = 30;
fprintf('P[N_k = 30] = %5.3f \n',poisspdf(expectedValue,30))

%Q1d
N = sum(reshape(pp,n/m,m)); %stem(N);
A = sum(N == K);
fprintf('A = %d \n',A)
fprintf('Frequency of occurrence for {N_k = 30} = %5.3f \n',A/m)


figure
[f, x] = hist(N,0:2*K); plot(x,f/m,'kx'); xlim([0,2*K]);
hold on
plot(x,poisspdf(expectedValue,x),'ko')
title('The pmf of N'); xlabel('k'); grid on
legend('From histogram','From poisspdf formula')


%Q1e
W = diff(find(pp==1))*T/n; 
f = sum([W > 2/60])/length(W);
fprintf('Frequency of occurrence for {W_k > 2 mins} = %5.3f \n',f)

%Q1f
fprintf('P[W_k > 2 mins] = %5.3f \n',1-expcdf(1/lambda,2/60))
%exp(-1)

%% ===============================================
%%Q3
%Q3a
I = find(pp==1); V = length(I);
fprintf('V = %d \n',V)
%Q3b
mu = 60/2;
D = ceil(exprnd(1/mu,1,V)/delta);
figure
averageGeo = (1/mu)*n/T; %slots;
r = 1-mu*delta;
fprintf('D is a geometric r.v. with mean = %5.3f and parameter r = %5.3f \n',averageGeo,r)

K = floor(averageGeo);
[f, x] = hist(D,1:2*K); plot(x,f/V,'kx'); 
hold on
plot(x,geopdf(x-1,1-r),'ko') %Need to shft x by 1 because MATLAB geo r.v. start from 0.
xlim([0,2*K]); ylim([0 1.2*mu*delta]);
title('The pmf of D'); xlabel('d'); grid on
legend('From histogram','From geopdf formula')

%Q3c
m = 2;
st = I; et = I+D;
mm = [1]; %Vector mm keeps the index of calls that are occupying the channels
s = ones(1,V); %Vector s indicates whether the call request is successful
%The first call is always successful because no one is using the
%channel(s).
for i = 2:V %for each new arrival
    mm = mm(et(mm) >= st(i)); %eliminate the calls that have already ended
    if (length(mm) < m) %if there is still some available channels
        mm = [mm i]; s(i) = 1;
    else
        s(i) = 0;
    end
end
B = V-sum(s);
fprintf('B = %d \n',B)
fprintf('The proportion of call requests that were blocked is B/V = %5.3f \n',B/V)
fprintf('From Erlang B formula, the blocking probability is %5.3f \n',erlangb(m,lambda/mu))

%% ===============================================



