function [] = fit_cdf()
%
%###########################################################################
% Failure Trace Archive (FAT) scripts
% Program to fit distibution for availability and unavailability intervals 
% 
% Exponential
% Weibull
% Log-Normal
% Gamma
% Pareto
% Output: QQ and PP plots

% Author: Bahman Javadi
%        School of Computing, Engineering and Mathematics
%        Western Sydney University, Australia 
%        b.javadi@westernsydeny.edu.au
%        http://fta.scem.uws.edu.au/
%
%%###########################################################################


set(0,'DefaultLineLineWidth',3)
set(0,'DefaultAxesFontName','Helvetica')
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultTextFontName','Helvetica')
set(0,'DefaultTextFontSize',16)


load fta_data.mat;

%dat = event_availability_interval;
dat = event_unavailability_interval;



xdat = dat;
xdat(xdat<=0) = [];

%Maximum x and y for qq plot
mux = 100;

[f,x]= ecdf(xdat);
n = numel (xdat);
l = 1;

while l
 disp('Q-Q and P-P plot ');
 disp('1. Exponential');
 disp('2. Weibull');
 disp('3. Pareto');
 disp('4. Log-Normal');
 disp('5. Gamma');
 disp('6. Quit');
    
 q=input('Select (1-6): ');
  
if q==1

expMLE = expfit(xdat)
y = exprnd(expMLE,n,1);
z = expcdf(x,expMLE);

h1 = figure;
qqplot2(xdat,y);
axis([0 mux 0 mux]);
xlabel('Emprical Quantiles (hrs)');
ylabel('Fitted Quantiles (hrs)');

h2 = figure;
qqplot2(f,z);
axis([0 1 0 1]);
xlabel('Emprical Quantiles');
ylabel('Fitted Quantiles');
dist = 'exp';
end


if q==2
wblMLE = wblfit(xdat)
y = wblrnd(wblMLE(1),wblMLE(2),n,1);
z = wblcdf(x,wblMLE(1),wblMLE(2));

h1 = figure;
qqplot2(xdat,y);
axis([0 mux 0 mux]);
xlabel('Emprical Quantiles (hrs)');
ylabel('Fitted Quantiles (hrs)');
h2 = figure;
qqplot2(f,z)
axis([0 1 0 1]);
xlabel('Emprical Quantiles');
ylabel('Fitted Quantiles');
dist = 'wbl';
end

if q==3
gpMLE = gpfit(xdat)
y = gprnd(gpMLE(1),gpMLE(2),gpMLE(2),n,1);
z = gpcdf(x,gpMLE(1),gpMLE(2),gpMLE(2));

h1 = figure;
qqplot2(xdat,y);
axis([0 mux 0 mux]);
xlabel('Emprical Quantiles (hrs)');
ylabel('Fitted Quantiles (hrs)');
h2 = figure;
qqplot2(f,z)
axis([0 1 0 1]);
xlabel('Emprical Quantiles');
ylabel('Fitted Quantiles');
dist = 'gp';
end

if q==4
lognMLE = lognfit(xdat)
y = lognrnd(lognMLE(1),lognMLE(2),n,1);
z = logncdf(x,lognMLE(1),lognMLE(2));

h1 = figure;
qqplot2(xdat,y);

axis([0 mux 0 mux]);

xlabel('Emprical Quantiles (hrs)');
ylabel('Fitted Quantiles (hrs)');
h2 = figure;
qqplot2(f,z)
axis([0 1 0 1]);
xlabel('Emprical Quantiles');
ylabel('Fitted Quantiles');
dist = 'logn';
end

if q==5
    gamMLE = gamfit(xdat)

y = gamrnd(gamMLE(1), gamMLE(2),n,1);
z = gamcdf(x,gamMLE(1),gamMLE(2));

h1 = figure;
qqplot2(xdat,y);
axis([0 mux 0 mux]);
xlabel('Emprical Quantiles (hrs)');
ylabel('Fitted Quantiles (hrs)');
h2 = figure;
qqplot2(f,z);
axis([0 1 0 1]);
xlabel('Emprical Quantiles');
ylabel('Fitted Quantiles');
dist = 'gam';
end

if q==6
    l=0;
end

if q~=6

end
    end

