# Fast creation of state transitions and probabilities?

Normally I'd wait until I'd fleshed/thought things out further but I find this interesting enough that I think others might.

Take an ordered list of $N$ slots, integers from $[0,2]$. This represents the counts for $N$ categories. The actual categories/order matter not - we only care that we know we have $\{m_{1},m_{2},...,m_{n}\}$ hits in the categories, the result is kept ordered for convenience.

At each iteration of the experiment, we select one of the $N$ slots to be incremented (uniformly randomly). After the increment, we have:

1. No slot has $>2$ hit count, the iteration is complete.
2. Some slot ends up incremented to $3$ hits. That slot's hit count is decremented by $2$, and a slot is then selected to be incremented (again, uniformly randomly). We repeat case $2$ so long as any slot exceeds a hit count of 2.

I'm interested in a function that when fed the initial state returns a list of the initial and ending state of the list along with the probability associated with that path.

For example, given the initial state $\{0,1,2,2\}$ the output will be something like

{{{0, 1, 2, 2} -> {0, 1, 2, 2}, 1/4},

{{0, 1, 2, 2} -> {0, 1, 1, 2}, 3/32},

{{0, 1, 2, 2} -> {1, 1, 1, 1}, 1/32},

{{0, 1, 2, 2} -> {1, 1, 1, 2}, 1/8},

{{0, 1, 2, 2} -> {0, 2, 2, 2}, 1/4},

{{0, 1, 2, 2} -> {1, 1, 2, 2}, 1/4}}

A very quick-and-dirty (and sloppy) scratchpad implementation I barfed out (is that a decent disclaimer?) is:

transX[ss_] :=
Module[{nc = Length@ss, ssx = Sort@ss, ends, probs, moves, tmp, cur,
fns},
fns[nc_, cs_] :=
Map[{cs, Sort[cs + #]} &, Permutations[UnitVector[nc, nc]]];

ends =
FixedPoint[
Join @@ Map[
If[#[[-1, -1]] != 3, {#},
With[{tmp = #[[-1]] - 2 UnitVector[nc, nc], cur = #},
Join[cur, #] & /@ fns[nc, tmp]]] &, #] &, fns[nc, ssx]];

probs = (1/nc^((Length /@ ends)/2));
moves =
GatherBy[Transpose[{Rule @@@ ends[[All, {1, -1}]], probs}], First];
Transpose[{moves[[All, 1, 1]], Tr /@ moves[[All, All, -1]]}]];


does the job, but terribly inefficiently (I'd like to be able to handle state lists of ~dozens of categories, this dies at 10-ish).

I'm pondering it presently, but ideas invited and welcomed...

• Permutations[UnitVector[nc, nc]] is just Reverse[IdentityMatrix[nc]], no? – J. M.'s ennui Mar 27 '17 at 7:54
• @J.M. - yes, the use of UnitVector is a remnant, but it's certainly not any kind of performance issue there. – ciao Mar 27 '17 at 8:05
• I know, just tossed out something obvious, as I'm not at a computer right now. – J. M.'s ennui Mar 27 '17 at 8:39
• @J.M.- I always figured you are a computer ;-} In any case, I think there's a very simple way to do this that just popped into head, but it's goofy hour of morning already, need to hit the hay... s/b interesting to see what if anything transpires here. I'll ponder my simple idea later and add to OP if it pans out. – ciao Mar 27 '17 at 10:22

Yes, it's possible to compute these explicitly. Let's look at exactly what's happening. Only the individual number of occurrences of $0$, $1$ and $2$ matter, which we'll call $c_0$, $c_1$, $c_2$. Let's also use $n_c = c_0+c_1+c_2$ as the total length of the list.

• Initially, between $0$ and $c_2$ (inclusive) of the $2$s are turned into $1$s (indirectly, by incrementing them to $3$ and then subtracting two). We can never have this reduction case more than $c_2$ times, because to get there we'd need to create additional $2$s which terminates the process.
• Afterwards, either a $0$ is turned into a $1$ or a $1$ is turned into a $2$.

That means there are at most $2\cdot(c_2+1)$ cases to look at. If $c_0$ is zero, we can't turn a $0$ into a $1$, so we'll lose $c_2+1$ cases, and if there is no $1$, we'll need to reduce at least one $2$ before we're able to increment a $1$ and lose one case. Either way, the number of cases we have to consider is very manageable. These cases will also always result in different tallies at the end, because whether we end up incrementing $0$ or $1$ changes the tally of $0$s and then the number of $2$s we decrement changes the split between $1$s and $2$s, so if we look at these cases individually there's no need to recombine results later.

Looking at the problem from the above perspective also gives us easy ways to compute the probabilities:

• If we reduce $n_2$ of the $2$s, we'll have to pick $n_2+1$ of the slots to increment (the $n_2$ twos and the final zero or one). So each of these individual cases has a probability of $1/n_c{}^{n_2+1}$.
• But of course there are many equivalent cases, since we have multiple $0$s, $1$s or $2$s to choose from. In particular, if we end up choosing a zero, there are $c_2!/(c_2-n_2)!$ ways to pick the $2$s and (independently) $c_0$ ways to choose the zero, so the probability of this event is $c_2!/(c_2-n_2)!\cdot c_0/n_c{}^{n_2+1}$.
• If we end up choosing a zero, the number of ways to pick $2$s is unchanged, but we can pick one of $c_1+n_2$ ones, since the $2$s we changed are $1$s as well now. That gives us a probability of $c_2!/(c_2-n_2)!\cdot (c_1+n_2)/n_c{}^{n_2+1}$.

And that's it. Here is a quickly thrown together code which computes these, working with the tallies all the time. I found it easiest if I just filter out the invalid results later, because we can easily detect them due to a negative amount of $0$s or $1$s:

transY[ss_] :=
Module[{nc = Length@ss, ssx = Sort@ss, c0 = Count[ss, 0], c1 = Count[ss, 1], c2 = Count[ss, 2], common},
common = 1/(c2 + 1);
{ssx -> Join @@ ConstantArray @@@ #, #2} & @@@ Select[
Join @@ Table[
common *= (c2 - n2 + 1)/nc;
{
{{{0, c0 - 1}, {1, c1 + 1 + n2}, {2, c2 - n2    }}, c0*common},
{{{0, c0    }, {1, c1 + n2 - 1}, {2, c2 + 1 - n2}}, (c1 + n2)*common}
}
,
{n2, 0, c2}
],
FreeQ[{_, _?Negative}]
]
]


This solves random inputs of length 10,000 in about 0.3 seconds on my machine, but 100,000 uses way too much memory, which isn't too surprising since the output format grows with the square of the input length.

Indeed, if we change the I/O format to simple triplets of $c_0$, $c_1$, $c_2$:

transY[in : {c0_, c1_, c2_}] :=
Module[{nc = c0 + c1 + c2, common = 1/(c2 + 1)},
Select[
Join @@ Table[
common *= (c2 - n2 + 1)/nc;
{
{in -> {c0 - 1, c1 + 1 + n2, c2 - n2}, c0*common},
{in -> {c0, c1 + n2 - 1, c2 + 1 - n2}, (c1 + n2)*common}
},
{n2, 0, c2}
],
FreeQ[{_, _?Negative}]
]
]


Then the same input as above only takes 0.06 seconds and an input where the three inputs sum to 100,000 is easily computed in 4 seconds.

• As a slight improvement, one could just compute the quantity c2!/(c2 - n2)!/nc^(n2 + 1) == FactorialPower[c2, n2]/nc^(n2 + 1) once. – J. M.'s ennui Mar 27 '17 at 16:11
• @J.M. Ah thanks for letting me know about FactorialPower. I first computed it as a Binomial times a factorial but that seemed equally inelegant. – Martin Ender Mar 27 '17 at 16:12
• @J.M. Managed to speed it up a lot more, by calculating the common factor incrementally as I go along. – Martin Ender Mar 27 '17 at 16:27
• "incrementally" - yeah, that's often what one does in code involving factorials/binomial coefficients in traditional languages. I guess I shouldn't be surprised it applies to Mathematica, too. – J. M.'s ennui Mar 27 '17 at 16:30
• Just woke up, saw your response.The latter is exactly what I had in mind when I nodded off in my comment to @J.M. earlier. Made me smile. +1 and accept. – ciao Mar 27 '17 at 20:41