Deriving Explicit Formulas from Markov Chains

Once we formulate the markov model correctly, we can obtain the generating function for each entry in the matrix, where there's a possibility of getting the explicit formula. Let's take a look at one such problem:

A six faced unbiased die is rolled :math:`n` times. What is the probability that we get to see all the six numbers in the sequence?

Setting up a markov chain is easy:

\begin{equation*} \displaystyle A= \begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1}{6} & \frac{5}{6} & 0 & 0 & 0 & 0\\ 0 & 0 & \frac{1}{3} & \frac{2}{3} & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{2}{3} & \frac{1}{3} & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{5}{6} & \frac{1}{6}\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \end{equation*}

The states indicate the number of faces shown up. E.g. the row above the last row indicates that when 5 faces are seen, there's a probability of \(5/6\) remaining in the same state and \(1/6\) moving to the final state.

So, \(A^n[0,6]\), gives the required answer. But we can also find the generating function for that entry by computing \((I-x\, A)^{-1}\), which is

\begin{equation*} \displaystyle G(x) = \frac{10 \, x^{6}}{{\left(5 \, x - 6\right)} {\left(2 \, x - 3\right)} {\left(x - 1\right)} {\left(x - 2\right)} {\left(x - 3\right)} {\left(x - 6\right)}} \end{equation*}

and on partial fractions it's

\begin{equation*} \displaystyle G(x) = \frac{36}{5 \, x - 6} - \frac{45}{2 \, x - 3} - \frac{1}{x - 1} + \frac{40}{x - 2} - \frac{45}{x - 3} + \frac{36}{x - 6} + 1 \end{equation*}

and the probability can be written by extracting \([x^n]G(x)\) as

\begin{equation*} \displaystyle \mathbb{P}(n) = 1-\frac{6}{6^n}+\frac{15}{3^n}-\frac{20}{2^n}+15\left(\frac{2}{3}\right)^n-6\left(\frac{5}{6}\right)^n \end{equation*}

which can be verified by a simulation in J:

1
2
3
n=:10
sim=: 3 : '6=+/~:?n#6'
(+/%#)(sim"0)1e5#0 NB. about 0.27

and

\begin{equation*} \mathbb{P}(10) = \dfrac{38045}{139968} \approx 0.271812128486511 \end{equation*}

Multisets and multivariate generating functions

Consider a multiset, \(S = \{11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15\}\). How many combinations of 8 elements can be made from the set so that the sum of those 8 elements is equal to 105, when the numbers are picked without replacement?

Here, if the sum is not asked, the problem can be solved using an ordinary generating function as

\begin{equation*} G(x) = (1+x+x^2+x^3+x^4)^2\, (1+x+x^2+x^3+x^4+x^5)\, (1+x+x^2+x^3+x^4+x^5+x^6)\, (1+x+x^2+x^3) \end{equation*}

and finding \([x^8]G(x)\)

But when the sum is also there as a constraint, we require one more variable to keep track of the sum. So, we may use \(x\) to know the number of elements chosen, and \(y\) for the sum of those numbers. Hence, the required bivariate generating function can be written as

\begin{align*} G(x,y) = {\left(x^{6} y^{84} + x^{5} y^{70} + x^{4} y^{56} + x^{3} y^{42} + x^{2} y^{28} + x y^{14} + 1\right)}\\ {\left(x^{5} y^{65} + x^{4} y^{52} + x^{3} y^{39} + x^{2} y^{26} + x y^{13} + 1\right)} {\left(x^{4} y^{48} + x^{3} y^{36} + x^{2} y^{24} + x y^{12} + 1\right)}\\ {\left(x^{4} y^{44} + x^{3} y^{33} + x^{2} y^{22} + x y^{11} + 1\right)} {\left(x^{3} y^{45} + x^{2} y^{30} + x y^{15} + 1\right)} \end{align*}

Now, the answer to the question would be \([x^8 y^{105}] G(x,y)\) in the expansion of the product.

If the question was to find the number of combinations with replacement, the generating function can be represented as

\begin{equation*} G(x,y) = \dfrac{1}{\left(1-x\, y^{11}\right)\left(1-x\, y^{12}\right)\left(1-x\, y^{13}\right)\left(1-x\, y^{14}\right)\left(1-x\, y^{15}\right)} \end{equation*}

Now, let us focus our attention to a related probability problem:

From the multiset S, what is the probability of choosing 8 elements such that the sum of those 8 elements is equal to 105, when the numbers are picked without replacement?

This looks simple and we may be tempted to say that the answer is \([x^8 y^{105}]G(x,y) / [x^8] G(x)\), but note that some combinations of numbers are more probable to be picked since the number of each element are not the same. E.g. If the set contains \(\{11, 11, 12\}\), the probability of choosing \(\{11, 11\}\) will be more than the probability of choosing \(\{11, 12\}\).

Anyway, it can still be solved using a generating function:

\begin{align*} P(x,y) = \left(1+\binom{4}{1}\, (x\, y^{11})^{1}+\binom{4}{2}\, (x\, y^{11})^{2}+\binom{4}{3}\, (x\, y^{11})^{3}+\binom{4}{4}\, (x\, y^{11})^{4}\right)\\ \left(1+\binom{4}{1}\, (x\, y^{12})^{1}+\binom{4}{2}\, (x\, y^{12})^{2}+\binom{4}{3}\, (x\, y^{12})^{3}+\binom{4}{4}\, (x\, y^{12})^{4}\right)\\ \left(1+\binom{5}{1}\, (x\, y^{13})^{1}+\binom{5}{2}\, (x\, y^{13})^{2}+\binom{5}{3}\, (x\, y^{13})^{3}+\binom{5}{4}\, (x\, y^{13})^{4}+\binom{5}{5}\, (x\, y^{13})^{5}\right)\\ \left(1+\binom{6}{1}\, (x\, y^{14})^{1}+\binom{6}{2}\, (x\, y^{14})^{2}+\binom{6}{3}\, (x\, y^{14})^{3}+\binom{6}{4}\, (x\, y^{14})^{4}+\binom{6}{5}\, (x\, y^{14})^{5}+\binom{6}{6}\, (x\, y^{14})^{6}\right)\\ \left(1+\binom{3}{1}\, (x\, y^{15})^{1}+\binom{3}{2}\, (x\, y^{15})^{2}+\binom{3}{3}\, (x\, y^{15})^{3}\right) \end{align*}

and the required probability is

\begin{equation*} \mathbb{P} = \dfrac{[x^8 y^{105}] P(x,y)}{\dbinom{22}{8}} = \dfrac{343}{2805} \approx 0.122281639928699 \end{equation*}

And what is the probability if we do it with replacement? In this case, the probability can be found by using an exponential generating function, which is written as:

\begin{equation*} E(x,y) = e^{x\, \left(4\, y^{11}+4\, y^{12}+5\, y^{13}+6\, y^{14}+3\, y^{15} \right)} \end{equation*}

and the probability is given by \([x^8 y^{105}]E(x,y)\dfrac{8!}{22^8} = \dfrac{5621995920}{22^8} \approx 0.102449319851133\).

The above probabilities can also be verified by monte-carlo simulations in J, for the without replacement case:

1
2
3
lst=:(4#11,12),(5#13),(6#14),3#15
sim=: 3 : '105=+/(8?#lst){lst'
(+/%#)(sim"0)1e5#0

and for the with replacement case:

1
2
3
lst=:(4#11,12),(5#13),(6#14),3#15
sim=: 3 : '105=+/(?8##lst){lst'
(+/%#)(sim"0)1e5#0

Getting started with experimental mathematics

Here is one nice problem to describe how arrive at a formula experimentally.

To rephrase the problem:

There are \(m\) people with one ball each, and \(n\) boxes. In a round, each of them picks one box randomly (uniformly and independently) and drops the ball in it. Whichever box is not empty is removed, and the next round starts. How many rounds, on an average, will it take till no boxes are left?

Obtaining a formula directly by combinatorial arguments without computing any values and getting it right is quite difficult, and prone to errors.

So, let us obtain it experimentally.

What does experimental math involve?

  • Brute force through the problem
    • Write a program which describes the problem
    • obtain the first few values
    • take it to either oeis or a sequence guessing routine
    • Then we may be able to find a formula
  • If it's a problem on probability, do a simulation to cross verify with the formula
  • Having a lot of fun, doing both math and programming at the same time!

Back to our problem, how many rounds can we expect for the game to last? Let us do the simulation by computing answers for small values, in J:

   'm n'=:5 3
   sim=: 3 : 0
a=:m
c=:0
while. (a>0) do.
b=:+/~:?n#a
a=:a-b
c=:c+1
end.
c return.
)
   (+/%#)(sim"0)1e5#0

Running the above gives a value of about \(2.554\)

Next, we will try to compute some numbers: How many ways is it possible for 3 balls to be placed 5 boxes such that everybody chooses the same box?

\(abc,0,0,0,0\)

\(= 5\) ways

How many ways is it possible for 3 balls to be placed 5 boxes such that two boxes are selected?

Do some casework:

One box may contain two balls, one box with one ball and one empty box.

\(ab,c,0,0,0\)

\(ac,b,0,0,0\)

\(bc,a,0,0,0\)

\(= 3\cdot 5!/3! = 60\) ways

How many ways is it possible for 3 balls to be placed 5 boxes such that 3 boxes are selected?

\(a,b,c,0,0\)

\(= 5!/2! = 60\) ways

And we see that the total turns out to be \(5 + 60 + 60 = 125\), which is 53, the number of ways of arranging the balls in boxes without any restriction.

To calculate the expected value, we have

\begin{equation*} E[n] = p_1 * E[n-1] + p_2 * E[n-2] + \cdots + p_m * E[n-m] \end{equation*}

where \(p_m\) means the probability of choosing \(m\) different boxes from \(n\) boxes.

For \(m=3\) and \(n=5\), \(E[5] = 5/125*E[4]+60/125*E[3]+60/125*E[2] + 1\)

Then, calculate similarly for \(n=4\) and \(m=3\) to get \(E[4]\) and so on. The boundary condition is \(E[1]=1\), since obviously the game would end in one round if there was a single box.

\(E[5]\) would be \(\dfrac{18391}{7200} = 2.5543\) which is close to the simulation. Hence, we can proceed with our experimentation for conjecturing a formula.

Let us calculate the number of ways to partition a number \(n\) of length \(3\) (number of people fixed at \(m=3\)), using sage:

def afun(aa,bb,cc):
    nn = aa
    nbac = nn
    mm = bb
    ll = cc
    summ = 0
    alst = Partitions(nn,length=ll).list()
    for a in alst:
        b = a+([0]*(mm-ll))
        tot = 1
        for c in a:
            tot *= binomial(nn,c)
            nn = nn-c
        summ += Permutations(b).cardinality()*tot
        nn = nbac
    return summ
print [afun(i,3,3) for i in range(3,16)]

What this function afun does is that for \(n\) and \(m\), it computes the number of partitions having length \(l (\le m)\), and we compute the list of values for \(l=m=3\) and varying \(n\).

Insert that list to oeis, and bingo! The second answer shown looks promising: \(A(k,3)\) where \(A(k,n)= \sum_{m=1}^k (-1)^{m+1}\cdot \binom{n}{m} \cdot m^k\).

It seems to be related to the stirling numbers of the second kind.

After some trial and error, the equation turns out to be:

\begin{align*} \displaystyle E_{n,m} &= \left(\sum_{j=1}^{n-1} \left\lbrace {m \atop j} \right\rbrace \dfrac{n!}{(n-j)!} \dfrac{E_{n-j,m}}{n^m}\right)+1\\ E_{1,m} &= 1 \end{align*}

In maxima (which will cache the values to speed up recurrence computation), it can be written as:

1
2
3
4
m:3$
E[1]:1$
E[n]:=sum(stirling2(m,j)*factorial(n)/factorial(n-j)*E[n-j]/n^m,j,1,n-1)+1$
float(E[5]);

Jolla and sailfish OS

Last week, I bought a Jolla phone after looking at some of the amazing features it offers. What are those amazing features? We'll see..

Most of the reviews scattered on the net miss the fact that it is a fully functional linux machine! We can install gcc, python etc. from the terminal (after the developer mode is activated). Let me list some of the things I did after getting my hands on the Jolla phone:

  • I didn't have an access point for wifi, so used hostapd to create one.

  • Create a Jolla account to download the required softwares. (Phone's very much usable even without a sim card)

  • After creating the account, we'll be able to enable the developer mode.

  • Now, the terminal will be displayed.

  • Open the terminal, become a root by typing in devel-su and the password from the settings->system->developer-mode

  • We can also login via ssh, which will make our typing a lot easier. I connected via WLAN's IP address (ssh nemo@ip-addr, and type the password)

  • gcc, g++, gdb, python and strace can be installed by issuing commands like pkcon install gcc

  • Since we are also interested to do some math, what are our options?

    • I downloaded sagemath for arm \(-\) sage-5.13-armv7l-Linux, and tried to run it. Too bad, it crashed saying ImportError: libcblas.so.3: cannot open shared object file: No such file or directory (Update: The sage 6.0 binary available for RPi works for jolla too, woohoo! Only issue I found till now is that plotting did not work.)
    • For FriCAS, Axiom, maxima etc., binaries are not available for ARM. To compile them, they require a lisp compiler, which isn't available in the repos yet (ecl is available in openrepos.net, did not try it though).
    • But anyway, if we look into openrepos, we find sympy, numpy, IPython, and matplotlib, which are more than sufficient for a mobile phone!
  • Why use IPython? Because, it provides us with a notebook interface which can be opened remotely, and also tab completion feature

  • We can also install several useful softwares like emacs, system monitor, cpufrequtils etc. from openrepos.net. Hence, we can have the tools we need with us all the time! (Emacs calc mode itself can serve as a mini CAS!)

  • The history of events like sms and calls are stored in a database. We can use SQL commands to extract the event of interest. E.g.

    sqlite3 .local/share/commhistory/commhistory.db
    
  • Then, to display all the sms messages, run commands like

    select * from events where type=2;
    
  • If we want to know what each column is about, run

    .headers ON
    
  • And to top it all, though an ARM binary is not available, we can easily compile J to run it from the console! I think it's perfectly suitable language that can be used on a phone. (J download for RPi works too)

To conclude, it's totally developer friendly. If anyone likes to program but hates all those manifest files seen in other mobile OSes, it feels like bliss! No need of bloated SDKs, no need of emulators (though it's there if a UI is required), just a simple ssh to connect to the phone and start coding!

It's friendly towards "end users" too, don't believe the "reviews" that say it has a "steep" learning curve because of its gesture based controls. The tutorial it displays when the phone is started is more than enough to understand the Sailfish OS. I won't talk about its UI, since pretty much all the other reviews are about only that and apps, pfff.

Nokia made a big mistake and died by dumping MeeGo and embracing m$.

I wish Jolla team a big success, but I also wish that they remove the mandatory signing up of a Jolla account for enabling developer mode.

1. Creating AP using hostapd

2. Run IPython remotely

3. Sage 6.0 for RPi download