%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    l_deriv.m     
%
%    Attila Vidacs
%
%    30/05/99
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  function y = l_deriv(h, METHOD)
%
%  Performs the minimization of the log-likelihood function
%  using its derivative.
%
%  Input :  h:  looks for the minimum near 'h'
%
%           METHOD = 'levinson' :  exact calculation of the inverse covariance
%                                  matrix using Levinson's algorithm.
%
%                  = 'app_lin'  :  using linear approximation for the inverse
%                                  of the covariant matrix .
%
%  Output:  The value of 'H' that minimizes the log-likelihood function

function y = l_deriv(h, METHOD)

global n t H x alpha

for i = 1:n
	for j = 1:n
		if( i ~= j )
			dS0(i,j) = log(t(i))*t(i)^(2*h) + log(t(j))*t(j)^(2*h) ...
				 - log(abs(t(i)-t(j)))*(abs(t(i)-t(j)))^(2*h);
		else
			dS0(i,j) = log(t(i))*t(i)^(2*h) + log(t(j))*t(j)^(2*h);

		end
	end
end

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 )';

switch (METHOD)
case('levinson')

	iG0 = lev(g0); % Inverse of G0 using Levinson's algorithm

case('app_lin') % approximate inverse with linear approximation of g(x)

	gamma = alpha^(1-h);
	iL(n,n) = 1/gamma;

	for i = 1:n-1
		iL(i,i+1) = -1;
		iL(i,i  ) = gamma + 1/gamma;
		iL(i+1,i) = -1;
	end

	iL(1,1) = 1/gamma;
	iG0 = iL * ( 1 / (1/gamma - gamma) );
end

iS0 = D*iG0*D;

diS0 = (-iS0*dS0*iS0);

u = (t'*iS0*t)*x;
v = (t'*iS0*x)*t;

y = trace(iS0*dS0)/n + ((v+u)'*(iS0*dS0*iS0)*(v+u))/((v+u)'*iS0*(v-u));

