Kasta tärning, räkna ut chans

n = 10; % Antal tärningar.
m = 6; % Antal sidor.
r = 6; % Antal kast.
N = nchoosek(n+m-1, m-1); % Antal tillstånd (exklusive starttillståndet).

% Generera matris där varje rad är ett tillstånd.
states = nan(N+1, m);
states(1, :) = 0;
states(2, :) = [n, zeros(1, m-1)];
for i = 3:N+1
  states(i, :) = states(i-1, :);
  k = find(states(i, :) != 0, 1);
  x = states(i, k);
  states(i, k) = 0;
  states(i, k+1) += 1;
  states(i, 1) = x-1;
endfor

% Generera övergångsmatris där p_ij är sannolikheten att gå från i till j.
P = nan(N+1, N+1);
P(:, 1) = 0;
% Första raden.
for j = 2:N+1
  p = factorial(n);
  for k = 1:m
    p /= factorial(states(j, k));
  endfor
  p /= m^n;
  P(1, j) = p;
endfor
% Resterande rader.
for i = 2:N+1
  [count, mode] = max(states(i, :));
  for j = 2:N+1
    if count == n
      % Kasta en tärning om alla är lika.
      if states(j, mode) >= n-1
        P(i, j) = 1/m;
      else
        P(i, j) = 0;
      endif
    elseif count == n-1
      % Kasta ingen tärning om alla utom en är lika.
      if j == i
        P(i, j) = 1;
      else
        P(i, j) = 0;
      endif
    else
      % Kasta de n-count tärningar som inte visar det mest förekommande
      % antalet prickar.
      if states(j, mode) < count
        P(i, j) = 0;
      else
        % Vektorn roll är tillståndet hos de kastade tärningarna.
        roll = states(j, :);
        roll(mode) -= count;
        p = factorial(n-count);
        for k = 1:m
          p /= factorial(roll(k));
        endfor
        p /= m^(n-count);
        P(i, j) = p;
      endif
    endif
  endfor
endfor

% Startfördelning.
w = [1, zeros(1, N)];

% Beräkna fördelning efter r kast.
for i = 1:r
  w = w*P;
endfor

% Summera sannolikheterna att få n-1 lika.
p = sum(w(any(states == n-1, 2)));
disp(p);