IT Lecture 7 - Channel Capacity for Noisy typewriter.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "My investigations of material presented in [Mackay: Information Theory](http://www.inference.phy.cam.ac.uk/itprnn/book.html) - [Lecture 7](http://videolectures.net/mackay_course_07/)\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 286
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Entropy Functions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "norm = lambda(X): X / sum(X)\n",
      "\n",
      "def H(X): \n",
      "    return sum(i*log2(1.0/i) for i in X if i > 0)\n",
      "\n",
      "def CondHs(C):\n",
      "    return [H(c) for c in C.transpose()]\n",
      "\n",
      "def HofYgivenX(X, C):\n",
      "    return dot(X, CondHs(C))\n",
      "\n",
      "def MutInf(X, C):\n",
      "    return H(dot(X,C)) - HofYgivenX(X,C)\n",
      "\n",
      "def report(X,C, show = False):\n",
      "    np.set_printoptions(precision=3)\n",
      "    if show:\n",
      "        print \"\"\"\n",
      "X:\n",
      "%s\n",
      "C:\n",
      "%s\n",
      "\"\"\" % (X,C)\n",
      "    print \"\"\"\n",
      "H(X): %f\n",
      "Mutual information: %f\n",
      "\"\"\" % (H(X), MutInf(X,C))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 298
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Channels"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Noisy Typewriter channel with failure param $0 \\le f \\le 1$\n",
      "\n",
      " * `noisy(2, f)` is the binary symmetric channel\n",
      " * `noisy(n, 2./3.)` is the noisy typewriter\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def noisy(n, f = .5):\n",
      "    return identity(n) * (1. - f) + array([roll(a,1) + roll(a,-1) for a in identity(n) * (f/2.)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 287
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Trivial Test Case: Reliable Binary Channel"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 2\n",
      "report(norm(ones(2)), identity(2), True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "X:\n",
        "[ 0.5  0.5]\n",
        "C:\n",
        "[[ 1.  0.]\n",
        " [ 0.  1.]]\n",
        "\n",
        "\n",
        "H(X): 1.000000\n",
        "Mutual information: 1.000000\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 299
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Noisy Typewriter"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n=27\n",
      "report(norm(ones(n)),noisy(n, 2./3.))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "H(X): 4.754888\n",
        "Mutual information: 3.169925\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 300
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Channel capacity\n",
      "\n",
      "$C = \\max_{X}I(X;Y)$\n",
      "\n",
      "where \n",
      "\n",
      "$\\begin{align}\n",
      "I(X;Y)& = H(X) - H(X|Y) \\\\\n",
      "      & = H(Y) - H(Y|X)\n",
      "\\end{align}$\n",
      "\n",
      "What is an (the?) optimal input distribution?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Average conditional entropy\n",
      "\n",
      "Any input $x_i$ produces three outputs $x_{i-1}$, $x_{i}$, $x_{i+1}$  with equal probability.\n",
      "\n",
      "$\\forall_{x_i \\in X} H(Y|X=x_i) = \\log_2(3) $\n",
      "\n",
      "Therefore the average $H(Y|X) = \\log_2(3) $.\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "CondHs(noisy(3, 2./3.))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 305,
       "text": [
        "[1.5849625007211559, 1.5849625007211559, 1.5849625007211561]"
       ]
      }
     ],
     "prompt_number": 305
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since $H(Y|X)$ doesn't depend on $X$, we must maximize $H(X)$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "report(norm(zebra(n)),noisy(n, 2./3.))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "H(X): 3.700440\n",
        "Mutual information: 3.033773\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 294
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def zebra(n):\n",
      "    return array([zeros(n/2),ones(n/2)]).transpose().ravel()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 239
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "zebra(27)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 240,
       "text": [
        "array([ 0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,\n",
        "        1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.])"
       ]
      }
     ],
     "prompt_number": 240
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "MutInf(norm(ones(2)), binSym(.1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 217,
       "text": [
        "0.53100440641071867"
       ]
      }
     ],
     "prompt_number": 217
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "AvgCondH([.5,.5], array([[.5,.5],[.5,.5]]))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 160,
       "text": [
        "1.0"
       ]
      }
     ],
     "prompt_number": 160
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n=9\n",
      "C = noisy(n)\n",
      "unif = norm(np.ones(n))\n",
      "stag = norm(array([1,1,0,1,0,0,1,0,0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 186
    }
   ],
   "metadata": {}
  }
 ]
}