% Demo: List all multinomial count vectors
% Three categories ABC with probabilities .5 .3 .2
% Sequence length 10
%
function multidemo
  n = 10;
  pqr = [.5 .3 .2];
  
  % Initialize list of count vectors and probabilities
  L = [];
  
  % Counts are na,nb,nc
  for na=0:10
    for nb=0:10-na
      nc = 10-na-nb;
      counts = [na nb nc];
      
      % Multinomial likelihood
      lik = multicoeff(n,counts) * prod(pqr.^counts);
      
      L = [L; lik counts];
    end
  end
  
  % Sort by decreasing likelihood
  L = flipud(sortrows(L));
  nrows = size(L,1);
  
  fprintf('nrows = %d\n', nrows);
  
  
  for i=1:nrows
    if L(i,1) < 1e-3
      fprintf('(%d,%d,%d) %.6f\n', L(i,[2:4 1]));
    else
      fprintf('(%d,%d,%d) %.4f\n', L(i,[2:4 1]));
    end
  end
end


function result=multicoeff(n, kk)
  result = factorial(n);
  for i=1:length(kk)
    result = result / factorial(kk(i));
  end
end