%% Temple university CIS 2033
% A modern Introduction to Probability and Statistics
% Chapter 5: Continuous random variables
%
% Author: Djordje Gligorijevic
% email: gligorijevic@temple.edu

%% Example 1 (a)
%
% This excercise demonstrates ways of creating a probability density 
% function using MATLAB's makedist and pdf functions.
%
% Plot probability density function for continuous random variables
% X ~ U(0, beta), beta = 0.5, 1, 2 
distnames = {'Uniform'};
colors = {'r','b','k','m','c'};
betas=[0.5, 1, 2];
xs = -1:0.01:3;
fig = figure;
for i=1:length(betas)
    uds = makedist(distnames{1},'lower',0,'upper',betas(i));
    ys = pdf(uds,xs);
    stairs(xs,ys,colors{i});
    hold on;
end
axis([-1 3 -0.1 3]);
legend('U(0, 0.5)','U(0, 1)', 'U(0, 2)');
saveas(fig,'D1a.eps','epsc');

% X ~ Par(alpha), alpha = 0.5, 1, 2
distnames = {'GeneralizedPareto'};
alphas = [0.5, 1, 2];
xs = 0:0.01:5;
fig = figure;
for i=1:length(alphas)
    k1 = 1/alphas(i);
    sigma1 = 1/alphas(i);    
    theta1 = 1;
    uds = makedist(distnames{1},'k',k1,'sigma',sigma1,'theta',theta1);
    ys = pdf(uds,xs);
    plot(xs,ys,colors{i});
    hold on;
end
axis([0 5 -0.2 2.5]);
legend('Par(0.5)','Par(1)', 'Par(2)');
saveas(fig,'D1b.eps','epsc');

%% Example 1 (b)
%
% This excercise does not requre creating a new probability function but 
% using MATLAB's functions.
%
% Plot probability density function for continuous random variables
% X ~ U(0, beta), beta = 0.5, 1, 2 
colors = {'r','b','k','m','c'};
betas=[0.5, 1, 2];
xs = -1:0.01:3;
fig = figure;
for i=1:length(betas)
    ys = unifpdf(xs,0,betas(i));
    stairs(xs,ys,colors{i});
    hold on;
end
axis([-1 3 -0.1 3]);
legend('U(0, 0.5)','U(0, 1)', 'U(0, 2)');
saveas(fig,'D1a.eps','epsc');

% X ~ Par(alpha), alpha = 0.5, 1, 2
alphas = [0.5, 1, 2];
xs = 0:0.01:5;
fig = figure;
for i=1:length(alphas)
    k1 = 1/alphas(i);
    sigma1 = 1/alphas(i);    
    theta1 = 1;
    ys = gppdf(xs,k1,sigma1,theta1);
    plot(xs,ys,colors{i});
    hold on;
end
axis([0 5 -0.2 2.5]);
legend('Par(0.5)','Par(1)', 'Par(2)');
saveas(fig,'D1b.eps','epsc');

%% Example 2 (a)
% 
% Plot distribution functions for continuous random variables. 
% X ~ U(0, 2)
distnames = {'Uniform'};
colors = {'b','c','k','m','r'};
betas=[2];
xs = -1:0.01:3;
fig = figure;
for i=1:length(betas)
    uds = makedist(distnames{1},'lower',0,'upper',betas(i));
    ys = cdf(uds,xs);
    stairs(xs,ys,colors{i});
    hold on;
end
median1 = uds.median;
plot(median1,0.5,'r*');
axis([-1 3 -0.2 1.2]);
legend('U(0, 2)','(q_{.5}, .5)');
saveas(fig,'D2a.eps','epsc');

% X ~ Par(2);
distnames = {'GeneralizedPareto'};
alphas = [2];
xs = 0:0.01:5;
fig = figure;
for i=1:length(alphas)
    k1 = 1/alphas(i);
    sigma1 = 1/alphas(i);    
    theta1 = 1;
    uds = makedist(distnames{1},'k',k1, 'sigma',sigma1,'theta',theta1);
    ys = cdf(uds,xs);
    plot(xs,ys,colors{i});
    hold on;
end
axis([0 5 -0.2 1.2]);
median1 = uds.median;
plot(median1,0.5,'r*');
legend('Par(2)','(q_{.5},.5)');
saveas(fig,'D2b.eps','epsc');

%% Example 2 (b)
% 
% Plot distribution functions for continuous random variables. 
% X ~ U(0, 2)
colors = {'b','c','k','m','r'};
betas=[2];
xs = -1:0.01:3;
fig = figure;
for i=1:length(betas)
    ys = unifcdf(xs,0,betas(i));
    stairs(xs,ys,colors{i});
    hold on;
end
median1 = 1; % You have to compute the median by yourself;
plot(median1,0.5,'r*');
axis([-1 3 -0.2 1.2]);
legend('U(0, 2)','(q_{.5}, .5)');
saveas(fig,'D2a.eps','epsc');

% X ~ Par(2);
alphas = [2];
xs = 0:0.01:5;
fig = figure;
for i=1:length(alphas)
    k1 = 1/alphas(i);
    sigma1 = 1/alphas(i);    
    theta1 = 1;
    ys = gpcdf(xs,k1,sigma1,theta1);
    plot(xs,ys,colors{i});
    hold on;
end
axis([0 5 -0.2 1.2]);
median1 = sqrt(2); % You have to compute the median by yourself;
plot(median1,0.5,'r*');
legend('Par(2)','(q_{.5},.5)');
saveas(fig,'D2b.eps','epsc');

%% Example 3
%
% Plots histogram of samples drawn from a Gaussian with 
% mean zero and std one N ~ (0,1)
% and plots the Gaussian with mean zero and std one
% the values of histogram are normalized so that the sum of the area of the
% histgoram bars (rectangles) is equal to one
% This way both plots are directly comparable, since the area under the
% Gaussian curve is also one.

x = -5:0.1:5;
g = normpdf(x,0,1);
sampleSize = 1000;

samples=randn(sampleSize,1);
% samples = normrnd(0,1,sampleSize,1);

[hs, hx] = hist(samples,x);
hs=(hs/sampleSize)*(101/10);

figure; 
bar(x,hs,'r'); %plots normalized histogram
hold on
plot(x,g,'LineWidth',3); %plots the Gaussian with mean zero and std one
title('Gaussian N(0,1) and the histogram of samples from N(0,1)')