{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Information theory: the noisy binary symmetric channel"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Florent Leclercq,
\n",
"Imperial Centre for Inference and Cosmology, Imperial College London,
\n",
"florent.leclercq@polytechnique.org"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.special import binom\n",
"from matplotlib import pyplot as plt\n",
"from PIL import Image\n",
"np.random.seed(123456)\n",
"%matplotlib inline\n",
"plt.rcParams.update({'lines.linewidth': 2})\n",
"plt.rcParams.update({'font.size': 16})\n",
"plt.rcParams.update({'axes.titlesize': 16})\n",
"plt.rcParams.update({'axes.labelsize': 16})\n",
"plt.rcParams.update({'xtick.labelsize': 16})\n",
"plt.rcParams.update({'ytick.labelsize': 16})\n",
"plt.rcParams.update({'legend.fontsize': 16})\n",
"plt.rcParams.update({'figure.titlesize': 16})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The signal"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Open the image and convert it to an array of black/white pixels\n",
"im_data = Image.open('data/correlation.png')\n",
"im_grey = im_data.convert('L')\n",
"signal = np.array(im_grey, dtype=np.int)\n",
"signal[np.where(signal>1)]=1"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAADqCAYAAABZY8eSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAG4hJREFUeJzt3UmofcldwPFfJS0GJ1BEEQfctPbCARQJughBcKMG3URUUBoHJCpB3KhZaISAgoTWjQgu3AkmQoMaUAzZiJDggGBUlNBoXKQTp0SjtigeF/dd+7x6NVedc+r3q+8HmuR/37lnquFXVadOXbdtmwAAAJ1ec/UJAACAdgRyAAAUI5ADAKAYgRwAAMUI5AAAKEYgBwBAMQI5sBjn3Budc5tz7o0TnMvzD+fypVefC6DVM1efAIDT/amIfL2I/OXVJwKgH4EcWMy2bf8qIu+/+jwAjMHQOmCQc+7LnHMvOuc+5px7xTn3Yefcu51zz4SG1p1zr3XOvcM59xHn3H84597nnHvuYbu377Z7+8Nnzzrn3uOc+6Rz7u+ccz/tnHvNbrvXOedecM598GGbl51zv+2ce+7cOwHYR48csOl3ROTjIvIWEflHEflCEflmiTfef1ZE3iYivyAi7xWRrxGR30rs/0UR+TUReUFE3vTw/b9/+ExE5FNF5DNF5B0i8hER+RwR+WEReb9z7rlt215uvTAAjxHIAWOcc58rIs+KyLdt27YPxr/+8Hd/+88WkR8TkV/Ztu0nHj7+fefcf4vIOyOHeee2bfeg/V7n3DeKyHfJQyDftu0TIvIDu2O8VkR+T0Q++rDdC80XCOARhtYBe/5JRF4SkZ93zv2gc+7ZzPZfKSKfLiLv9j7/zcR33uP9+4Mi8iX7D5xz3+Gc+4Bz7uMi8j8i8u8i8hki8uWZ8wFQgUAOGLPdftLwm0Tkj0Xk50Tkb5xzLznn3hL5yhc8/O/HvM8/mjjMP3v//i8Red39H865N4nIb4jIX4nId4vI60Xk60TkH/bbAejH0Dpg0LZtL4nI97rbOPpXi8iPisgvO+f+VkT+09v8Iw//+3ki8he7zz+/4xS+U0Q+tG3b8/cPnHOfIrdn5QAGokcOGLbd/JmI/PjDR18R2OzP5Tbs/Wbvc//fNT5NbsPpe98jIq/t2CeAAHrkgDHOua8SkV+S29D2h+QWPJ+XW2B9n9xmk/+/bdv+xTn3iyLyNufcv8mrs9a//2GT/204jd8VkW93zr0gtxn0Xysib5XbTHoAAxHIAXteFpEPy60X/kUi8orcet3fum3bn0SWZv0ZEXFyC95vFZEPyC34/6GIfKLhHH5VRL5YRL5PRH5IRP5Ibq+pvdiwLwAJ7jYvBgAec869WUTeJSJv2LbtD64+HwBhBHIA4px7vYh8i9x64q/IbSj8J0Xkr0XkGzYqCmBaDK0DEBH5pIi8QUR+REQ+S26vor1LRH6KIA7MjR45AACK8foZAACKEcgBAFBM5TNy5xzPA3CKmkdP/o+RYE3kGRxl27ZghqFHDgCAYip75HdM1MNRenpK5Ms1kWdwlFzeokcOAIBiBHIAABQjkAMAoBiBHAAAxQjkAAAopnrWek5uph8zReP8e8e9GqNmZvNR9/x+DmekqbYy6Jyb7pzuZj43izTVgeYC+RULLOQqxtA5bdsWPdfQfmqv68pFKfb7mznzz6zkvqXSrfT7qe3ObHSUVpqx7UYGudny7/58Zjs3X6yu06ClHozlu7PTyVwgTym9obHA7H/uJ3xNRkhtO6JSaq0YW46xb5SE7lnpfcxtZ8H+WkKFvSRftOSz2ns4qnGXu97YsVL5Jsb/zoj8c3UvuDQdUgGl5vPQsWN5rjZ9So5Tcg65uuyIvJ5r8KbK5Rl5yFwgH5W5Unr36Sd66Jz9xM/14GsrudB3Q+d19AhHrBLHY/d7UlPp9C5QclTvKhfE/b+l8mWq8qwVa2BeHcxFnt6z0D1M3aPc9v4x9koC1Jk98VDgzOWp1LmU1K2to1dn5RtzgbxXbriytmUaS8iWIJk7Tm4IPxXMj5QqBP55a3oudbWj701pwB19LP94ofKU6wHF9lt7HlcqCQ4tDZnYvStpsITKq/+dGeZe1OrtAJaMGBzJ/Kz1mQpmrd4e1b7Sq/neVQjix5lxlKP2GX9JwBj5bHL2/Odf677MH8k/zqi81bqfUL7w/xthxjJ0R4+8U+lzppJte5w5uaI0Q5eMBqS+u7JYevbc01kdcQ2tz0lnqqyvSOua48Qeyx3RIL966Dp3HlcjkHc6shdZuq8zgnhPa7m2MrIQnEYqeY4skn+Moznwaz3vXiXzAUbfm9Z80jNk36Lkmf4V+eaKkUXzQ+tHKZ0YcuaQV8zVrchVK+Ec7ktcSZ49Kl/PUGb3zjqP1IS6mn2MvH/+aNTeWfVaqoE8C3rkntJJNPdta/5+tZ7zyT2PLFVaQc92785WUmm13KOreyshsVnVtUK915JrnPGeiDy+hpHPpEsnwuWOc3UH4X4OZwbaWUcUlwrk+9cWSgv4fTs/0/a8GhM6jq/2eXJsSKt2f0f2cnL73t/rWSvXmYRew7mLTQzrbXz1NFZD51TyXNXPD7F97L/XEtBT21zRuPQfidTMUA/dh1RdFTpu7N+htEptX3Pf/Do3F6RDeSOV18+aZ3D2cc0G8lymz2WSXAVz5OSYnkQ/qlEwQum+QxWHxWBeE2Ry+wgJ3e9cQ/bsXlbrtaeGXO+f1UzKHLHNEfbBPHQ+ufIeq8vu29QG89h5tHRgYlKNhtx5jjj+fb+pf/fs6whOYwXpnNtE+oZ7a4fcQt8p6Z2kKuvRvc5RQ/uhoc6R+aRm/1f1zFuC7IjAnNp3rld9RB4avX/LjswzqXJQso8R5aimvqs9zog8fWS5yB03VDZH15kP+wxepOlArvHa8NQVaTpbIMf8yDM4Si6Qmxtap0DY0/N8FwCs4/UzqEEjDQCeIpADAKAYgRwAAMUI5AAAKGZushugwVWLVQCwx3wg5/UOHVZ6bTC1WFHLCnwl7w9bv6crGLVOxFlK1ofILT3buiJl6Tv1++01N67NBvLQaj9HLdSxV7vQzFUZ5YrV0vzCWxLQcgvqlLyadvWSjXslad+yKEvu1bzShYlKlZ5HahGlnkVQVhNKo5aGX2n65wJsbF/+d0rOMZX//KVlS8+/tG7InYOWVSWXekZeupTffSlC/7/QNiXHyO2n9NxHvUe9L2Q19yS3v5JjlqrZviSItVzXFVruZeu9as1TsQoyl0d68nrNd1uPocG21f+ymJ/msb+FPm8t9zMEv5Y0Dt1fDfnEZI881JKPZeaalluuZeq3Qmsrn5JeTuk+Y63YVOu1phXt7y/Vcs6NjsR6ba0FqLQX2HK9o9QcL7Zt7n7n0qDkOLG8mUub3NrbLSMLuR56Lj1rerW5Mnmm2LmUjFrF9heqn85YAjhn1JrmI/Jca514BZOB/C4VLGq+K1KWuUMJXzI82NvbLf2efy4lQ1qter5f0wjyv1fz9xkLaq6hMWIIMFaB9wjl89B+exowNeXEvx+1DZfZtKZRrMyn8kBvupWez1HHWJG5QN7Sqz6jEu/df6wiTzUWYj35fUXb2/tNyY2M7LerDd4lDZUSrY2GVvfj1fQYYvvZ76/HyPx/9TPxXGM6tu/e9DjKqLzZe+970qb1sciIxkvo89pyNmvjbm/ZZ+Sh5yAlPeaaAOHvo+WZTc0xSs6t5ju159pzbS0t9Rkq2pFqGqGt+Sn2fLUnb5amXek5146s5I7pKzm+tbwV0poXZg1s/nmPHmmambkeea1Yr/TIobjSHlDsHI7OXC3HSH2npUCNPoeZ1D7f7emZlDSC/Hzf0xtOfZ46zlGNwNx3Zg7YI0ZaUvu6p0EqHWIjaTX55IpHGDOn6xHMBvJUQqaelbbsr3bblkx2VYBqrdhrJhCNLHRHDftdLddIuvdGetKqJyCOOk5vXsg1wEMB6b6N36BfLRiIxOeNjGp0+Ub3oEeO3sXmF8yYL8wG8hqjgnhN4+EILRV5qCLz/311sBtd2I/eb8kxeyqcI861ZQj6TLH71jJqE/uu39A4Y+5MSqqD0dKAin03dp33488UuPbzHlLzg3qPsadhpI9AXqBmKPyIzOQXqJ4h7JL9tDojw4/a/xVB3D9+SzqMfO4XOn7t/kuDzVmPhGLHr+lNzVRx38/lynJ1RJ0RCsal+27pYJTWmbHG3eyWCOSpFmhq+KTk2WGoFZ/6fu0kktDQVotUgUn1Dq8cTqqtUHMt6dJ5D0dqSYeeY6TU5IMSNXk7ts+a897vo+eRSqzc78v2jD3zWN1Tur+aEcTWAD5qlPLIxlVo37HHr7MGdrOBvGfIzf+sZKgtlhlygeRMJRm25Du97vfx6GHkkfMeRipNh5pAVTPPo+b4JT2Z3HYlPf6WxkdLvsxdd+wZ+tWVeGujuqRxlkqjXCOi5dih7WrLQ6nUSEJtQ3DWIC4i4mY+uRjn3CZS16ovSdBUhVM7NJN6Fhf7Xug7GtKnppI7okI8ImC33H9NaXaU3nswogzGziUVLGJ/q3nsQJ7BUXb5JFgITAbyh22yfw/tI1Ww/CF5jfcOZaiU1xMr06VlnTyDo+QCudmh9ZKhyNq/9TyTAzC33jkCwFWWWtmtBIUWAKAJgRwAAMXMDq0DaFM7qexqPGfG6gjkAIpf/aqZ9OXrmZxaggmoWBWBHDhZrgd5RA8zFeRqFhIpOU7ub7HXvnpnh+de+QxtB1iwxDPyqxdiGcG5Y35mEtepXbWuJw+k7BeK8f///r9RUiuR1S744m9fco8oD7DGdI88VCn4rfZRK4fF3j8t/fteaiGM1PvvqWMw7DinknSJrTZ217rql7+P/b5H5Jdcfoxp6aXnzqNmYRftmDOwHpM98ppWea6SjH3e+u+SXkNtr6KkcoROtb3S2bUsvZkbJg+NJOS2AywxF8hDvZXSiiIX3HsDZsv3/J5E7jlgbGjyrArsqOFfy2qeUVsLSLkA3LKv3GeANeYC+Z6/ElvvsGPsb6kAG5qQE/pOroJmVTl7atJx1COgo9X0rluWMgXwlKlAvn82VLrOeuwHFWqOmeqF+s/4aic43c8zt81+O3rFeoSeTd/VpmEq/WsbAiN+9KTmb7G/t5Qd8j5WYyqQlwjNcq1VO8RZ+tx9BHrrerXmh9TrXL3n0lNGWrW+VobHZrpX+w7PkefVc4zS7+23a5n7dMT1m5m1PkOmneEcoI8/qzr22d5+tCY1d6PnVaye2eP77Wsm7I1siNKoPUYuL5Smt58/RqdXajQ092ZQyZsgqTlVubw/4m2MPTM98hHPv3sXzAjt+4iJSbkhyNR10NjQozQ/pnqwR59LSmyI338EUFKR914jeX+M0fcwFehGHis2gTl2DrG/1zzeyb1tMZKZHnmJ3ESbXA8k1IL09zkq8/k9rVALr+Z4Pb0rHK+0B+5/p3V/JfsYKXeckkB9v6bafRHAx6upQ66sa2pGpFL1e6z8hRqhvfNTWpjpke+FWl+lrSM/QVK9hZoWV29i+tdQU5mFZshjPvtXsfzPYtvu/137JsQRSo9fO88ktP/a76yS72cLnCXPhkedc20jLyXVWCkZLTszHUz1yHMtqP12tfuNHcffpjbApnr+sV5Va0t4lYpMs1jQ02TG853xnCyL1V0lPfkRkzRD51Dz6CalpFN3dn4zFchF0kE2973QPmpbZKXHKN3PfrjHQiUP0gx6jXx8eIbSxzG5jl/tNacmvB1R/s0FcpExN2q2yna28wGwtiN6u0eIPb8e2SAJBe2e0dRaJgM5AOBcGnvq+3/v1b4O6u/n7GflJie7AQCu0fvWwNk9+NTxzmicjHjVjh45htHUIgdWcNarpj0LD6X+nuoZXz1k36p2saQSBHIMp7WAAcjzg6nGBnzunGvrsNo1IEY3sAjkAIBqo16Tjc3qPrLnPfKto1gQL2nkjLo2AjkAGHN0Tzm2tsaoAJn6vGVfPfssPWao4XHWowACOYbQ/twK7Uh7iJyT/jMP449skNRi1jqAZjNXrLixmEY0Gh8jkAPoRsWKs5DXniKQoxtDq8B8LJZHi9c0AoEcAADFCOToYvH5G+qRD+ZG+thGIMcQDHmtiXTHSkYsp3oEAjkAGEVDaw28R7643pWHqCgAHc5ad92qGXvidwTyBdVkSGakA8CrZqwLGVpfSM/znZlbo5gDeWROMwYejEWPfAGpCrakkN+/z9AcAMyHQF4p9Us3M+r5VaL99gRzxGj9KcuVUGb7zJ6/CeQFShJxv80shSb284Aj9jvLNWIu5A3gfDwjz2hpic3QehsdxKmckUL+wApmzef0yBNqny3vt5+phz7q+AyhosSoXnlpXru6fAFXI5AH5AJ47ezvs4cbz2hEMIQKX29Dr/eNCvIjjqCh88LQeiUNiQpcrbaclCxMlAvUsy6fCRtmbijSI/ccVRFc0YOdOePBppFvOIS+m3ukNeK4wJ2WhiE98p2WRJutwtCS8WDXvky0vPFR0vv2j+dvTznASLPV8z565I3uCUuFgSvNOiHsirUHmIyJkfMlNM29IJA/KP3xkH3llBP6joZMUcLKdWhTE6j2+e+oNQVKjp8yemImixetIVcOehbumumNo1IE8kK1Lf2rM8DVx8f1UnnWcpCjZ16ndwnns9WsVhmaP5HavmSfMyKQH0BTBqhB5Xgty/ffapmZVWleumIkJ6X2fGKNulhDVtNw+h6BXMZUkNZaeJhLydoGLfu8Mk+e0TC5+hpn1Puu/1X3s3XIu2bypda8wqz1TrUzbLXS+NxoFZZ76hhrRF7RnN9q36jQgkDeaPYAbimTrm5kWs6cZ3Gskfno7PqFjkQaQ+sJZ64XrQWFSK8V027Fa27BBEHd6JEDyrQEp5UCGgHpqSPmWFxhpXxcg0B+EjIgRhhR6ZIXIfLq48GWtQnOpKWRcSUCeQTD6q/S+kqGRS3rGZBu8LWsTGklH1m5jj0CuTGjZ2USxHVKBXArDcwQJkU9FXr3OvVLcbkfq7nqvpKecUx2O8gsFQrv0dpEmqJV7fKm1vKaxYYsPfKIUYl9RSEY8UtQ/LbzvKxUrEeNHmEcK3nNOgJ5QmvFMEOF0loAr/qBDWAk8mw/7qEeywfy1l/RKdn+6oLg93hyP6JR82ME0G+GtB3VK5+h8WxFboLkmfeadC2z/DPyktcvchO+Zg6A/vXN+vvVOMfsFWPLnA5GkPpwv/RbPpCLlL+bqzUItrx7zCS5eVlLm1Bjs/Z3o/f7Qp6W+3TkYjVa7kEJAvkD6wttxCbAWf0RAaTNlldTI0eleXS2a5qR5ns0ogFrtY4jkA+gsXCU9ng0XtsKRgxBzya2SMns5w1cjUDuSb26RVDDTGqCuaZ8XDI6pml9cEuuWCBqn9Y9nQtNZaAWgXwn9/615R4qlaJOVvOkxWvCdazXb8u/foY06wVgdiMmYVrsiYSumbyap/UVv565PCu81UCPHDCgpHKzWIHhfDPko5Ih/pXeaiCQR8QqRqtDmSkrXrM1ltKPnnef2vI8y/0OzYtY8dfbQgjkiLZcZynAK/OHFDX8fnQrfmnvPFp7tC0/v7oCAnkDeqi4wj3PpZ4XWsiXPQGdsllHa6Mwt8rmagGfQB6wSuLn0Cuf24pps+I1X22mAJ6j6VxHYtb64qgYbdGcnqtWwjO6/3CK5jRZaWEhAnkFzZk6x/K1WWS1UmI51mPU3C/urT4E8oRZfsrvKLyyBE1KA7uFsjlSaG5FaBvtPfC91fIAgdyzWgYoxX3BWQjKx7ESqFtZzU8E8kqrFwRcL/dakNXK6o4yCDxGII9o/T1kLSwu2wk7/IYJz8vHSC10ZYWlaylFIN9ZMQOE+BUjFSVmkXvOC+yF8oTFep5A3kB7hdHaG7dYADCvnnJGXn2Ke2IXgTxgxeBm5TqsWy2d9mVxpXkAaJPqpGjvgKUQyB/UVgxaK5XeZ+OarnU1WvNkTs1rUVbvQa9Yubd6v3JLtFq6VhEC+RMrteLurGVqq5igiCNYyUsr/+gOgVzKA1nqByo0BEMCAbC2XLnXUI+1KFkURzMCOWCItsYl5qA9wK3cGxchkD/Skgm0VJw1vfHVFxyBbqtW5jG15VVb+W49X23XmUIg3wklrPVf0LF2PcAe+ds20veGQN4g9Kxcyw8O9Jyj9edMM6sZUbE+emLxms5gPd+U/siVxfqLQL5jMYHvRjU0LN8jwJpcQNYYsO9GBGfN179HIJeyIOdvYyUD7Fm/vlVZT0vmdLRLvW47832z3sOuRSCvRKbB2VoqLav5dObgMjNL+aG0PKyUVwjkHVbKKJiL9UmYI1gKXi1K8sSqoxnWrpFA3sBiBWHxmlZneXJizciExeuHzXzdikCeYKGllrJqa1yr3vf/tUutrIg4/z6VlOuZ6wEm7j5FIO80WyaHPfQq8yiHj9Xcj9i2ml6r9ZWsCSJip2wRyAuEEtpC4gMaUfaAxwjkA2jvDVAx6kXaISWVP1bJOytcJ4F8UanGx8zPx6D/V/haWL2uWXB/dSOQd1ihpQdAHyvPfluk3jO32mAhkEdY/8UgQKPUSmTIo56yiUC+IJY3tGGfdisMr+eset0tLK8xsCICeSerBYHAMKcV02LFa8ZYVuvpOwJ5Rk0G0FbhWM/cAMppq79CSkYbLVynj0C+GIuZGABWRiAPqA12Vnu2Vq/LgtwiRVYabMznGMdKnsBTBPLBZi4sVIq6zZy3oMeqM/8tXyeBHIBqlitotGn5gSHNCOQJNRWE1crkfl3WMj50KS1f5FOsiEDuGVERUJlgtNLHIpaek2s//5mU3kvuuU4E8kXwfBxAzCp1gtXrJJAPZDWTMLyO2VkaiTiS1Tpqb8XHMATyiN4MbymTYB4rVMS+Fa/5SlrrLq3nPQKBfGeFjEClqE9PvtSap7We94y4l/YRyAebcYhv5HnMck2IW7WxNmPZ02SVfGPxOgnkALAAy78bcWcxSJcgkB/AYs/A4jVps2IlteI1o15rnWSlLiOQAxOzUtHUWPGaj8LvRqyBQB4wIjPPWCB6z4le+XVq027VtFr1unNayr6W+9dyntbyCYH8QNu2TRHQZzgH1LNQwdRi4aJxVsw/qyKQ78wSeI8w6rru98jqfYJ+1npbZ1ulbFu6TgI5YJSligrtevKBpobQKtcZQiAHAGO0BybUeebqEwAQZnXSZYmREzNX1noftm1T0RjoPcf9dTrn1OYbAjmAaWitSGczIhCTFnowtA4ABhGIy1i4TwRyAABE79wCAjkAQDULveoeBHIAgEqj1rTQ3hBgshsAYHmagzk9cgAAFCOQAwCgGIEcAADFCOQAACimerKb1nf+YBv5ErXIM+hBjxwAAMWc5in3AACsjh45AACKEcgBAFCMQA4AgGIEcgAAFCOQAwCgGIEcAADFCOQAAChGIAcAQDECOQAAihHIAQBQjEAOAIBiBHIAABQjkAMAoBiBHAAAxQjkAAAoRiAHAEAxAjkAAIoRyAEAUIxADgCAYgRyAAAUI5ADAKAYgRwAAMUI5AAAKEYgBwBAMQI5AACKEcgBAFCMQA4AgGIEcgAAFPs/mJSMKdy3RnMAAAAASUVORK5CYII=\n",
"text/plain": [
"