Determine the limit distribution for a simple Markov chain based on the transition matrix

TransitionMatrixMarkovChain.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following matrix is the transition matrix of a random walk on a circle with N=4 points and transition probability p = 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.2,  0.4,  0. ,  0.4],\n",
       "       [ 0.4,  0.2,  0.4,  0. ],\n",
       "       [ 0. ,  0.4,  0.2,  0.4],\n",
       "       [ 0.4,  0. ,  0.4,  0.2]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K = np.asarray([[0.2, 0.4, 0, 0.4], [0.4, 0.2, 0.4, 0], [0, 0.4, 0.2, 0.4], [0.4, 0, 0.4, 0.2]])\n",
    "K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a first approach to calculating the limit distribution, we take a high power of that matrix - the rows should then be close to the limit distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.25000914,  0.24999086,  0.25000914,  0.24999086],\n",
       "       [ 0.24999086,  0.25000914,  0.24999086,  0.25000914],\n",
       "       [ 0.25000914,  0.24999086,  0.25000914,  0.24999086],\n",
       "       [ 0.24999086,  0.25000914,  0.24999086,  0.25000914]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# \n",
    "# Be careful: K**20 would not be the right thing as this would take the power element-wise\n",
    "# We could either use np.matmul in a loop or np.linalg.matrix_power\n",
    "#\n",
    "P = np.linalg.matrix_power(K,20)\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No let us try to find the eigenvalues and eigenvectors of the transposed matrix (which, as K is symmetric, is in this case simply K itself)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "w,v = np.linalg.eig(K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the [scipy documentation](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig), the first return value w is an array with the eigenvalues and the second value is an array such that column i is an eigenvector for the eigenvalue w[i]. So we expect that one of the entries in w is 1.0. We need to pick that column in v, and renormalize it so that the sum of its components is one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.6,  1. ,  0.2,  0.2])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.25,  0.25,  0.25,  0.25])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigenvector = v[:,1] / np.sum(v[:,1])\n",
    "eigenvector"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}