function odds
% Calculate odds for n dice for various outcomes.
n = 5
tot_num_throws = 6^n
if tot_num_throws > 10^6
  return  % Too many possibilities.
end

% Create all possible n dice outcomes.
% Each row of matrix is one possible dice throw.
combs = (1:6)';
for index = 2:n
  % The add_combs function takes a set of vectors such as [1 1], [1 2], ...
  %  and appends every possible next element to every vector:
  %  [1 1 1], [1 1 2], ..., [1 2 5], [1 2 6], ..., .
  combs = add_combs(combs,1:6);
end

% Take off semicolons for diagnostics:
combs = combs;
size(combs);

% For each dice throw, count how many 1's, 2's, ..., n's
num_num_matrix = [];
for dice_num = 1:6
  % Create an empty array of zeros the same size as all possible dice throws.
  position_dice_num = zeros(size(combs,1),6);

  % Put a 1 in the matrix wherever a dice shows number = dice_num.
  position_dice_num(find(combs == dice_num)) = 1;

  % Count how many dice showed dice_num in each roll.
  num_num = sum(position_dice_num')';

  % num_num_matrix contains (for each possible shake) [#1's, #2's, ..., #6's].
  num_num_matrix = [num_num_matrix num_num];
end

num_num_matrix = num_num_matrix;

% Calculate prob of n of a kind.
for num_of_kind = 1:n
  % Calcualte prob of num_of_kind.
  num_of_kind = num_of_kind;

  % Count how many times num_of_kind_occurs.
  % This counts two pairs in one roll as two separate pairs.
  % This also counts more than one m-of-a-kind in one roll as separate
  % m-of-a-kind's.  This is only a problem for 2*m < n.
  tot_num_of_kind(num_of_kind) ...
    = length(find(num_num_matrix == num_of_kind));

  % Create matrix with 1's where num_of_kind occurs.
  % The column is the number showing on the dice.
  % The row is each possible roll.
  num_kind_matrix = (num_num_matrix ==  num_of_kind);

  % Count how many num_of_kind are in each shake.
  % Sum each row.  Result is one-column matrix, one row for each possible roll.
  rep_num_of_kind = sum(num_kind_matrix')';

  % Count how many times num_of_kind_occurs once, twice, ..., in each roll.
  for m = 1 : floor(n/num_of_kind)
    % Count how many rolls have m pairs, or m three-of-a-kind, etc.
    tot_num_of_kind(num_of_kind,m) ...
      = length(find(rep_num_of_kind == m));

    % Calculate probs of m num_of_kind in throw.
    prob_num_of_kind(num_of_kind,m) ...
    = tot_num_of_kind(num_of_kind,m)/tot_num_throws;

    % Print results for num_of_kind.
    fprintf('# times %d %d-of-kind occur in one roll = %6d      ', ...
            m,num_of_kind,tot_num_of_kind(num_of_kind,m))
    fprintf('prob of %d %d-of-kind occur in one roll = %6g\n', ...
            m,num_of_kind,prob_num_of_kind(num_of_kind,m))
  end
  fprintf('\n')

end

%disp('tot 1-of-kind, 2-of-kind, ...')
%tot_num_of_kind = tot_num_of_kind
%prob_num_of_kind = prob_num_of_kind

return

function comb_matrix = add_combs(prev_comb_matrix,vec)
% Add each element of vec to end of each row of prev_comb_matrix.
% Creates combinations of existing rows with elements of vec on the end.
rep_comb = repmat(prev_comb_matrix,length(vec),1);
new_col_comb = sort(repmat(vec',size(prev_comb_matrix,1),1));
comb_matrix = [new_col_comb, rep_comb];
return
