restart:

with(Statistics):

Dice1 := RandomVariable(DiscreteUniform(1, 6));
Dice2 := RandomVariable(DiscreteUniform(1, 6));

g := (Dice1, Dice2) -> -4*(1-Heaviside(Dice1-3/2)*Heaviside(Dice2-3/2)) + abs(Dice2-Dice1)*Heaviside(Dice1-3/2)*Heaviside(Dice2-3/2)

Mean(g(Dice1, Dice2));
Variance(g(Dice1, Dice2));

S := Sample(g(Dice1, Dice2), 10^3):
Tally(S);
Histogram(S);

Nb_of_all_different_sequences := 6^6;

# also given by

combinat:-numbperm([seq(k$6, k=1..6)], 6)

# The probability to get the sequence S = [1, 2, 3, 4, 5, 6] is

P := 1 / Nb_of_all_different_sequences

with(Statistics):

# The "world" is made of the 46656 different sequences.
# An experiment consist of randomly picking one sequence from this "world".
# The occurrence of the sequence S can then modelled by a Bernoulli RV
# with parameter P.
# The number of sequences we must pick to obtain N occurences of S is then
# described by a begative binomial RV with parameters N and P

X := N -> RandomVariable(NegativeBinomial(N, P))

# Median of the number of sequences one must pick from the "world" before
# to get the sequence S

Quantile(X(1), 0.5, numeric);

restart:

with(GraphTheory):
with(SpecialGraphs):
with(combinat, powerset):

graph := 1;

if graph = 1 then
  G := OctahedronGraph();
  Cycles := CycleBasis(G);
  N := numelems(Cycles);
 
elif graph = 2 then
  G := CompleteGraph(4);
  Cycles := CycleBasis(G);
  N := numelems(Cycles);

elif graph = 3 then
  G := Graph({{1, 2}, {2, 3}, {3, 4}, {4, 1}, {2, 5}, {5, 6}, {6, 3}});
  Cycles := CycleBasis(G);
  N := numelems(Cycles);

elif graph = 4 then
  G := Graph({{1, 2}, {2, 3}, {3, 4}, {4, 1}, {5, 6}, {6, 7}, {7, 8}, {8, 5}, {3, 8}});
  Cycles := CycleBasis(G);
  N := numelems(Cycles);

elif graph = 5 then
  G := Graph({{1, 2}, {2, 3}, {3, 4}, {4, 1}, {5, 6}, {6, 7}, {3, 7}, {3, 5}});
  Cycles := CycleBasis(G);
  N := numelems(Cycles);

end if;


try
  DrawGraph(G, style=planar):
catch:
  DrawGraph(G, style=spring);
end try;

# close the list of vertices by adding the first vertex to its right

for n from 1 to N do
   C||n := [Cycles[n][], Cycles[n][1]];
end do:

# Create the sets of edges for each cycle

for n from 1 to N do
   S||n := map(m -> {(C||n)[m], (C||n)[m+1]}, {$1..numelems(C||n)-1});
end do;
 

S0 := {seq(S||n, n=1..N)}:
PS := powerset(S0):
NS := numelems(PS);

SD := proc(L::{list, set})
local N, n, sd:
N:= numelems(L):
if N = 1 then
   return L[1]
else
   if L[1] intersect L[2] <> {} then
     sd := ( L[1] union L[2] ) minus ( L[1] intersect L[2] );
     for n from 3 to N do
       if sd intersect L[n] <> {} then
         sd := ( sd union L[n] ) minus ( sd intersect L[n] )
       else
         return {}
       end if:
     end do:
     return eval(sd)
   else
     return {}
   end if:
end if:
end proc:
 

AllCycles := NULL:
for ns from 2 to NS do
   SymDiff      := SD([seq(PS[ns][k], k=1..nops(PS[ns]))]):
   vertices   := sort(ListTools:-Flatten(op~([SymDiff[]]))):
   NbMaxOccur := 0:
   for v in {vertices[]} do
     NbMaxOccur := max(NbMaxOccur, numelems([ListTools:-SearchAll(v, vertices)]))
   end do:
   if NbMaxOccur < 3 then
     AllCycles := AllCycles, SymDiff
   end if;
end do:
AllCycles := {AllCycles} union {seq(S||n, n=1..N)} minus {{}};
NAC := numelems(AllCycles);

# close the list of vertices by adding the first vertex to its right

for n from 1 to NAC do
   CAC||n := Graph(convert~(AllCycles[n], set));
end do:

p := proc(m)
  local n, CG:
  n := round(m):
  CG := CopyGraph(G):
  HighlightEdges(CG, Edges(CG), 'blue');
  HighlightEdges(CG, Edges(CAC||n), 'red');
  try
    DrawGraph(CG, style=planar):
  catch:
    DrawGraph(CG, style=spring);
  end try;
 
end proc:

plots:-animate(p, [m], m=1..NAC, frames=NAC)