function [ll, A, B] = lyngby_pde_binmix_ll(N, P1, P2, lambda)

% lyngby_pde_binmix_ll  - Loglikelihood of mixture of binomials
%
%       [ll, A, B] = lyngby_pde_binmix_ll(N, P1, P2, lambda)
%
%       Computes the log-likelihood of the binomial mixture.
%
%       Ref: John Uebersax, MIXBIN.
%
%       Example: 
%         N = [ 13560 877 168 66 42 28 23 39 64];
%         x=0:0.01:1;
%         for n=1:length(x) 
%           ll(n) = lyngby_pde_binmix_ll(N, x(n), 0.01, 0.0195); 
%         end
%         figure, plot(x, ll)
%
%         y=0.001:0.001:0.03
%         for n=1:length(x) 
%           for m=1:length(y)
%             ll(m,n) = lyngby_pde_binmix_ll(N, x(n), y(m), 0.0195); 
%           end
%         end
%         figure, mesh(x, y, ll), xlabel('x'), ylabel('y')
%
%       See also LYNGBY, LYNGBY_PDE_BINMIX, LYNGBY_PDE_BINMIX_DER. 
%
% $Id: lyngby_pde_binmix_ll.m,v 1.1 2004/03/08 21:42:42 fnielsen Exp $


    if 1
      P1 = min(max(P1,0),1);
      P2 = min(max(P2,0),1);
      lambda = min(max(lambda,0),1);
    else
      % With link function 
      P1 = 1 / ( 1+exp(-P1) );
      P2 = 1 / ( 1+exp(-P2) );
      lambda = 1 / ( 1+exp(-lambda) );
    end
      
    M = length(N)-1;

    ll = 0;
    llnorm = 0;
    for n = 1:length(N)
      m = n - 1;
      A(n) = P1^m * (1-P1)^(M-m); 
      B(n) = P2^m * (1-P2)^(M-m);
      ll = ll + N(n) * log( lambda*A(n) + (1-lambda)*B(n) + realmin);
      
      % Normalization factor
      llnorm = llnorm + N(n) * log( factorial(M) / factorial(m) / ...
	  factorial(M-m) );
    end

    ll = ll + llnorm;







