% Written by Vanja Dukic, 2000.
%
% Please cite: 
% Dukic and Gatsonis, “Meta-analysis of Diagnostic Test Accuracy Assessment Studies with Varying Number of Thresholds”, Biometrics, 59(4):936-946
% Any questions or comments, please contact me: vanja.dukic@colorado.edu
%
% Prof. Vanja Dukic
% Department of Applied Mathematics
% University of Colorado-Boulder
% Box 526; Boulder, CO 80309-0526
% tel: 303-735-0919; fax: 303-492-4066


nzcum=cumsum(nz);

% maximize likelihood using Newton's method
epsilon=std(d2l(:))/1000;

iter=2; Lt=[]; Lt(1)=10; Lt(2)=L; 
while abs(Lt(iter)-Lt(iter-1)) > tol,
   iter=iter+1
   nt=nt+1;
   devold=dev;
   tb(:,nt)=tb(:,nt-1)- inv(d2l)*dl;

% likelihood at the current tb values
   Lk=zeros(1,K); devk=zeros(1,K);

   for k=1:K,
     eval([ 'y = y' num2str(k) ' ; ']);
     eval([ 'x = x' num2str(k) ' ; ']);

     if k==1, 
     theta = tb( 1 : nzcum(k),nt); 
     else 
     theta = tb( (nzcum(k-1)+1) : nzcum(k) , nt);  
     end;
     eval([ 'theta' num2str(k) ' = theta; '])
     beta = tb( (nzcum(K)+1) :(nzcum(K)+1)  ,nt);
     scale = tb( (nzcum(K)+2) : (nzcum(K)+2) ,nt);
     sc1 = tb(size(tb,1), nt);

     tbt=[]; tbt=[theta; beta; scale; sc1];
     eval([ 'tb' num2str(k) '(:,nt) = tbt; '])
     [Lk(k), devk(k)] = lik(y,x,x,tbt,pr(k));
   end;
     Lt(iter)=sum(Lk); dev=sum(devk);

    if (dev-devold)/(dl'*(tb(:,nt)-tb(:,(nt-1)))) < 0,
      epsilon=epsilon/decr;
   else,
      while (dev-devold)/(dl'*(tb(:,nt)-tb(:,(nt-1))   )) > 0,
         epsilon=epsilon*incr;
         if epsilon>1e+15,
            disp('epsilon too large');
            return;
         end;
         tb(:,nt)=tb(:,nt-1)-(d2l-epsilon*eye(size(d2l)))\dl;
      
for k=1:K,
    eval([ 'y = y' num2str(k) ' ; ']);
    eval([ 'x = x' num2str(k) ' ; ']);

     if k==1, 
     theta = tb( 1 : nzcum(k),nt); 
     else 
     theta = tb( (nzcum(k-1)+1) : nzcum(k) , nt);  
     end;
     eval([ 'theta' num2str(k) ' = theta; '])
     beta = tb( size(tb,1)-2 ,nt);
     scale = tb( size(tb,1)-1 ,nt);
     sc1 = tb(size(tb,1), nt);


     tbt=[]; tbt=[theta; beta; scale; sc1];
     eval([ 'tb' num2str(k) '(:,nt) = tbt; '])
    [Lk(k), devk(k)] = lik(y,x,x,tbt,pr(k));
end;
     Lt(iter)=sum(Lk); dev=sum(devk);
         disp('Epsilon'); disp(epsilon);
      end;
    end;

    


% derivatives at the current values
dlt=[]; d2lt=[]; dl=[]; d2l=[]; db=0; dbb=0; db2=[];

for k=1:K,
    eval([ 'y = y' num2str(k) ' ; ']);
    eval([ 'x = x' num2str(k) ' ; ']);

    dk=[]; d2k=[];
    eval([ 'tbt = tb' num2str(k) '(:,nt); '])
    [dk,d2k] = derivs(y,x,x,tbt, pr(k));
    eval([ 'dlk' num2str(k)  '  =  dk; ' ]);
    eval([ 'd2lk' num2str(k) '  =  d2k; ' ]);



% first derivatives

    dlt = [dlt;dk(1:nz(k))];		% d wrt cutpoints

    db = db + dk( (nz(k)+1) : (np(k)) ); % d wrt loc and scale 

% second derivatives

    zert1=[]; zert2=[]; zert3=[];
    zert1 = zeros( size(d2lt,1), yrange(k) ); 
    zert2 = zeros( yrange(k), size(d2lt,2) );
    zert3 = d2k(1:yrange(k),1:yrange(k));

    d2lt = [ [d2lt, zert1]; [zert2, zert3]  ];  % d2 wrt cut/cut


    zert4=[]; 
    zert4 = d2k( (yrange(k)+1): np(k) , 1:yrange(k) );

    db2 = [db2, zert4];	    % d2 wrt (beta, scale, sc1)/cut


    zert5=[]; 
    zert5 = d2k( (yrange(k)+1): np(k) ,  (nz(k)+1):np(k) );
    
    dbb = dbb + zert5;	% d2 wrt (b, sc, sc1)/(b, sc, sc1)



end;

dl  = [ dlt;  db ];
d2l = [ [d2lt, db2'] ; [db2 , dbb] ];

 
      disp('Iteration'); disp(iter-2);
      disp('Deviance'); disp(dev);
      format compact
      disp('First derivative'); disp(dl');
%     disp('Eigenvalues of second derivative'); disp(eig(d2l)');

end;






