function D = lyngby_pde_binmix_der(N, A, B, P1, P2, lambda)

% lyngby_pde_binmix_der - Derivative of the binomial mixture
%
%       D = lyngby_pde_binmix_der(N, A, B, P1, P2, lambda)
% 
%       Output: D      Derivative as [ dp1 dp2 dlambda ]
%
%       Derivative of the log likelihood of the mixture of two
%       binomial distribution. 
%
%       See also LYNGBY, LYNGBY_PDE_BINMIX, LYNGBY_PDE_BINMIX_LL. 
%
% $Id: lyngby_pde_binmix_der.m,v 1.1 2004/03/08 21:41:25 fnielsen Exp $


    if nargin ~= 6
      error('Wrong number of input arguments')
    end

    N = N(:);
    A = A(:);
    B = B(:);
    
    P1 = min(max(P1,0),1);
    P2 = min(max(P2,0),1);
    lambda = min(max(lambda,0),1);

    
    % Number of replications
    M = length(N)-1;
    
    if 1

      % Derivative of lambda
      dlambda = sum( N .* (A-B) ./ ( (A-B)*lambda + B + realmin ));
      
      % Derivative of P1
      dp1 = 0;
      for n = 1:length(N)
	m = n - 1;
	dp1 = dp1 - N(n) * ...
	    lambda * P1^m *(1-P1)^(M-m) * (-m + P1*M) / ...
	    ( P1 * (-1+P1) * (-lambda*P1^m*(1-P1)^(M-m) - B(n) + ...
	    B(n)*lambda) + realmin ); 
      end
      
      % Derivative of P2
      dp2 = 0;
      for n = 1:length(N)
	m = n - 1;
	dp2 = dp2 + N(n) * ...
	    (1-P2)^(M-m) * P2^m * (-1+lambda) * (-m + P2*M) / ...
	    ( P2*(-1+P2) * ( -lambda * A(n) - P2^m * (1-P2)^(M-m) + P2^m ...
	    * (1-P2)^(M-m) * lambda) + realmin );
      end
      
      D = [ dp1 dp2 dlambda ];

    else
      % With link function
      
      P1 = 1 / ( 1+exp(-P1) );
      P2 = 1 / ( 1+exp(-P2) );
      lambda = 1 / ( 1+exp(-lambda) );
     
      
      e_y = (1-lambda) / lambda;
      dy = sum( N .* (e_y*(A-B))  ./ ( (A+e_y*B)*(1+e_y) + realmin) );
      
      dx1 = 0;
      for n = 1:length(N)
	m = n - 1;
	dx1 = dx1 - N(n) * ...
	    (m*(1-P1)/P1 - M + m) * (1-P1)^(M-m) * P1^m * lambda / ...
	    ( (lambda * P1^m * (1-P1)^(M-m) + B(n) - B(n)*lambda) * P1 );
      end

      dx2 = 0;
      for n = 1:length(N)
	m = n - 1;
	dx2 = dx2 + N(n) * ...
	    (m*(1-P2)/P2 - M + m) * (1-P2)^(M-m) * P2^m * (-1+lambda) ...
	    / ...
	    ( (-lambda * A(n) - P2^m * (1-P2)^(M-m) + P2^m * ...
	    (1-P2)^(M-m) * lambda) * P2 + realmin);
	  end
      
      D = [ dx1 dx2 dy ];
      
    end