%% Part 2 of the program that starts with "run_fitting.m" Read the program
%% notes for that file before running the program. 
%% 
%% Uses Maximum Likelihood methods for Weibull distribution(to be discussed in class).
%%
%% Written by Dhanoop Varghese, Sept. 2006. 
%%

%clear all; 
%close all; 
%clc;

%% Read in density function 
fden = load('fden.dat');
n = fden(:,1); 
f=fden(:,2);

%% Plot density function
figure(1);
loglog(n, f, 'ob', 'linewidth', 2);
set(gca, 'ylim', [1e-5 max(f)*5]);
hold on;

%% Find mean and standard deviation from read in data
mu0 = sum(f.*n)/sum(f);
sig0 = sqrt(sum(f.*(n-mu0).^2)/sum(f));
fprintf('Data statistics from input file:\n');
fprintf('mu0 = %e, sig0 = %e\n', mu0, sig0);

%% Read in distribution function
Fdist = load('Fdist.dat');
N = Fdist(:,1);

%% MLE fitting
fprintf('\n\nMLE fitting..\n');

%% Define grid for alpha and beta
alpha = mu0*10.^[-2:0.1:2];
beta = 10.^[-1:0.1:1.3];
fprintf('Grid boundaries for MLE:\n');
fprintf('alpha_min = %e, alpha_max = %e\n', min(alpha), max(alpha));
fprintf('beta_min = %e, beta_max = %e\n', min(beta), max(beta));

%% Evaluate joint probability
%% For each alpha
for i=1:size(alpha, 2)
    %% For each beta
    for j=1:size(beta, 2)
        prob(i, j)=0;
        %% Area under density function should be 1
        %% Evaluate normalizing factor 'norm' for each (alpha, beta)
        c1 = beta(j)/alpha(i);
        norm = quadl(@(x) c1*(x/alpha(i)).^(beta(j)-1).*exp(-(x/alpha(i)).^beta(j)), min(n), mu0+100*sig0);
        %fprintf('norm = %e\n', norm);
        for k=1:size(N, 1)
           if (norm <1e-3)
               prob(i, j) = -Inf;
               break;
           else
               peval=(c1/norm)*(N(k)/alpha(i))^(beta(j)-1)*exp(-(N(k)/alpha(i))^beta(j));
           end
           if (peval == 0)
                prob(i, j) = -Inf;
                break;
           else
                %fprintf('peval = %e\n', peval);
                prob(i, j) = prob(i, j)+ log(peval);
            end
        end
    end
end

% Plot joint probability for the defined (alpha, beta) grid
figure(2)
mesh(log(beta), log(alpha), prob);

% Find (<alpha>, <beta>) for which joint probability maximizes
[I, J] = find(prob==max(max(prob)));
fprintf('\nParameters estimated using MLE:\n');
fprintf('<alpha> = %e, <beta> = %e\n', alpha(I), beta(J));

% Check if <alpha>, <beta> lies on grid boundary
if ((I==1) || (I==size(alpha,2)))
    fprintf('\nWARN: <alpha>=%e found at grid boundary..\n', alpha(I))
    fprintf('WARN: The estimated value might be wrong!\n');
    fprintf('WARN: Please resize the alpha grid and re-run..\n');
    %return;
end

if ((J==1) || (J==size(beta,2)))
    fprintf('\nWARN: <beta>=%e found at grid boundary..\n', beta(J));
    fprintf('WARN: The estimated value might be wrong!\n');
    fprintf('WARN: Please resize the beta grid and re-run..\n');
    %return;
end

% Evaluate density function using estimated parameters
Nstep = 500;
Nmin = min(n);
Nmax = max(n);

c1 = beta(J)/alpha(I);
norm = quadl(@(x) c1*(x/alpha(I)).^(beta(J)-1).*exp(-(x/alpha(I)).^beta(J)), min(n), mu0+100*sig0);
        
for i=1:Nstep+1
    Neval(i) = Nmin + (Nmax-Nmin)*(i-1)/Nstep;
    feval(i)=(c1/norm)*(Neval(i)/alpha(I))^(beta(J)-1)*exp(-(Neval(i)/alpha(I))^beta(J));
end

figure(1);
loglog(Neval, feval, 'b', 'linewidth', 2);

%% LSQ fitting..
%clear alpha; clear beta; 
fprintf('\n\nLSQ Method..\n');

%% Define grid for mean and standard deviation
alpha = mu0*10.^[-2:0.1:2];
beta = 10.^[-1:0.1:1];
fprintf('Grid boundaries for LSQ:\n');
fprintf('alpha_min = %e, alpha_max = %e\n', min(alpha), max(alpha));
fprintf('beta_min = %e, beta_max = %e\n', min(beta), max(beta));

%% Evaluate LSQ error 
%% for each alpha
for i=1:size(alpha, 2)
    %% For each beta
    for j=1:size(beta, 2)
        error(i, j)=0;
        %% Area under density function should be 1
        %% Evaluate normalizing factor 'norm' for each (alpha, beta)
        c1 = beta(j)/alpha(i);
        norm = quadl(@(x) c1*(x/alpha(i)).^(beta(j)-1).*exp(-(x/alpha(i)).^beta(j)), min(n), mu0+100*sig0);
        for k=1:size(n, 1)
            if (f(k)==0)
                continue;
            end
            if (norm==0)
                error(i, j) = Inf;
                break;
            else
                peval=(c1/norm)*(n(k)/alpha(i))^(beta(j)-1)*exp(-(n(k)/alpha(i))^beta(j));
                error(i, j) = error(i, j)+ (f(k)-peval)^2/f(k)^2;
            end
        end
    end
end

%% Plot LSQ error for the defined grid
figure(3)
mesh(log(beta), log(alpha), log(error));

%% Find (<alpha>, <beta>) for which LSQ error minimizes
[I, J] = find(error==min(min(error)));
fprintf('\nParameters estimated using LSQ:\n');
fprintf('<alpha> = %e, <beta> = %e\n', alpha(I), beta(J));

% Check if <alpha>, <beta> lies on grid boundary
if ((I==1) || (I==size(alpha,2)))
    fprintf('\nWARN: <alpha>=%e found at grid boundary..\n', alpha(I))
    fprintf('WARN: The estimated value might be wrong!\n');
    fprintf('WARN: Please resize the alpha grid and re-run..\n');
    %return;
end

if ((J==1) || (J==size(beta,2)))
    fprintf('\nWARN: <beta>=%e found at grid boundary..\n', beta(J));
    fprintf('WARN: The estimated value might be wrong!\n');
    fprintf('WARN: Please resize the beta grid and re-run..\n');
    %return;
end

%% Evaluate density function for estimated (<alpha>, <beta>)
Nstep = 500;
Nmin = min(n);
Nmax = max(n);

c1 = beta(J)/alpha(I);
norm = quadl(@(x) c1*(x/alpha(I)).^(beta(J)-1).*exp(-(x/alpha(I)).^beta(J)), min(n), mu0+100*sig0);
        
for i=1:Nstep+1
    Neval(i) = Nmin + (Nmax-Nmin)*(i-1)/Nstep;
    feval(i)=(c1/norm)*(Neval(i)/alpha(I))^(beta(J)-1)*exp(-(Neval(i)/alpha(I))^beta(J));
end

figure(1);
loglog(Neval, feval, '--r', 'linewidth', 2);
set(gca, 'linewidth', 2, 'fontsize', 14);
xlabel('N -->');
ylabel('f -->');
legend('MC', 'MLE', 'LSQ');

figure(2);
set(gca, 'linewidth', 2, 'fontsize', 14);
xlabel('log(\sigma)');
ylabel('log(\mu)');
zlabel('log(Prob)');

figure(3);
set(gca, 'linewidth', 2, 'fontsize', 14);
xlabel('log(\sigma)');
ylabel('log(\mu)');
zlabel('LSQ Error');

figure(1);