% ################################################################
% Supporting file for top level program percolation.m.
% ###############################################################
% This program reads in the data from 'weibull.dat' file 
% and estimate alpha, beta using LSQ and MLE.
% #################################################################

format long;

% Initialize variables
maxerror = 1e-6;

% Read in data from input file
weibull=load('weibull.dat');
X = weibull(:,1);
F = weibull(:,2);
param = polyfit(log(X), log(-log(1-F)), 1);
fprintf('================\nLSQ Fitting\n===============\n');
fprintf('Beta from LSQ = %e\n', param(1));
fprintf('Alpha from LSQ = %e\n\n', exp(-param(2)/param(1)));
n = size(X, 1);

% Find bmax, amax for which L is maximized
% Start with an initial guess
bmax = param(1);
delta = bmax;
while abs(delta/bmax) > maxerror
    delta = -f(X, bmax, n)/fprime(X, bmax, n);
    bmax = bmax + delta;
end
amax = (sum(X.^bmax)/n)^(1/bmax);
fprintf('================\nMLE Fitting\n===============\n');
fprintf('bmax = %e\n', bmax);
fprintf('amax = %e\n', amax);

% To correct for the systematic error in bmax, amax
% b11 and a11 corresponding to alpha=beta=1 
% Compute the Y vector corresponding to F values
Y = (-log(1-F));

% Find b11 and a11
b11 = 1;
delta = b11;
while abs(delta/b11) > 1e-6
    delta = -f(Y, b11, n)/fprime(Y, b11, n);
    b11 = b11 + delta;
end
a11 = (sum(Y.^b11)/n)^(1/b11);
fprintf('b11 = %e\n', b11);
fprintf('a11 = %e\n', a11);

% Compute back alpha and beta
beta = bmax/b11;
alpha = amax/(a11^(b11/bmax));
fprintf('\nBeta from MLE = %e\n', beta);
fprintf('Alpha from MLE = %e\n', alpha);

% Evaluate Weibull distribution using the estimated parameters
Xmin = min(X);
Xmax = max(X);
Xstep = (Xmax-Xmin)/100;
Xeval = [Xmin:Xstep:Xmax];
Feval = 1-exp(-(Xeval/alpha).^beta);

figure;
loglog(X, -log(1-F), 'rd', 'linewidth', 2);
hold on;
loglog(Xeval, -log(1-Feval), 'r', 'linewidth', 2);
set(gca, 'linewidth', 2, 'fontsize', 16);
xlabel('T_B_D (s)');
ylabel('-log(1-F)');

% Compute variance for alpha, beta parameters
% varb = sum(X.^beta).^2/(sum(X.^beta)*sum(beta*X.^beta.*log(X).^2+X.^beta.*log(X))-sum(beta*X.^beta.*log(X))*sum(X.^beta.*log(X)));
% vara = 1/(beta*alpha^beta-1);

varb = 1/(length(X)/beta^2+sum((X/alpha).^beta.*log(X/alpha).*log(X/alpha)));
vara = 1/(beta*length(X)/alpha^2+sum(X.^beta)*beta*(beta+1)/alpha^(beta+2));

fprintf('\nVariance (Beta) = %e; Variance (Alpha) = %d\n', varb, vara);
fprintf('std (Beta) = %e; std (Alpha) = %d\n', sqrt(varb), sqrt(vara));