{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Maximum entropy principle" ] }, { "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": [ "from scipy.optimize import fsolve\n", "from math import log, pow\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dice example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The constraints are $\\sum_{n=1}^{6} n p_n = 4$ (contraint on the mean) and $\\sum_{n=1}^{6} p_n = 1$ (normalization constraint). The Lagrangian is\n", "\\begin{equation}\n", "\\mathcal{L}[\\{p_n\\},\\lambda,\\mu] = - \\sum_{n=1}^{6} p_n \\ln p_n - \\lambda \\left( \\sum_{n=1}^{6} n p_n - 4 \\right) - \\mu \\left( \\sum_{n=1}^{6} p_n - 1 \\right).\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "mean=4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\dfrac{\\partial \\mathcal{L}}{\\partial p_n} =0$ fixes $p_n = \\dfrac{\\mathrm{e}^{-\\lambda n}}{Z} = \\dfrac{x^{-n}}{Z}$ (we define $x \\equiv \\mathrm{e}^\\lambda$)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def p1(x):\n", " return pow(x,-1)/Z(x)\n", "def p2(x):\n", " return pow(x,-2)/Z(x)\n", "def p3(x):\n", " return pow(x,-3)/Z(x)\n", "def p4(x):\n", " return pow(x,-4)/Z(x)\n", "def p5(x):\n", " return pow(x,-5)/Z(x)\n", "def p6(x):\n", " return pow(x,-6)/Z(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The normalization constraint $\\sum_{n=1}^{6} p_n = 1$ fixes\n", "\\begin{equation}\n", "Z = \\sum_{n=1}^{6} \\mathrm{e}^{-\\lambda n} = \\dfrac{1-\\mathrm{e}^{-6\\lambda}}{\\mathrm{e}^\\lambda-1} = \\dfrac{1-x^{-6}}{x-1}.\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def Z(x):\n", " return (1-pow(x,-6))/(x-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The constraint on the mean is obtained by noting that\n", "\\begin{equation}\n", "- \\frac{\\mathrm{d} \\ln Z}{\\mathrm{d} \\lambda} = -\\frac{1}{Z} \\frac{\\mathrm{d}Z}{\\mathrm{d}\\lambda} = \\sum_{n=1}^6 n \\frac{\\mathrm{e}^{-\\lambda n}}{Z} = \\sum_{n=1}^6 n \\, p_n = 4.\n", "\\end{equation}\n", "This gives an equation for $\\lambda$: $\\mathrm{e}^\\lambda/(\\mathrm{e}^\\lambda-1) - 6/(\\mathrm{e}^{6\\lambda}-1)= 4$. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.8397685748659793" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# we now solve the equation giving x=exp(lambda)\n", "def f(x):\n", " return x/(x-1) - 6/(pow(x,6)-1) - mean\n", " # This is -dlnZ/dZ-mean, which should be 0\n", "x0=fsolve(f,2)[0]\n", "x0" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.10306524522362508,\n", " 0.12273053351641969,\n", " 0.14614804267474113,\n", " 0.17403371244043664,\n", " 0.20724008691110057,\n", " 0.246782379233677)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# using the solution for x, we have the probability assignment:\n", "pn = (p1(x0), p2(x0), p3(x0), p4(x0), p5(x0), p6(x0))\n", "pn" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGu5JREFUeJzt3X+cVXW97/HXW8DAH/iLsUzEwROpFALjQAk5/iLBh4pZdkEr8ZpxtGPpw4ea2ZWM8h4zq6Ne88dVTtrhHDhp1NQhFQNBNHVAEK8/SETUEa8iFKABOvC5f+wFdzsMs9Ywe82e2byfj8d+zF5rf9fanzU+nDff71rruxQRmJmZtWa3chdgZmadn8PCzMxSOSzMzCyVw8LMzFI5LMzMLJXDwszMUjkszMwslcPCzMxSOSzMzCxV93IXUCp9+vSJ6urqcpdhZtalLFy48J2IqEprVzFhUV1dzYIFC8pdhplZlyLp1SztPAxlZmapHBZmZpYq17CQNEbSUknLJF3VwueXSXpe0hJJf5J0aNFnmyUtTl71edZpZmaty+2chaRuwK3A54FGoEFSfUQ8X9RsEVAbEX+XdBFwAzAu+WxDRAxpTw0ffPABjY2NbNy4sT27sTbo2bMnffv2pUePHuUuxcxKKM8T3MOBZRGxHEDSNOAMYFtYRMScovZPAF8tZQGNjY3svffeVFdXI6mUu7YWRASrV6+msbGR/v37l7scMyuhPIehDgZeL1puTNbtyNeBPxYt95S0QNITkr6wMwVs3LiRAw44wEHRQSRxwAEHuCdnVoHy7Fm09Be6xcfySfoqUAscV7S6X0SslHQYMFvSsxHxcrPtJgITAfr169dyEQ6KDuXft1llyrNn0QgcUrTcF1jZvJGkUcD3gLERsWnr+ohYmfxcDjwCDG2+bUTcGRG1EVFbVZV6T4mZme2kPHsWDcAASf2BN4DxwDnFDSQNBe4AxkTE20Xr9wP+HhGbJPUBRlI4+d0u1Vf9V3t38SErrj81tc3NN9/MbbfdRk1NDVOnTm2xzYIFC7j33nu5+eab2/T95513HqeddhpnnXUWF1xwAZdddhkDBw5s0z7MzLLILSwioknSxcCDQDdgSkQ8J2kysCAi6oGfAHsBv06GL16LiLHAkcAdkrZQ6P1c3+wqqi7jF7/4BX/84x9bPeFbW1tLbW3tduubmpro3j3bf6K77rprp2s0sx2447j0Np3BP87N/Styvc8iImZGxCcj4h8i4rpk3aQkKIiIURHx0YgYkrzGJusfj4hBETE4+Xl3nnXm5cILL2T58uWMHTuWn//85zz11FOMGDGCoUOHMmLECJYuXQrAI488wmmnnQbAtddey8SJEzn55JM599xzP7S/iODiiy9m4MCBnHrqqbz99rbOGMcff/y26U4eeOABampqGDx4MCeddBIA7733Hueffz7Dhg1j6NCh/O53v+uIX4GZVYiKmRuqM7r99tt54IEHmDNnDn369GHdunXMmzeP7t278/DDD3P11Vdz//33b7fdwoULmT9/Pr169frQ+hkzZrB06VKeffZZ3nrrLQYOHMj555//oTarVq3iG9/4BvPmzaN///6sWbMGgOuuu44TTzyRKVOm8Le//Y3hw4czatQo9txzz/x+AWZWMRwWHWjt2rVMmDCBl156CUl88MEHLbYbO3bsdkEBMG/ePM4++2y6devGxz/+cU488cTt2jzxxBPU1dVtG/baf//9AXjooYeor6/nxhtvBAqXFb/22msceeSRpTo8M6tgDosOdM0113DCCScwY8YMVqxYwfHHH99iu9b+tZ92aWpEtNgmIrj//vs5/PDD21SzmRl4IsEOtXbtWg4+uHBf4i9/+cs2b19XV8e0adPYvHkzb775JnPmzNmuzTHHHMPcuXN55ZVXALYNQ40ePZpbbrmFiMKtLosWLdrJozCzXdEu1bPIcqlrnq688komTJjAz372sxaHkNKceeaZzJ49m0GDBvHJT36S447b/kqNqqoq7rzzTr74xS+yZcsWDjzwQGbNmsU111zDpZdeylFHHUVEUF1dzR/+8IdSHJaZ7QK09V+aXV1tbW00f/jRCy+84DH5MvDv3SrGLnDprKSFEbH9tfvNeBjKzMxSOSzMzCyVw8LMzFI5LMzMLJXDwszMUjkszMws1S51n0XJL4PLcLnaXnvtxbvvvtvur7r22mvZa6+9uPzyyzNvk+W7i/c7adIk6urqGDVqVHvLNbMKs2uFhbVq8uTJ5S7BzDopD0N1kIjgiiuu4NOf/jSDBg1i+vTpALz77rucdNJJ1NTUMGjQoA9NHX7ddddx+OGHM2rUqG3TmQO8/PLLjBkzhqOPPppjjz2WF198EYBXXnmFY445hmHDhnHNNdfssJYd7fe8887jvvvuA6ChoYERI0YwePBghg8fzvr169m8eTNXXHEFw4YN46ijjuKOO+4o6e/IzDov9yw6yG9+8xsWL17MM888wzvvvMOwYcOoq6ujqqqKGTNm0Lt3b9555x0++9nPMnbsWJ5++mmmTZvGokWLaGpqoqamhqOPPhqAiRMncvvttzNgwACefPJJvvnNbzJ79mwuueQSLrroIs4991xuvfXWFutYuHDhDve71fvvv8+4ceOYPn06w4YNY926dfTq1Yu7776bffbZh4aGBjZt2sTIkSM5+eSTW32wk5lVBodFB5k/f/626cU/+tGPctxxx9HQ0MApp5zC1Vdfzbx589htt9144403eOutt3j00Uc588wz2WOPPYDCtOVQ6Ik8/vjjfPnLX962702bCo8uf+yxx7Y9H+NrX/sa3/nOd7arY0f7LbZ06VIOOugghg0bBkDv3r2BwjTnS5Ys2db7WLt2LS+99JLDwmwX4LDoIDuag2vq1KmsWrWKhQsX0qNHD6qrq9m4cSPQ8nTkW7ZsYd9992Xx4sUt7i9tCvMsbVqb5vyWW25h9OjRqd9hZpXF5yw6SF1dHdOnT2fz5s2sWrWKefPmMXz4cNauXcuBBx5Ijx49mDNnDq+++uq29jNmzGDDhg2sX7+e3//+90DhX/n9+/fn17/+NVD4A/7MM88AMHLkSKZNmwYUQmhHdbS032JHHHEEK1eupKGhAYD169fT1NTE6NGjue2227Y9tOkvf/kL7733Xgl/S2bWWe1aPYsOeKj5jpx55pn8+c9/ZvDgwUjihhtu4GMf+xhf+cpXOP3006mtrWXIkCEcccQRANTU1DBu3DiGDBnCoYceyrHHHrttX1OnTuWiiy7iRz/6ER988AHjx49n8ODB3HTTTZxzzjncdNNNfOlLX2qxjtb2u9Xuu+/O9OnT+da3vsWGDRvo1asXDz/8MBdccAErVqygpqaGiKCqqorf/va3+fzCzKxT8RTlVnL+vVvF8BTl23gYyszMUjkszMwsVcWHRaUMs3UV/n2bVaaKDouePXuyevVq/wHrIBHB6tWr6dmzZ7lLMbMSq+irofr27UtjYyOrVq0qdym7jJ49e9K3b99yl2FmJVbRYdGjRw/fXWxmVgIVPQxlZmal4bAwM7NUDgszM0vlsDAzs1QVfYLbzDrYLjA9xq7KPQszM0uVa1hIGiNpqaRlkq5q4fPLJD0vaYmkP0k6tOizCZJeSl4T8qzTzMxal1tYSOoG3AqcAgwEzpY0sFmzRUBtRBwF3AfckGy7P/B94DPAcOD7kvbLq1YzM2tdnj2L4cCyiFgeEe8D04AzihtExJyI+Huy+ASw9dbf0cCsiFgTEX8FZgFjcqzVzMxakWdYHAy8XrTcmKzbka8Df9zJbc3MLEd5Xg3V0oOeW5zRT9JXgVpg66UUmbaVNBGYCNCvX7+dq9LMzFLl2bNoBA4pWu4LrGzeSNIo4HvA2IjY1JZtI+LOiKiNiNqqqqqSFW5mZh+WZ1g0AAMk9Ze0OzAeqC9uIGkocAeFoHi76KMHgZMl7Zec2D45WWdmZmWQ2zBURDRJupjCH/luwJSIeE7SZGBBRNQDPwH2An4tCeC1iBgbEWsk/ZBC4ABMjog1edVqZmaty/UO7oiYCcxstm5S0ftRrWw7BZiSX3VmZpaV7+A2M7NUDgszM0vlsDAzs1QOCzMzS+WwMDOzVA4LMzNL5bAwM7NUDgszM0vlsDAzs1QOCzMzS+WwMDOzVA4LMzNL5bAwM7NUDgszM0vlsDAzs1QOCzMzS+WwMDOzVLk+Kc/MWnHHceWuIJt/nFvuCqwTcM/CzMxSOSzMzCyVw8LMzFI5LMzMLJXDwszMUvlqqET1Vf9V7hIyWXH9qZnaVdrxQOUd05I31uZcSWkc1Ya2lXZMlXY87eGehZmZpXJYmJlZKoeFmZmlcliYmVkqh4WZmaXKFBaSuuVdiJmZdV5ZexbLJP1E0sBcqzEzs04pa1gcBfwFuEvSE5ImSuqdY11mZtaJZAqLiFgfEf87IkYAVwLfB96UdI+kT+RaoZmZlV3mcxaSxkqaAdwE/BQ4DPg9MLOV7cZIWippmaSrWvi8TtLTkpokndXss82SFiev+jYdlZmZlVTW6T5eAuYAP4mIx4vW3yeprqUNkpPitwKfBxqBBkn1EfF8UbPXgPOAy1vYxYaIGJKxPjMzy1HWsDg3IuYXr5A0MiIei4hv72Cb4cCyiFietJ8GnAFsC4uIWJF8tqWthZuZWcfJeoL75hbW3ZKyzcHA60XLjcm6rHpKWpCcUP9CG7YzM7MSa7VnIekYYARQJemyoo96A2n3XqiFddGG2vpFxEpJhwGzJT0bES83q28iMBGgX79+bdi1mZm1RVrPYndgLwqhsnfRax1wVivbQaEncUjRcl9gZdbCImJl8nM58AgwtIU2d0ZEbUTUVlVVZd21mZm1Uas9i4iYC8yV9MuIeLWN+24ABkjqD7wBjAfOybKhpP2Av0fEJkl9gJHADW38fjMzK5G0Yah/iYhLgf8labshpIgYu6NtI6JJ0sXAgxSGrKZExHOSJgMLIqJe0jBgBrAfcLqkH0TEp4AjgTuSE9+7Adc3u4rKdkH1u3+v3CVklP2BTmZdRdrVUL9Kft64MzuPiJk0uw8jIiYVvW+gMDzVfLvHgUE7851mZlZ6acNQC5OfczumHDMz64zShqGepZUrmCKiIx79amZmZZY2DHVah1RhZmadWtowVFuvgDIzswrU6n0WkuYnP9dLWtf8Z8eUaGZm5ZbWs/hc8nPvjinHzMw6o6wTCSKpBvgchRPe8yNiUW5VmZlZp5L1eRaTgHuAA4A+wC8l/Y88CzMzs84ja8/ibGBoRGwEkHQ98DTwo7wKMzOzziPrFOUrgJ5Fyx8BXm65qZmZVZq0m/JuoXCOYhPwnKRZyfLngfmtbWtmZpUjbRhqQfJzIYUJ/7Z6JJdqzMysU0q7dPaejirEzMw6r0wnuCUNAP4ZGEjRuYuIOCynuszMrBPJeoL7X4HbgCbgBOBe/v/05WZmVuGyhkWviPgToIh4NSKuBU7MrywzM+tMst5nsVHSbsBLydPv3gAOzK8sMzPrTLL2LC4F9gC+DRwNfA2YkFdRZmbWuWTqWSSPPyXpXXw7ItbnWpW1m59XbWallHVuqNrkqXlLgGclPSPp6HxLMzOzziLrOYspwDcj4lEASZ+jcIWUH6tqZrYLyHrOYv3WoACIiPmAh6LMzHYRaXND1SRvn5J0B/AfFOaGGoen/DAz22WkDUP9tNny94veR4lrMTOzTiptbqgTOqoQMzPrvLJeDbWPpJ9JWpC8fippn7yLMzOzziHrCe4pFE5o/7fktY7C1VBmZrYLyHrp7D9ExJeKln8gaXEeBZmZWeeTtWexIbm3AgBJI4EN+ZRkZmadTdaexYXAvUXnKf6K54YyM9tlpIZFMh/U4RExWFJvgIhYl3tlZmbWaaQOQ0XEFuDi5P06B4WZ2a4n6zmLWZIul3SIpP23vnKtzMzMOo2sYXE+8E1gLrCg6NUqSWMkLZW0TNJVLXxeJ+lpSU2Szmr22QRJLyUvnx8xMyujrCe4B1IIi89RmObjUeD21jaQ1A24Ffg80Ag0SKqPiOeLmr0GnAdc3mzb/SlMLVKbfN/CZNu/Zqy3zfz8BzOzHcvas7gHOBK4GbgleX9PyjbDgWURsTwi3gemAWcUN4iIFRGxBNjSbNvRwKyIWJMExCxgTMZazcysxLL2LA6PiMFFy3MkPZOyzcHA60XLjcBnMn5fS9senHFbMzMrsaw9i0WSPrt1QdJngMdStlEL67LOVJtpW0kTt85XtWrVqoy7NjOztsoaFp8BHpe0QtIK4M/AcZKelbRkB9s0AocULfcFVmb8vkzbRsSdEVEbEbVVVVUZd21mZm2VdRhqZ84XNAADJPUH3gDGA+dk3PZB4H9K2i9ZPhn47k7UYGZmJZApLCLi1bbuOCKaJF1M4Q9/N2BKRDwnaTKwICLqJQ0DZgD7AadL+kFEfCoi1kj6IYXAAZgcEWvaWoOZmZVG1p7FTomImcDMZusmFb1voDDE1NK2UyhMjW5mZmWW9ZyFmZntwhwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapcg0LSWMkLZW0TNJVLXz+EUnTk8+flFSdrK+WtEHS4uR1e551mplZ67rntWNJ3YBbgc8DjUCDpPqIeL6o2deBv0bEJySNB34MjEs+ezkihuRVn5mZZZdnz2I4sCwilkfE+8A04Ixmbc4A7kne3wecJEk51mRmZjshz7A4GHi9aLkxWddim4hoAtYCBySf9Ze0SNJcScfmWKeZmaXIbRgKaKmHEBnbvAn0i4jVko4GfivpUxGx7kMbSxOBiQD9+vUrQclmZtaSPHsWjcAhRct9gZU7aiOpO7APsCYiNkXEaoCIWAi8DHyy+RdExJ0RURsRtVVVVTkcgpmZQb5h0QAMkNRf0u7AeKC+WZt6YELy/ixgdkSEpKrkBDmSDgMGAMtzrNXMzFqR2zBURDRJuhh4EOgGTImI5yRNBhZERD1wN/ArScuANRQCBaAOmCypCdgMXBgRa/Kq1czMWpfnOQsiYiYws9m6SUXvNwJfbmG7+4H786zNzMyy8x3cZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpbKYWFmZqkcFmZmlsphYWZmqRwWZmaWymFhZmapHBZmZpYq17CQNEbSUknLJF3VwucfkTQ9+fxJSdVFn303Wb9U0ug86zQzs9blFhaSugG3AqcAA4GzJQ1s1uzrwF8j4hPAz4EfJ9sOBMYDnwLGAL9I9mdmZmWQZ89iOLAsIpZHxPvANOCMZm3OAO5J3t8HnCRJyfppEbEpIl4BliX7MzOzMsgzLA4GXi9abkzWtdgmIpqAtcABGbc1M7MO0j3HfauFdZGxTZZtkTQRmJgsvitpaZsqzF8f4J2S7vGHLf1qOkylHQ9U3jFV2vFA5R1TZzueQ7M0yjMsGoFDipb7Ait30KZRUndgH2BNxm2JiDuBO0tYc0lJWhARteWuo1Qq7Xig8o6p0o4HKu+Yuurx5DkM1QAMkNRf0u4UTljXN2tTD0xI3p8FzI6ISNaPT66W6g8MAJ7KsVYzM2tFbj2LiGiSdDHwINANmBIRz0maDCyIiHrgbuBXkpZR6FGMT7Z9TtJ/As8DTcA/RcTmvGo1M7PW5TkMRUTMBGY2Wzep6P1G4Ms72PY64Lo86+sAnXaIbCdV2vFA5R1TpR0PVN4xdcnjUWHUx8zMbMc83YeZmaVyWORA0hRJb0v6P+WupRQkHSJpjqQXJD0n6ZJy19QeknpKekrSM8nx/KDcNZWCpG6SFkn6Q7lrKQVJKyQ9K2mxpAXlrqcUJO0r6T5JLyb/Px1T7pqy8jBUDiTVAe8C90bEp8tdT3tJOgg4KCKelrQ3sBD4QkQ8X+bSdkoyS8CeEfGupB7AfOCSiHiizKW1i6TLgFqgd0ScVu562kvSCqA2Ikp7T0IZSboHeDQi7kquEt0jIv5W7rqycM8iBxExj8LVXRUhIt6MiKeT9+uBF+jCd9RHwbvJYo/k1aX/1SSpL3AqcFe5a7GWSeoN1FG4CpSIeL+rBAU4LKyNkpmBhwJPlreS9kmGbBYDbwOzIqJLHw/wL8CVwJZyF1JCATwkaWEyW0NXdxiwCvjXZLjwLkl7lruorBwWlpmkvYD7gUsjYl2562mPiNgcEUMozA4wXFKXHS6UdBrwdkQsLHctJTYyImoozFz9T8nwblfWHagBbouIocB7wHaPbuisHBaWSTK2fz8wNSJ+U+56SiUZBniEwlT4XdVIYGwyxj8NOFHSv5W3pPaLiJXJz7eBGXT9macbgcaiXux9FMKjS3BYWKrkhPDdwAsR8bNy19Nekqok7Zu87wWMAl4sb1U7LyK+GxF9I6KawiwIsyPiq2Uuq10k7ZlcTEEyVHMy0KWvLoyI/wu8LunwZNVJFGap6BJyvYN7VyXpP4DjgT6SGoHvR8Td5a2qXUYCXwOeTcb5Aa5O7tDvig4C7kkeqLUb8J8RURGXm1aQjwIzCv9OoTvw7xHxQHlLKolvAVOTK6GWA/+9zPVk5ktnzcwslYehzMwslcPCzMxSOSzMzCyVw8LMzFI5LMzMLJXDwqwdJH07mT10arlrMcuTL501awdJLwKnRMQr5a7FLE/uWZjtJEm3U5gcrl7SdyQ9nkwQ9/jWu3STCQtvTJ7LsETSt5L1R0uam0yS92AyDbxZp+WehVk7bH3mAvA+8PeIaJI0CrgoIr4k6SIK04mMSz7bH1gPzAXOiIhVksYBoyPi/DIdhlkqT/dhVhr7UJhCZACFqbV7JOtHAbdHRBNARKxJZrj9NDArmc6iG/Bmx5dslp3Dwqw0fgjMiYgzk2d+PJKsF9s/WEnAcxHRZR6paeZzFmalsQ/wRvL+vKL1DwEXSuoOkAxDLQWqtj5/WVIPSZ/qwFrN2sxhYVYaNwD/LOkxCsNKW90FvAYskfQMcE5EvA+cBfw4WbcYGNHRBZu1hU9wm5lZKvcszMwslcPCzMxSOSzMzCyVw8LMzFI5LMzMLJXDwszMUjkszMwslcPCzMxS/T8kQFR1gZlq/gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.bar(np.arange(1,7),1/6.,label=\"fair dice\")\n", "plt.bar(np.arange(1,7),pn,label=\"loaded dice\",alpha=0.8)\n", "plt.xlabel(\"face\")\n", "plt.ylabel(\"probability\")\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check normalization constraint\n", "p1(x0)+p2(x0)+p3(x0)+p4(x0)+p5(x0)+p6(x0)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.999999999999999" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check mean constraint\n", "1*p1(x0)+2*p2(x0)+3*p3(x0)+4*p4(x0)+5*p5(x0)+6*p6(x0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cosmic web example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cosmic web is described by four possible structures ($\\theta_0=$ void, $\\theta_1=$ sheet, $\\theta_2=$ filament, and $\\theta_3=$ cluster. We are looking for the prior probabilities $p_n$ ($n = 0, ..., 3$) that one should assign to each of these structures. Let us assume that we have measured in astronomical observations the (average) volume of structures:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Volume of the different structures (arbitrary normalization)\n", "# (based on numbers given in table II in Leclercq et al. 2015a, arXiv:1502.02690)\n", "V0 = 60\n", "V1 = 20\n", "V2 = 8\n", "V3 = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cosmological principle states that there cannot be cosmological structures at arbitrary large scales (the \"End of Greatness\"); let us say that this imposes the average volume of a random cosmic web element to be:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# The average volume of a random cosmic web structure\n", "V = 30" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just on the basis of these observations, the maximum entropy principle fixes the new prior probabilities $p_n$ that one should use. They are given in canonical notations as $p_n = \\mathrm{e}^{-\\beta V_n}/Z$ with $Z=\\sum_n \\mathrm{e}^{-\\beta V_n}$." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# x=exp(-beta) in the usual thermodynamic notations\n", "def Z(x):\n", " return pow(x,V0)+pow(x,V1)+pow(x,V2)+pow(x,V3)\n", "def p0(x):\n", " return pow(x,V0)/Z(x)\n", "def p1(x):\n", " return pow(x,V1)/Z(x)\n", "def p2(x):\n", " return pow(x,V2)/Z(x)\n", "def p3(x):\n", " return pow(x,V3)/Z(x)\n", "def dlnZdx(x):\n", " return p0(x)*V0+p1(x)*V1+p2(x)*V2+p3(x)*V3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solving $\\left\\langle V \\right\\rangle = -\\partial Z/\\partial \\beta = \\bar{V}$ fixes $\\beta$:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0134673199344182" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def f(x):\n", " return dlnZdx(x)-V\n", "x0=fsolve(f,0.9)[0]\n", "x0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solution is therefore:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.39392560584840686,\n", " 0.23068759870801528,\n", " 0.19647498278472048,\n", " 0.17891181265885733)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P0=p0(x0); P1=p1(x0); P2=p2(x0); P3=p3(x0)\n", "P0,P1,P2,P3" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check normalization constraint\n", "P0+P1+P2+P3" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "30.000000000001336" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check average volume constraint\n", "P0*V0+P1*V1+P2*V2+P3*V3" ] } ], "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.7.3" } }, "nbformat": 4, "nbformat_minor": 1 }