Random Image

{
 "metadata": {
  "name": "",
  "signature": "sha256:3ef7145c7aa9b79a3d9fc1565673db017213f2d09b0e1a3d849f80a0fd20242e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Binary images generated by a mixture of multivariate Bernoulli distributions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "import scipy.stats as st\n",
      "from IPython.html.widgets import interact\n",
      "from IPython.display import display"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "   $X=(X_1, X_2, \\ldots, X_d)\\sim Bern(p_1, p_2, \\ldots, p_d)$ is a multivariate Bernoulli random variable, whose components, $X_i$, are  independent Bernoulli random variables of parameters $p_i$, $i=\\overline{1,d}$.\n",
      "\n",
      "\n",
      "$X$  is a model for a binary image  of $d$ pixels.\n",
      "We take $d=r^2$,  and define a mixture of three  multivariate Bernoulli random variables:\n",
      "\n",
      "$$Y=\\pi_1X^{(1)}+\\pi_2X^{(2)}+\\pi_3 X^{(3)}, \\quad \\pi_i\\in(0,1)\\,\\,\n",
      "\\sum_{k=1}^3\\pi_i=1,$$ \n",
      "\n",
      "\n",
      "\n",
      "$$X^{(k)}=(X^{(k)}_1,X^{(k)}_2,\\ldots, X^{(k)}_{d})\\sim \\mbox{Bern}(p_{k1},p_{k2}, \\ldots, p_{k d})$$\n",
      "The parameters $p_{ki}$ of the r.v. $X^{(k)}_i$,  $k=\\overline{1,3}$, $i=\\overline{1,d}$, are uniformly distributed \n",
      "on a subinterval of $(0,1)$.\n",
      "\n",
      "In our example, the parameters associated to the random vectors  $X^{(k)}$, $k=1,2,3$,  are   uniformly distributed respectively\n",
      "    in the intervals, $(0.2, 0.8)$, $(0.35, 0.75)$, $(0.3, 0.9)$.\n",
      "\n",
      "Simulating  this mixture, we generate n binary images, and display them  progressively with `interact`.\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def  simDiscrete(pr):\n",
      "    k=0\n",
      "    F=pr[0]\n",
      "    u=np.random.random()\n",
      "    while(u>F):\n",
      "        k+=1\n",
      "        F=F+pr[k]\n",
      "    return k"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.rcParams['figure.figsize'] = 4, 4\n",
      "\n",
      "Pi=[0.37, 0.41, 0.22]\n",
      "\n",
      "d=2500\n",
      "r=50\n",
      "p=np.zeros((3,d))\n",
      "sint=[(0.2, 0.8), (0.35, 0.75), (0.3, 0.9)]# subintervals for parameters\n",
      "for k in range(3):\n",
      "    p[k]=st.uniform.rvs(size=d, loc=sint[k][0], scale=sint[k][1]-sint[k][0] )\n",
      "\n",
      "\n",
      "def  RandomImage(n): \n",
      "        fig, ax=plt.subplots()\n",
      "        k=simDiscrete(Pi)\n",
      "        pixels=map(lambda (pr): st.bernoulli.rvs(pr), p[k])     \n",
      "        img=np.array(pixels).reshape((r,r)) \n",
      "        ax.imshow(img, origin='upper',  cmap='gray', interpolation='nearest')           "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "interact(RandomImage, n=(1,20))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<function __main__.RandomImage>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD+CAYAAADxoQNSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFFFJREFUeJzt3U1sE0cbB/C/I5A4FBG1VZaKUKUfxBBCcVRabiXBcnIA\nQtRIERwiBFx7gEM/6IWoBzDlUAW1x7ayOEBzogi1qaiwaRHqh1AjSqFKD0G4qW2JhrQCQYuSeQ+v\n8JuYjXc8mZldv/P/SZZI2J15du0nuzOemY0JIQSIyBkNYQdARHYx6Ykcw6QncgyTnsgxTHoixzDp\niRyzqKQfHR3F2rVrsWbNGhw7dkxXTERkUEz1e/qZmRnE43F8/fXXWLVqFV555RWcOnUK69at0x0j\nEWmkfKX/4Ycf8OKLL6KlpQVLly7Frl278Pnnn+uMjYgMWKK64+TkJFavXl3+ubm5Gd9///28bTo7\nO3Hx4kX16IhI2UI38cpX+lgsFrjNxYsXIYTA4cOHFwygkhBC6WWqHJkYVY9LxUJxPzrHJo9L93sT\nFLOO96qWc6ijbFN0xqec9KtWrUI+ny//nM/n0dzc/Nh2Q0NDyOVyGBoaUq2KiHQSih4+fCief/55\nMTExIf755x+xceNGcf369XnbAJj3sqmybr9XNYcPH5Yu22TMsqrFK1uXyktXXTr2UaUaj0rZsvWb\nfG+U2/RLlizBhx9+iJ6eHszMzGD//v3/Vz33nZ2dYYdQk3qLl8Kj/JWdVOEV7X6DVQXW7Uc1HlPH\n5RezqXMmc35kyMSn8l6YfP8qyZx31fdG5rOi672otFB8HJFH5BgmPZFjlNv0subeYui6ZbN5G+xH\n5VZUF13HrnKbqXqObTYBdN2Wy1Bp5ulqJiwGr/REjmHSEzmGSU/kGONt+rlMte1Ut1Gl0r6SaUOH\nvY0Km+1aXWx+lWzzuGTxSk/kGCY9kWOY9ESOYdITOcZqR17UxlPXA5uDTYKoDoYJ2kd1P5sduro6\nKG0O5FoIr/REjmHSEzmGSU/kmNAH56i0cUxOWghzUoeu9p5qPKb6BmyeL5sDwHTso1pO0LFXO0Ze\n6Ykcw6QncgyTnsgxTHoix4Q+OMdW55Ff/aqdW7pmtZlaWNHUQBJd50KVTDmmBr/YnP1pegAPr/RE\njmHSEzmGSU/kmNAH58jQ1WbVUbds/bomWpiasKHSPtd1DDJ0DYpS3cbm5DCVclVjBnilJ3IOk57I\nMUx6Iscw6YkcY7wjr9bOnyjNVJJl8wmwYa7+Y/Oc2lzq2+aS4TJ0dJhylh0RlTHpiRzDpCdyjNVH\nVfsxtXKOzXZ21MqO2mo/NgdXmfocmIxHZpVfnRObeKUncgyTnsgxgUm/b98+eJ6HDRs2lH83NTWF\nVCqF1tZWdHd3Y3p62miQRKRPYNLv3bsXo6Oj836XTqeRSqUwPj6OZDKJdDptLEAi0kxImJiYEO3t\n7eWf4/G4KBaLQgghCoWCiMfjvvtVFg9A6aWLTLmm4lEtV9e5sHVcuraRiS/sc6oSj+p5r3Wfatso\ntelLpRI8zwMAeJ6HUqmkUgwRhWDRX9nFYrGqXxUMDQ0ttgoiCpDL5ZDL5eQ2lrm18Lu9LxQKQggh\n/vjjD97eK9QtW66uc2HruHRtIxNf2OdUJR7V817rPtW2Ubq97+3tRSaTAQBkMhn09fUtuO2jO4FY\nLAYhhNJrbhmyLz+V5apuI0MmnqB9VPfzI3NcOs6p334q28h8DlSpvMcy50P1860Ss8p7VT4WEVDr\n7t27cfHiRdy+fRue5+G9997Dzp07MTAwgFu3bqGlpQUjIyNobGz0PVGVgatQeYMXk7BBdcuUbWoZ\nY7+6bY14Uy23cj+b5zTsz4Epi1nOKzDpF4NJrx+TvjZhfw5MWUzSc0QekWOY9ESOsXp7L8Pm7ZjJ\nxwfZvKU1dVxh3r6q0rW6TtSaPzJ1yTbzeKUncgyTnsgxTHoix4S+co6MMPsGTLW7/fbz20alfl1t\nVJkydPVdqNSvWq7KfjLnVFe/jemvhHmlJ3IMk57IMUx6Iscw6YkcY/X59Kp0DT5ZzHhl3fv47adr\nmWXVbSqZGrykq0Mw7AFYQXWrbqOjrmrHzSs9kWOY9ESOYdITOSZyj6qWoTpQwhSTA1JU2qi6+gZ0\nTQQxRXUQkq4BRTbXVdDZV8ErPZFjmPREjmHSEzmGSU/kmLqcZVcPgzJUZk6pdkyFOQtRhs34THU+\nmhw4pWu1n2rbz8UrPZFjmPREjmHSEznG6uAck4MpgvaRFbX2sc0Vc3X0i4S96m/QACNVulbVVV0d\nSednjFd6Iscw6Ykcw6QncgyTnsgxdTE4R4VqZ5/NpbPDpNLpZHKgi6mlolXp6sRU2c/0o7h4pSdy\nDJOeyDFMeiLHRG41XF2TV2TL1sHmyq26Hn0lU7au825ytR8dTE6m0YWDc4hIWWDS5/N5dHV1Yf36\n9Whvb8eJEycAAFNTU0ilUmhtbUV3dzemp6eNB0tEixcTAfcNxWIRxWIRiUQCd+/excsvv4wzZ87g\n008/xdNPP4233noLx44dw507d5BOp+cXHovVPO9X1+29ycUqVeoyeXsvI+q39ypMPQ037Pp1neeF\nygm80q9cuRKJRAIA8MQTT2DdunWYnJzE2bNnsWfPHgDAnj17cObMGaVAiciuwCv9XDdv3sSWLVtw\n7do1PPvss7hz5w6A//5FefLJJ8s/lws3NJhC119mXXcDNpeGDnOAjGrdNp+/ruv9C/MuzFTH7CPS\nHXl3795Ff38/hoeHsXz58scCCLP3lYjkSSX9w4cP0d/fj8HBQfT19QEAPM9DsVgEABQKBTQ1NQWW\nk8vl1CMlIi0Ck14Igf3796OtrQ0HDhwo/763txeZTAYAkMlkyn8Mquns7FSPlIi0CGzTX7p0Ca+9\n9hpeeuml8i380aNH8eqrr2JgYAC3bt1CS0sLRkZG0NjYOL9wQz3zuth8ZLLNeGTL1kFX2zfM1X9M\ntsVN9cHI9L8stE1NHXm1YtKHE49s2Tow6auXE8Wk54g8Iscw6Ykcw6QnckzoK+fU42AKFbpWVdEl\nzJVzZPbzq0vXo8KCylUtO2pDvxfCKz2RY5j0RI5h0hM5xupjrfyotHl0rRZq87FIuuh6dJKu1Wx0\nPQYsanM3dK4+G7SP7f4oXumJHMOkJ3IMk57IMUx6IsdEbglsk0ytSGJyooUKXZ1iYa90VMnU+Qp7\nQJGpDsGF8EpP5BgmPZFjmPREjgl9wo2KKAxwqDUe2f2Cygn7cdEy5aiUa2ohEJsLt+gaECZjMf0k\nvNITOYZJT+QYJj2RY0L/nt7UwggmJ3Xo+l61ks1FIVTOj80nuMjUH/aipTKiNrYA4JWeyDlMeiLH\nMOmJHMOkJ3JM6B15lUyuKKoy0EW1HJtPbKnHgUAq+9h8kk/Un7wku58fXumJHMOkJ3IMk57IMVaf\nWuvHZns0zKeemhzcYevpwDafAKtatspx6mpD23xqbVA5fGotEZUx6Ykcw6QncgyTnsgxdblyjq6O\nKl0dL6qdNaY6nUwN8rE5y87mije6hNlxXQte6YkcUzXpHzx4gM2bNyORSKCtrQ2HDh0CAExNTSGV\nSqG1tRXd3d2Ynp62EiwRLV7VpF+2bBmy2SzGxsZw9epVZLNZXLp0Cel0GqlUCuPj40gmk0in07bi\nJaLFEpLu3bsnNm3aJK5duybi8bgoFotCCCEKhYKIx+O++1QWD0DpJUNlH5Vydcaj69hNxaNSrswx\n6NpGV3y6Xqps50Bgm352dhaJRAKe56Grqwvr169HqVSC53kAAM/zUCqVgoohoogI7L1vaGjA2NgY\n/vrrL/T09CCbzc77/1gspm06IxGZJ/2V3YoVK7Bt2zZcuXIFnuehWCxi5cqVKBQKaGpqWnC/oaEh\nHXESURW5XE5626oTbm7fvo0lS5agsbER9+/fR09PDw4fPoyvvvoKTz31FN5++22k02lMT0/7duZV\nDvoPe9KJinqYjCHD1KScqC0WUsnmXajqe2V7slHVpP/555+xZ88ezM7OYnZ2FoODg3jzzTcxNTWF\ngYEB3Lp1Cy0tLRgZGUFjY6NvELWeCJOPNQ5zVpTNeEx90FX/COhKVlvHqUrXe6zjD3213DM+tZZJ\nbz8eJn046iXpOSKPyDFMeiLHWL29N7nSrWw8tbLZORP2SrdBdZns1NTFVOeoyaaEqXh4e09EAJj0\nRM5h0hM5hklP5JjQH2ulq5PHZseZSt269jM5jkHHPiY7vHStwGOqLplydNWlOqYD4JWeyDlMeiLH\nMOmJHGO8Ta+jjWdzNVVd4+FVhN1GDdrG1GOgdZZjqt+mHuKRndHKKz2RY5j0RI5h0hM5hklP5JjQ\nB+eodF6pDHhQFfZiDrpmBuo4z/W4BJnqNlGb/anzM8crPZFjmPREjmHSEzkm9Da9Lrra/Srl2Gz/\n6RosFOYKMzaPwSaTi17qHDTGKz2RY5j0RI5h0hM5hklP5BjjHXm1PsvO5io0usoJ+9lspjofdQ0M\nMrX0uK5zEfbS7JVU4+EsOyLyxaQncgyTnsgxdTE4x9QEFxMrljxi6mmltp56KsvUZCiTg6JkYrY5\nccfm05QBXumJnMOkJ3IMk57IMUx6IsdYfT697D5BbD76SrUcGfVWl8mVc1SE/Zx7Gbre41oHe1XL\nPV7piRwjlfQzMzPo6OjAjh07AABTU1NIpVJobW1Fd3c3pqenjQZJRPpIJf3w8DDa2trKtxTpdBqp\nVArj4+NIJpNIp9NGgyQijUSAfD4vksmkuHDhgti+fbsQQoh4PC6KxaIQQohCoSDi8bjvvhLFCwBa\nXjLl6qpbV10qbNZls26V864rZl1MHnut56da3YFX+oMHD+L48eNoaPjfpqVSCZ7nAQA8z0OpVAoq\nhogiouow3HPnzqGpqQkdHR3I5XK+28Risaq9j0NDQ+V/d3Z2orOzUyVOIgowN9eqqfqV3bvvvouT\nJ09iyZIlePDgAf7++2+8/vrr+PHHH5HL5bBy5UoUCgV0dXXh119/fbxwia/sTH01YvIhFbrqCjo3\nfsJ8MqrJum0u1KnjXOisy0QOKH9ld+TIEeTzeUxMTOD06dPYunUrTp48id7eXmQyGQBAJpNBX1+f\nlqCJyLyaZtk9+ov0zjvvYGBgAB9//DFaWlowMjJS0/5zqVx9Tf2lVi3bbx9TM6dk4tMVj81zIfM5\nMHW+TH4uTc2AVP08AZZH5IW91JOt22C/snUlvc14TN26m7zltdlkkqlb53r1OuoCOCKPyDlMeiLH\nhD7hJsxVX1Tpul0Ns7nhx9a5D7tvx9TqvKp16/pWKajcR3ilJ3IMk57IMUx6Iscw6YkcY3UJ7KgN\nibQ5NFZ1UI1MXSrl6PqeXqVcP1Hr0LX5OTC1cs5CeKUncgyTnsgxTHoix4T+WCtTA110taFl9tP1\neG2bEzbCfHy0Lrr6baL2/pkew88rPZFjmPREjmHSEzmGSU/kmNA78nQxtWiFyYUjwlzXT4Wu2V82\nHxllcxUh1c+KyU5nP7zSEzmGSU/kGCY9kWOMt+lrbb/pahepto/DfESxzbadyYk7KmyuIqSrrjDX\n4eeEGyKSxqQncgyTnsgxTHoixxjvyJPtXPDbvpZtdHXAqQh7CeyorToTxOSAmTA7R2VEYaYir/RE\njmHSEzmGSU/kGKuDc2w+Gkh2v6ByVFZy9SvHZD+Eyccf18rm4CqZcvyY+hzIiEI/BK/0RI5h0hM5\nhklP5BgmPZFjjCe9EALZbFa64yEWiwW+hBCPvWT284ut8vVILpdbcJugcmWPy69s1boexbtQXSpk\n4q2k+t7oiNcvZj8y532hfR59lnUOlpE5X0HHWcs5s3Kln/uBrBf1FnO9xVuv/h/OM2/viRzDpCdy\njTBoy5YtAgBffPFl+bVly5YF8zImdPZIEFHk8faeyDFMeiLHGE360dFRrF27FmvWrMGxY8dMVqVs\n37598DwPGzZsKP9uamoKqVQKra2t6O7uxvT0dIgRPi6fz6Orqwvr169He3s7Tpw4ASC6cT948ACb\nN29GIpFAW1sbDh06BCC68c41MzODjo4O7NixA0B9xBzEWNLPzMzgjTfewOjoKK5fv45Tp07hxo0b\npqpTtnfvXoyOjs77XTqdRiqVwvj4OJLJJNLpdEjR+Vu6dCk++OAD/PLLL/juu+/w0Ucf4caNG5GN\ne9myZchmsxgbG8PVq1eRzWZx6dKlyMY71/DwMNra2sqDX+oh5kCmeu4vX74senp6yj8fPXpUHD16\n1FR1izIxMSHa29vLP8fjcVEsFoUQQhQKBRGPx8MKTcrOnTvF+fPn6yLue/fuiU2bNolr165FPt58\nPi+SyaS4cOGC2L59uxCi/j4bfoxd6ScnJ7F69eryz83NzZicnDRVnValUgme5wEAPM9DqVQKOaKF\n3bx5Ez/99BM2b94c6bhnZ2eRSCTgeV65aRLleAHg4MGDOH78OBoa/pcmUY9ZhrGkD3OhSp0WOxbc\npLt376K/vx/Dw8NYvnz5vP+LWtwNDQ0YGxvD77//jm+++QbZbHbe/0ct3nPnzqGpqQkdHR0LjrOP\nWsyyjCX9qlWrkM/nyz/n83k0Nzebqk4rz/NQLBYBAIVCAU1NTSFH9LiHDx+iv78fg4OD6OvrA1Af\nca9YsQLbtm3DlStXIh3v5cuXcfbsWTz33HPYvXs3Lly4gMHBwUjHLMtY0m/atAm//fYbbt68iX//\n/RefffYZent7TVWnVW9vLzKZDAAgk8mUkyoqhBDYv38/2tracODAgfLvoxr37du3y73c9+/fx/nz\n59HR0RHZeAHgyJEjyOfzmJiYwOnTp7F161acPHky0jFLM9lh8MUXX4jW1lbxwgsviCNHjpisStmu\nXbvEM888I5YuXSqam5vFJ598Iv7880+RTCbFmjVrRCqVEnfu3Ak7zHm+/fZbEYvFxMaNG0UikRCJ\nREJ8+eWXkY376tWroqOjQ2zcuFFs2LBBvP/++0IIEdl4K+VyObFjxw4hRP3EXA2H4RI5hiPyiBzD\npCdyDJOeyDFMeiLHMOmJHMOkJ3IMk57IMUx6Isf8B3+vKIjgD6JgAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xa8d3860>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "def  css_styling():\n",
      "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
      "    return HTML(styles)\n",
      "css_styling()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<style>\n",
        "    @font-face {\n",
        "        font-family: \"Computer Modern\";\n",
        "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
        "    }\n",
        "    div.cell{\n",
        "        width:800px;\n",
        "        margin-left:16% !important;\n",
        "        margin-right:auto;\n",
        "    }\n",
        "    h1 {\n",
        "        font-family: Helvetica, serif;\n",
        "    }\n",
        "    h4{\n",
        "        margin-top:12px;\n",
        "        margin-bottom: 3px;\n",
        "       }\n",
        "    div.text_cell_render{\n",
        "        font-family: Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
        "        line-height: 145%;\n",
        "        font-size: 130%;\n",
        "        width:800px;\n",
        "        margin-left:auto;\n",
        "        margin-right:auto;\n",
        "    }\n",
        "    .CodeMirror{\n",
        "            font-family: \"Source Code Pro\", source-code-pro,Consolas, monospace;\n",
        "    }\n",
        "    .prompt{\n",
        "        display: None;\n",
        "    }\n",
        "    .text_cell_render h5 {\n",
        "        font-weight: 300;\n",
        "        font-size: 22pt;\n",
        "        color: #4057A1;\n",
        "        font-style: italic;\n",
        "        margin-bottom: .5em;\n",
        "        margin-top: 0.5em;\n",
        "        display: block;\n",
        "    }\n",
        "    \n",
        "    .warning{\n",
        "        color: rgb( 240, 20, 20 )\n",
        "        }  \n",
        "</style>\n",
        "<script>\n",
        "    MathJax.Hub.Config({\n",
        "                        TeX: {\n",
        "                           extensions: [\"AMSmath.js\"]\n",
        "                           },\n",
        "                tex2jax: {\n",
        "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
        "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
        "                },\n",
        "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
        "                \"HTML-CSS\": {\n",
        "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
        "                }\n",
        "        });\n",
        "</script>\n"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<IPython.core.display.HTML at 0xaa4ba58>"
       ]
      }
     ],
     "prompt_number": 6
    }
   ],
   "metadata": {}
  }
 ]
}