clear all, close all, clc

% Estimates the PDF with the hystogram of uniform and normal distributions.
% then computes the exponential distribution from a uniform one

UniformLimits = [0 1];

NTMean = 0; % Normal Mean
NTVariance = 1; % Normal Variance

ETMean = 1; % Exponential Mean

N = 1e4; % Number of Samples
NBins = 100;

Legenda = ['Simulata';'Teorica '];

%% Uniform

UnifRV = rand(1,N); % Uniform Random Variable between 0 and 1
UnifRV = UniformLimits(1) + (UniformLimits(2)-UniformLimits(1))*UnifRV; % Translated between the Limits

UTMean = mean(UniformLimits);
UTVariance = (UniformLimits(2) - UniformLimits(1))^2/12;

USampleMean = mean(UnifRV);
USampleVariance = var(UnifRV);

fprintf('Uniforme con %d campioni\n\nMedia Teorica = %g\nVarianza Teorica = %g\n\nMedia Campionaria = %g\nVarianza Campionaria = %g\n\n',...
    N, UTMean, UTVariance, USampleMean, USampleVariance);

figHandler = figure('units','normalized','position',[0 0 1 1]);
[Counter, Centers] = hist(UnifRV,NBins);
Counter = Counter/sum(Counter)/(Centers(2) - Centers(1));

bar(Centers, Counter), hold on, grid on
x = min(UniformLimits)-0.1:.001:max(UniformLimits)+0.1;
y = ((x>=min(UniformLimits))&(x<=max(UniformLimits)))/abs((UniformLimits(2) - UniformLimits(1)));
plot(x, y, 'LineWidth', 2)
axis([min(Centers)-.1, max(Centers)+.1, 0, max(Counter)*1.1])

xlabel('x')
ylabel('f_X(x)')

title('P.d.f. Uniforme teorica e simulata')

hleg = legend(Legenda);
hleg.FontSize = 12;

%% Normal

NormalRV = randn(1,N); % Uniform Random Variable between 0 and 1
NormalRV = NTMean + sqrt(NTVariance)*NormalRV; % Magnified and Translated

GSampleMean = mean(NormalRV);
GSampleVariance = var(NormalRV);

fprintf('Gaussiana con %d campioni\n\nMedia Teorica = %g\nVarianza Teorica = %g\n\nMedia Campionaria = %g\nVarianza Campionaria = %g\n\n',...
    N, NTMean, NTVariance, GSampleMean, GSampleVariance);

figHandler = figure('units','normalized','position',[0 0 1 1]);
[Counter, Centers] = hist(NormalRV,NBins);
Counter = Counter/sum(Counter)/(Centers(2) - Centers(1));

bar(Centers, Counter), hold on, grid on
x = NTMean-4*sqrt(NTVariance):.001:NTMean+4*sqrt(NTVariance);
y = normpdf(x, NTMean, sqrt(NTVariance));
plot(x, y, 'LineWidth', 2)
axis([NTMean-4*sqrt(NTVariance), NTMean+4*sqrt(NTVariance), 0, max(Counter)*1.1])
xlabel('x')
ylabel('f_X(x)')

title('P.d.f. Normale teorica e simulata')
hleg = legend(Legenda);
hleg.FontSize = 12;

%% Exponential from a Uniform

EXPQuantile = @(x) -ETMean*log(1-x); % Pointer at command defined funtion

ExpRV = rand(1,N); % Uniform Random Variable between 0 and 1

ExpRV = EXPQuantile(ExpRV);

ETVariance = ETMean^2;

ESampleMean = mean(ExpRV);
ESampleVariance = var(ExpRV);

fprintf('Esponenziale con %d campioni\n\nMedia Teorica = %g\nVarianza Teorica = %g\n\nMedia Campionaria = %g\nVarianza Campionaria = %g\n\n',...
    N, ETMean, ETVariance, ESampleMean, ESampleVariance);

figHandler = figure('units','normalized','position',[0 0 1 1]);
[Counter, Centers] = hist(ExpRV,NBins);
Counter = Counter/sum(Counter)/(Centers(2) - Centers(1));

bar(Centers, Counter), hold on, grid on
x = -ETMean/10:.001:5*ETMean;
y = (exp(-x/ETMean)/ETMean).* (x >= 0);
plot(x, y, 'LineWidth', 2)
axis([-ETMean/10, 5*ETMean, 0, max(Counter)*1.1])

xlabel('x')
ylabel('f_X(x)')

title('P.d.f. Esponenziale teorica e simulata con l''inversione della C.D.F.')
hleg = legend(Legenda);
hleg.FontSize = 12;