%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    mle.m     
%
%    Attila Vidacs
%
%    30/05/99
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This function performs the maximum likelihood estimation of the parameters 
%  of fractional Brownian traffic in the time domain. The process is defined 
%  as follows:
%
%        X(t) = m*t + sqrt(a)*Z(t),
%
%  where X(t) represents the amount of traffic arrived in (0,t), 'm' is the 
%  mean arrival rate, 'a' is a variance parameter, and Z(t) is a normalized 
%  fractional Brownian motion with Hurst parameter 'H'.

%
%  Functions called directly:
%    
%       alph.m    :  calculates the parameter 'alpha' for geometrical sampling
%                    to get the pre-defined measurement (time) resolution.
%
%       lev.m     :  calculates the inverse and determinant of a Toeplitz
%                    matrix using Levinson's algorithm
%
%       l_deriv.m :  performs the minimization of the log-likelihood function
%                    using its derivative.
%
%  Inputs:  No input arguments are given, all the necessary parameters must be 
%           set by editing the file 'mle.m' header part.
%
%  Outputs: The MLE estimates of the model parameters 'm', 'a' and 'H' are 
%           given, together with their standard deviations calculated from 'l' 
%           simulations.
%
%  Usage:   Edit the file 'mle.m' header, and run 'mle'.
%
%  Should you have any questions, don't hesitate to contact:
%    Attila Vidacs (vidacs@ttt-atm.ttt.bme.hu)
%

clear all
global alpha n m a H D t delta x

n = 32;			% number of sample points
l = 100;                % number of simulations
delta = 10^(-6);	% measurement (time) resolution

m = 1.0;		% mean input rate
a = 1.0;		% variance parameter
H = 0.8; 		% Hurst parameter

SAMPLING = 'geom';	% sampling method: 'geom' or 'lin' or 'cantor'  implemented: |  lin  |  geom  | cantor |
			%                                                            +-------+--------+--------+
EXACT      = 1; 	% exact MLE using iG0                                        |   X   |    X   |    X   |
LEVINSON   = 0; 	% exact MLE using Levinson algorithm                         |   X   |    X   |
APP_LIN    = 0; 	% MLE using linear approximation                             |   -   |    X   |
APP_FIT    = 0; 	% MLE using advanced approximation (parameter fitting)       |   X   |    X   |
         p = 5;		% number of parameters in advanced approx. (PARN)
	
D_LEVINSON = 0;		% MLE using the derivative of the log-likelihood function    |   -   |    X   |   
D_APP_LIN  = 0;         %     ...with linear approximation                           |   -   |    X   |
       tol = 0.005;	% relative error tolerance for H (for fzero())
     trace = 0;		% fzero() displays information

SEED       = 0;		% SEED value for the random number generator
TRACE      = 1;		% '= 0' : compressed output

NOISE      = 0;		% '= 1': Gaussian noise added to the samples
     sigma = 0.01;	%        ...with variance sigma^2

disp('MLE estimation')
disp('--------------')
disp(sprintf('Sampling method           = %s', SAMPLING))
disp(sprintf('Number of samples (n)     = %d', n))
disp(sprintf('Number of simulations (l) = %d', l))
disp(sprintf('Model parameters: m       = %f', m))
disp(sprintf('                  a       = %f', a))
disp(sprintf('                  H       = %f', H))

simulation_start = clock; 

%%%%%%%%%%%%%%%%%%%%%%
%   generating fBm   %
%%%%%%%%%%%%%%%%%%%%%%
if ( TRACE ) disp('generating fBm...'), end;

switch ( SAMPLING )

  case 'geom'   % geometrical sampling: The sample points on (0,1) form a geometric sequence: t(i) = alpha^(i-1)
  
	alpha = min( fzero('alph', 0.9), delta^(1/(n-1)) );
	t = alpha.^(0:n-1)';
	disp(sprintf('alpha                     = %f', alpha))
	
  case 'lin'    % linear sampling: The sample points are evenly dsirtibuted on the sampled interval

	t = [1/n : 1/n : 1]';
	
  case 'cantor' % fractal-like sampling: The sample points are selected according to a (finite-resolution) Cantor set
  
	k = floor(log2(n));
	scale = delta^(-1/k)
	t = [1/scale; 1-1/scale];

	for i = 1:k-1
		t = t./scale;
		t(size(t,1)+1:2*size(t,1)) = 1-flipud(t);
	end	

	if (n ~= 2^k)	
		n = 2^k;
		disp(sprintf('n was not a power of 2 and has been changed!!!  n = %d', n))
	end
end

%  S0 is the covariance matrix of X, where X is the vector of samples
%  SB is the covariance matrix of a Brownian motion (for noise)

for i = 1:n
	for j = 1:n
		S0(i,j) = ( t(i)^(2*H) + t(j)^(2*H) - abs(t(i) - t(j))^(2*H) )/2;
		SB(i,j) = min(t(i), t(j));
	end
end

randn('seed', SEED);
Z = sqrtm(S0) * normrnd(0,1,[n,l]);  % matrix of samples of normalized fractional Brownian motion (n x l)


if ( NOISE )
	B = sqrtm(SB) * normrnd(0,1,[n,l]);  % samples of additive noise

%	for i = 1:l  X(:,i) = m.*t + sqrt(a).*Z(:,i) + normrnd(0, sigma, [n,1]); end % Gaussian white noise
	for i = 1:l  X(:,i) = m.*t + sqrt(a).*Z(:,i) + sigma.*B(:,i); end            % Brownian motion as noise
else
	for i = 1:l  X(:,i) = m.*t + sqrt(a).*Z(:,i); end  % matrix of samples of fractional Brownian traffic (n x l)
end

%%%%%%%%%%% 
%   MLE   %
%%%%%%%%%%% 
if ( TRACE  ) disp('MLE estimation...'), end;

ndx = 0;
tic;     % timer set

for ih = 10:80

	h = 0.5 + ih/200;
	if ( TRACE ) & ( EXACT | LEVINSON | APP_LIN | APP_FIT ) h, end;
	ndx = ndx + 1;
	h_list(ndx) = h;

	D = diag(t.^(-h));  % descaling matrix

	if ( EXACT  ) | ( LEVINSON  )
		switch ( SAMPLING )
		  case 'geom'
			g0 = ( alpha.^(-h*[0:n-1]) .* ( 1 + alpha.^(2*h*[0:n-1]) - (1-alpha.^([0:n-1])).^(2*h) )/2 )';

			if ( EXACT  )
				G0 = toeplitz(g0);  % covariance matrix of descaled samples
				iS0 = D*inv(G0)*D;
				log_det_S0 = h*n*(n-1)*log(alpha) + log(det(G0));
				
			elseif ( LEVINSON  )
				[iC,log_detC] = lev(g0);  % inverse covariance matrix and its determinant of descaled samples
				iC0 = D*iC*D; 
				log_det_C0 = h*n*(n-1)*log(alpha) + log_detC;
			end
			
		  case 'lin'
		  	s0 = zeros(n,1);
			for d = 0:n-1
				s0(d+1) = ( (abs(d-1)/n)^(2*h) + ((d+1)/n)^(2*h) - 2*(d/n)^(2*h) )/2;
			end
			
			if ( EXACT  )
				S0 = toeplitz(s0);  % covariance matrix of samples
				iS0 = inv(S0);
				log_det_S0 = log(det(S0));
				
			elseif ( LEVINSON  )
				[iC0, log_det_C0] = lev(s0);  % inverse covariance matrix and its determinant of samples
			end
		
		case 'cantor'
			for i = 1:n
				for j = 1:n
					S0(i,j) = ( t(i)^(2*h) + t(j)^(2*h) - abs(t(i) - t(j))^(2*h) )/2;  % covariance matrix
				end
			end
			iS0 = inv(S0);			
			log_det_S0 = log(det(S0));
		end
	end

	if (APP_FIT ) % approximate inverse with parameter fitting
	
		switch ( SAMPLING )
		  case 'geom'
		  	clear g0;
			g0 = alpha.^(-h*[0:(2*p)]) .* ( 1 + alpha.^(2*h*[0:(2*p)]) - (1-alpha.^([0:(2*p)])).^(2*h) )/2;
			F0 = toeplitz(g0);  % covariance matrix of descaled samples
			
		  case 'lin'
			clear c
			c = zeros(2*p+1,1);
			for d = 0:(2*p)
				c(d+1) = ( (abs(d-1)/n)^(2*h) + ((d+1)/n)^(2*h) - 2*(d/n)^(2*h) )/2;
			end

			F0 = toeplitz(c);  % covariance matrix of samples
		end
		
		% The inverse of the covariance matrix is approximated by a band matrix
		% (same for the linear and geometrical sampling!)

		B = zeros(2*p-1, 1);
		B(p,1) = 1;
	
		res = F0(1:2*p-1, 1:2*p-1)\B;

		iAp = zeros(n+p,n+p);
	
		for i = 1:n
			for j = 1:p
				iAp(i, i+j-1) = res(j+p-1);
				iAp(i+j-1, i) = res(j+p-1);
			end
		end
	
		iA = iAp(1:n,1:n);

		% corner correction: Elements at he upper-left and lower-right corners are corrected

		B = zeros(p-1, 1);
		B(1,1) = 1;

		res = F0(1:p-1,1:p-1) \ (B - iA(p,1)*F0(1:p-1,p));

		for i = 1:p-1
			iA(1,i) = res(i);
			iA(i,1) = res(i);
			iA(n,n-i+1) = res(i);
			iA(n-i+1,n) = res(i);
		end

		b_p = ( 1 - F0(2,1:p+1)*iA(1:p+1,2) ) / F0(2,2) + iA(2,2);
	
		iA(2,2) = b_p;
		iA(n-1,n-1) = b_p;

		% calculating the determinant of the approximate inverse
		
		switch ( SAMPLING )
		  case 'geom'
	
			g0 = alpha.^(-h*[0:4]) .* (1 + alpha.^(2*h*[0:4]) - (1-alpha.^([0:4])).^(2*h) )/2;
			F0 = toeplitz(g0);

			B = zeros(3, 1);
			B(2,1) = 1;
	
			res = F0(1:3, 1:3)\B;
	
			iA0 = D*iA*D;
			log_det_iA = (n-1)*log(-res(3));
			log_det_A0 = h*n*(n-1)*log(alpha) - log_det_iA;
		
		  case 'lin'  % The determinant is calculated for a simplified band-matrix only, which is co-diagonal!
		              % (the same approximation must be repeated first as before 
		              %  to calculate the elements of the co-diagonal matrix)

			B = zeros(3, 1);
			B(2,1) = 1;
	
			res = F0(1:3, 1:3)\B;

			iDp = zeros(n+2,n+2);
	
			for i = 1:n
				for j = 1:2
					iDp(i, i+j-1) = res(j+1);
					iDp(i+j-1, i) = res(j+1);
				end
			end

			iD = iDp(1:n,1:n);

			% corner correction

			B = 1;

			res = F0(1,1) \ (1 - iD(2,1)*F0(1,2));

			iD(1,1) = res;
			iD(1,1) = res;
			iD(n,n) = res;
			iD(n,n) = res;
	
			b_p = ( 1 - F0(2,1:3)*iD(1:3,2) ) / F0(2,2) + iD(2,2);
	
			iD(2,2) = b_p;
			iD(n-1,n-1) = b_p;
	
			% determinant

			clear K;
			K = zeros(n+1, n+1);
	
			for i = 1:n
				for j = 1:n
					K(i+1,j) = iD(i,j);
				end
			end
	
			K(1,1) = 1;
			K(n+1,n+1) = 1;
	
			US = zeros(n+1,1);
			B = zeros(n+1,1);
	
			B(1,1) = 1;
	
			US = K\B;

			S = US(n+1:n+1,1:1);
	
			cp_total = 0;
	
			for i = 1:n-1
				cp_total = cp_total + log(-(iD(i,i+1)));
			end
	
% 			...would be in the general case
%			log_det_iA =log((-1)^(n)*cp_total*det(S))

			iA0 = iA;
			log_det_iA =cp_total + log(-det(S));
			log_det_A0 = -log_det_iA;
		end
	end

	if (APP_LIN )  % approximate inverse with linear approximation of g(x)
		gamma = alpha^(1-h);
		c1 = gamma/(gamma^2-1);
		c2 = (1+gamma^2)/(1-gamma^2);
	end

%  The actual parameter estimation procedure:
%  Instead of minimizing the log-likelihood function for each sample set (resulting in 'l' minimizations), the log-likelihood
%  function is calculated for a range of H values (using sufficiently fine grid for H). Using this method, the determinant and the 
%  inverse of the covariance matrix need to be calculated only once for each trial value of H. The calculated matrices and 
%  determinants then can be used for all the 'l' sample sets.

	for k = 1:l  % the estimation is repeated 'l' times
		switch ( SAMPLING )
		
		  case 'geom'
			x = X(:,k);
			
		  case 'lin'  % for linear sampling, the process used is the _increment process of X_, which is stationary.
			x = zeros(n,1);
		
			x(1) = X(1,k);
			for i = 2:n
				x(i) = X(i,k) - X(i-1,k);
			end
			
			t = ones(n,1)./n;
			
		  case 'cantor'
		  	x = X(:,k);
		end;

		if (EXACT )	% exact formulas
			h_res_exact(ndx,k) = log_det_S0 + n*log(( x'*iS0*x - (( t'*iS0*x )^2 ) / (t'*iS0*t) ));

			m_res_exact(ndx,k) = ( t'*iS0*x ) / ( t'*iS0*t ) ;

			a_res_exact(ndx,k) = ( x'*iS0*x - (( t'*iS0*x )^2 ) / (t'*iS0*t)  )/(n-1);
		end

		if (LEVINSON )  % exact formulas using Levinson's algorithm (much faster and more stable!)
			h_res_levinson(ndx,k) = log_det_C0 + n*log(( x'*iC0*x - (( t'*iC0*x )^2 ) / (t'*iC0*t) ));

			m_res_levinson(ndx,k) = ( t'*iC0*x ) / ( t'*iC0*t ) ;

			a_res_levinson(ndx,k) = ( x'*iC0*x - (( t'*iC0*x )^2 ) / (t'*iC0*t) )/(n-1);
		end

		if (APP_LIN )% estimates using iL (linear approximation)
			xdrdx = 0;

			for i = 2:n-1
				xdrdx = xdrdx + 2*c1 * alpha^(-h*(2*i-1)) *x(i)*x(i+1) + c2 * alpha^(-2*h*(i-1)) * x(i)^2;
			end

			xdrdx = xdrdx + 2*c1 * x(1)*x(2)*alpha^(-h) + (x(1)^2)/(1-gamma^2)...
				+ (alpha^(-2*h*(n-1)))/(1-gamma^2) * x(n)^2;

			h_res_lin(ndx,k) = h*(n-1)*log(alpha) + ((n-1)/n) * log(1-gamma^2) + log(( xdrdx - x(1)^2 ));

			m_res_lin(ndx,k) = x(1);
	
			a_res_lin(ndx,k) = ( xdrdx - x(1)^2 )/(n-1);
		end

		if (APP_FIT ) % estimates using iA (advanced approximation)

			h_res_fit(ndx,k) = log_det_A0 + n*log(( x'*iA0*x - (( t'*iA0*x )^2 ) / (t'*iA0*t) ));
			
			m_res_fit(ndx,k) = ( t'*iA0*x ) / ( t'*iA0*t ) ;
		
			a_res_fit(ndx,k) = ( x'*iA0*x - (( t'*iA0*x )^2 ) / (t'*iA0*t)  )/(n-1);
		end
	end
end

% The graphs of the log-likelihood functions are plotted as a function of H

if ( EXACT     ) plot(h_list, h_res_exact,    'b-'); end;
if ( LEVINSON  ) plot(h_list, h_res_levinson, 'b-'); end;
if ( APP_LIN   ) plot(h_list, h_res_lin,      'b-'); end;
if ( APP_FIT   ) plot(h_list, h_res_fit,      'b-'); end;

% The actual minimization is done in the following

for k = 1:l

	if ( EXACT  )
		[value, index] = min(h_res_exact(:,k));
		h_est_exact(k) = h_list(index);
		m_est_exact(k) = m_res_exact(index,k);
		a_est_exact(k) = a_res_exact(index,k); 
	end

	if ( LEVINSON  )
		[value, index] = min(h_res_levinson(:,k));
		h_est_levinson(k) = h_list(index);
		m_est_levinson(k) = m_res_levinson(index,k);
		a_est_levinson(k) = a_res_levinson(index,k); 
	end

	if( APP_LIN  )  
		[value, index] = min(h_res_lin(:,k));
		h_est_lin(k) = h_list(index);
		m_est_lin(k) = m_res_lin(index,k);
		a_est_lin(k) = a_res_lin(index,k); 
	end

	if ( APP_FIT  )  
		[value, index] = min(h_res_fit(:,k));
		h_est_fit(k) = h_list(index);
		m_est_fit(k) = m_res_fit(index,k);
		a_est_fit(k) = a_res_fit(index,k); 
	end

% The following methods use the derivative of the log-likelihood function. The minimization is performed by 'fzero'.

	if ( D_LEVINSON  )
		x = X(:,k);
		h = fzero('l_deriv', [0.55, 0.95], tol, trace, 'levinson');
		h_est_d_levinson(k) = h;
		
		D = diag(t.^(-h));
		g0 = ( alpha.^(-h*[0:n-1]) .* ( 1 + alpha.^(2*h*[0:n-1]) - (1-alpha.^([0:n-1])).^(2*h) )/2 )';
		iC0 = D*lev(g0)*D;

		m_est_d_levinson(k) = ( t'*iC0*x ) / ( t'*iC0*t ) ;
		a_est_d_levinson(k) = ( x'*iC0*x - (( t'*iC0*x )^2 ) / (t'*iC0*t) )/(n-1);
	end
	
	if ( D_APP_LIN  )
		x = X(:,k);
		h = fzero('l_deriv', [0.55, 0.95], tol, trace, 'app_lin');
		h_est_d_lin(k) = h;

		gamma = alpha^(1-h);
		c1 = gamma/(gamma^2-1);
		c2 = (1+gamma^2)/(1-gamma^2);
		xdrdx = 0;

		for i = 2:n-1
			xdrdx = xdrdx + 2*c1 * alpha^(-h*(2*i-1)) *x(i)*x(i+1) + c2 * alpha^(-2*h*(i-1)) * x(i)^2;
		end

		xdrdx = xdrdx + 2*c1 * x(1)*x(2)*alpha^(-h) + (x(1)^2)/(1-gamma^2)...
			+ (alpha^(-2*h*(n-1)))/(1-gamma^2) * x(n)^2;

		m_est_d_lin(k) = x(1);
		a_est_d_lin(k) = ( xdrdx - x(1)^2 )/(n-1);
	end
	
	if ( TRACE & ( D_LEVINSON | D_APP_LIN ) )
		tmore = fix( (toc/k)*(l-k) );
		if (tmore > 60) disp(sprintf('%d (cca. %d min more...)', l-k, fix(tmore/60) )), end;
		if (tmore <= 60) disp(sprintf('%d (cca. %d sec more...)', l-k, tmore )), end;
	end
end

% The results are displayed. The MLE estimates of the three model parameters are printed, together their sample variance
% calculated from the 'l' simulations.

if ( EXACT  )
	disp(sprintf('H(exact) = %f, (std = %f)', mean(h_est_exact), std(h_est_exact)))
	disp(sprintf('m(exact) = %f, (std = %f)', mean(m_est_exact), std(m_est_exact)))
	disp(sprintf('a(exact) = %f, (std = %f)', mean(a_est_exact), std(a_est_exact)))
	disp(sprintf('---')) 
end

if ( LEVINSON  )
	disp(sprintf('H(levinson) = %f, (std = %f)', mean(h_est_levinson), std(h_est_levinson)))
	disp(sprintf('m(levinson) = %f, (std = %f)', mean(m_est_levinson), std(m_est_levinson)))
	disp(sprintf('a(levinson) = %f, (std = %f)', mean(a_est_levinson), std(a_est_levinson)))
	disp(sprintf('---')) 
end

if ( APP_LIN  )
%	disp(sprintf('H(lin) = %f, (std = %f)', mean(h_est_lin), std(h_est_lin)))
%
%       The 95% confidence interval is displayed in this case:

	HL = mean(h_est_lin) - 1.96*std(h_est_lin);  
	HR = mean(h_est_lin) + 1.96*std(h_est_lin);
	disp(sprintf('H(lin) = %f, (%f), [%f, %f]', mean(h_est_lin), std(h_est_lin), HL, HR))

	disp(sprintf('m(lin) = %f, (std = %f)', mean(m_est_lin), std(m_est_lin)))
	disp(sprintf('a(lin) = %f, (std = %f)', mean(a_est_lin), std(a_est_lin)))
	disp(sprintf('---')) 
end

if ( APP_FIT )
	disp(sprintf('H(fit) = %f, (std = %f)', mean(h_est_fit), std(h_est_fit)))
	disp(sprintf('m(fit) = %f, (std = %f)', mean(m_est_fit), std(m_est_fit)))
	disp(sprintf('a(fit) = %f, (std = %f)', mean(a_est_fit), std(a_est_fit)))
	disp(sprintf('---')) 
end

if ( D_LEVINSON  )
	disp(sprintf('H(d_levinson) = %f, (std = %f)', mean(h_est_d_levinson), std(h_est_d_levinson)))
	disp(sprintf('m(d_levinson) = %f, (std = %f)', mean(m_est_d_levinson), std(m_est_d_levinson)))
	disp(sprintf('a(d_levinson) = %f, (std = %f)', mean(a_est_d_levinson), std(a_est_d_levinson)))
	disp(sprintf('---')) 
end

if ( D_APP_LIN  )
	disp(sprintf('H(d_lin) = %f, (std = %f)', mean(h_est_d_lin), std(h_est_d_lin)))
	disp(sprintf('m(d_lin) = %f, (std = %f)', mean(m_est_d_lin), std(m_est_d_lin)))
	disp(sprintf('a(d_lin) = %f, (std = %f)', mean(a_est_d_lin), std(a_est_d_lin)))
	disp(sprintf('---')) 
end

simulation_end = clock;
ttime = fix(etime(simulation_end, simulation_start));
disp(sprintf('Total simulation time: %d minute(s) %d second(s)', fix(ttime/60), ttime - 60*fix(ttime/60) ))

