Distribution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook I show how to generate a set of numbers under hard constraints on their distributional properties (namely mean and variance)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting dist.lp\n"
     ]
    }
   ],
   "source": [
    "%%file dist.lp\n",
    "\n",
    "#const n=100.\n",
    "\n",
    "lo(0).\n",
    "hi(10).\n",
    "avg(5).\n",
    "var(2).\n",
    "\n",
    "n(n).\n",
    "\n",
    "sample(0..N-1) :- n(N).\n",
    "1 { val(I,Lo..Hi) } 1 :- sample(I), lo(Lo), hi(Hi).\n",
    "    \n",
    ":- avg(A), n(N), not A*N = #sum { V,I:val(I,V) }.\n",
    ":- avg(A), var(SS), n(N), not SS*N = #sum { (V-Mean)**2,I:val(I,V),Mean=A}.\n",
    "    \n",
    "#show val/2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture out\n",
    "n = 100\n",
    "!clingo dist.lp -c n={n} --outf=2 --sign-def=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.0, 1.4142135623730951, 2.0, 1.0, 10.0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFRRJREFUeJzt3V+MXOd53/HvI5EEp6JJicFKhiN7x27i0LFNy7xw0rqF\nh7YVKw4iJUBb2EFbx2mBGDAjIS2aSM2FeFXUAdLGbUoUaTaqbGQpwXYCy4CDyIK0F0Zgl65MrWNR\njlJnKSmJuZv6D0qDsBj56cXMUuPZ2dmZOWd2Zt79foAFZ95555znvHP4m7Nnds4bmYkkqSzXTbsA\nSVL9DHdJKpDhLkkFMtwlqUCGuyQVyHCXpALtGO4RsRQRlyJitavtNyPiQkScj4hPRcThyZYpSRrF\nMEfuDwDv6Wl7FHhjZt4GPAvcV3dhkqTx7Rjumfl54Fs9bY9l5vc7d78A3DqB2iRJY6rjnPsvAX9c\nw3IkSTWpFO4R8RvA1cxcrqkeSVIN9o37xIj4APBe4J079PPiNZI0hsyMcZ877JF7dH7adyLuAH4N\nuDMzv7fTkzNzbn/uv//+qdewV+uf59qtf/o/815/VcP8KeQy8KfA6yPiuYj4IPBfgUPA5yLiyYg4\nU7kSSVJtdjwtk5m/0Kf5gQnUIkmqid9Q3UGr1Zp2CZXMc/3zXDtY/7TNe/1VRR3ndgauICInvQ5J\nKk1EkLvwgaokaY4Y7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK\nZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCG\nuyQVyHCXpAIZ7pJUoB3DPSKWIuJSRKx2td0UEY9GxNci4k8i4shky5QkjWKYI/cHgPf0tN0LPJaZ\nPwY8DtxXd2GSpPHtGO6Z+XngWz3NdwEPdm4/CPxczXVpj9nY2ODcuXNsbGxMuxSpCOOec785My8B\nZOY3gIX6StJec/bswywuHuP22z/E4uIxzp59eNolSXMvMnPnThGLwGcy83jn/jcz82jX4/83M39o\nm+fmMOvQ3rSxscHi4jGuXHkCOA6s0mic5OLFZ1hY8JhBe1dEkJkx7vP3jfm8SxFxS2ZeiohXAuuD\nOp8+ffra7VarRavVGnO1Ks3a2hoHDjS5cuV4p+U4+/cvsra2ZrhrT1lZWWFlZaW25Q175N6kfeT+\n5s79jwDfzMyPRMSvAzdl5r3bPNcjd23LI3epv6pH7sP8KeQy8KfA6yPiuYj4IPAfgdsj4mvAuzv3\npZEtLCywtHSGRuMkhw+foNE4ydLSGYNdqmioI/dKK/DIXUPY2NhgbW2NZrNpsEtUP3I33CVpBk38\ntIwkaf4Y7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDfY45e5Gk7Rju\nc8rZiyQN4oXD5pDXQJfK54XD9qDN2YvawQ7dsxdJEhjuc6nZbPLii2vAaqdllatXL9JsNqdXlKSZ\nYrjPIWcvkrQTz7nPMWcvksrlTEySVCA/UJUkbWG4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ\n7pJUIMNdkgpkuEtSgSqFe0T8akT8WUSsRsQfRMSBugqTJI1v7HCPiFcBvwKcyMzjwD7gfXUVprLU\nPSXgPE0xOGqtvf2rPn83zNPrsWdk5lg/wKuAi8BNtIP9M8C7+/RL7W3Lyw9lo3E0jxw5kY3G0Vxe\nfmimljdJo9ba2//UqXsqPX83xmaeXo950snO8TO60pPhbuD/AZeAj2/TZ7IjoJm2vr6ejcbRhKcS\nMuGpbDSO5vr6+kwsb5JGrXVr/ycSGhWeP/mxmafXY95UDfd94x7xR8SNwF3AIvAd4JMR8QuZudzb\n9/Tp09dut1otWq3WuKvVnNmcEvDKla1TAo5zDfq6lzdJo9a6tf8NwKvpN53icM+f/NjM0+sx61ZW\nVlhZWalvgeO+KwD/BPgfXff/BfA7ffpN9N1Ns80jd4/cNR6mdVoGeBvwFeAgEMD/BD7cp9+kx0Az\nbvOc7OHDb631nHtdy5ukUWvt7X/q1N2Vnr+b59zn4fWYJ1XDvdJMTBFxP+2/kLkKfBn415l5tadP\nVlmHylD3lIDzNMXgqLX29q/6/N0wT6/HvHCaPUkqkNPsSZK2MNwlqUCGuyQVyHCXpAIZ7pJUIMNd\nkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWp\nQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAJVCveIOBIRn4iI\nCxHx1Yj4iboKkySNr+qR+0eBz2bmG4C3ABeql6RJ2djY4Ny5c2xsbOyJOia9njqWPyuviQqUmWP9\nAK8A/s8Q/VLTt7z8UDYaR/PIkRPZaBzN5eWHiq5j0uupY/mz8ppoNnWyc/yMHvuJ7SP1LwIPAE8C\nvws0+vSb9BhoB+vr69loHE14KiETnspG42iur68XWcek11PH8mflNdHsqhru+yoc9O8DTgAfzswv\nRcRvA/cC9/d2PH369LXbrVaLVqtVYbUa1draGgcONLly5Xin5Tj79y+ytrbGwsJCcXVMej11LH9W\nXhPNjpWVFVZWVupb4LjvCsAtwNe77v8j4DN9+k303U07m5WjRI/cd69GzT8qHrmP/YFqZl4Cno+I\n13ea3gU8XeWNRpOxsLDA0tIZGo2THD58gkbjJEtLZ3b9CHG36pj0eupY/qy8JipXtN8gxnxyxFuA\n3wP2A18HPpiZ3+npk1XWofpsbGywtrZGs9mcaojsVh2TXk8dy5+V10SzJyLIzBj7+ZMOXsNdkkZX\nNdz9hqokFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklSgKtdz\nL97mRZ0OHTrE5cuXaTabAD/QNuixUfuPuqx+F5vqvhDVoP67vW071dpb105tVV/LOrZtUD2j1jyp\ni5ANahtmvxjUp99Y7LT/9bbVvY/11rOnVble8DA/zOn13DenQGs0XpfQyEbjzbl//yvywIEj19r2\n73/Nto+N2n/UZTUab94yNVv3tG2D+u/2tu1Ua29dO7VVfS3r2LZB9Yxa86Sm7BvUNsx+0b3d/ZbZ\nOxanTt0zcP+b9D42zGszT5jWNHtDr2AOw/3liRSeSNicUGE94aautkGPjdp/1GVtneDhByd/2L7/\nwYM37vK27VRrb12D28afFKPObdu+nlEn4ZjUxB+Dx3OY/WJQn35j8URCo8Z9Ztz+5Ux+UjXcPefe\nx+YUaHAD0ASOA2vAa7vaBj02av9Rl7V1araXax7c//rrb+a66169i9u2U629dQ1u21zGsLa+lnVs\n2/b19Nu2QTWP2n/YZQwez2H2i0F9+o3FDUDvY7u5j/XWM/o4lsZw76PZbPLii2vAd2nvNKu0d6K/\n7Gob9Nio/Udd1mqn0lWuXr1Is9nsqnlw/5deWuf7339+F7dtp1p76xrctrmMYW19LevYtu3r6bdt\ng2oetf+wyxg8nsPsF4P69BuL7wK9j+3mPtZbz+jjWJwqh/3D/DCHp2UyXz6HefBgM9vn9N6U+/cf\nygMHjlxr27//1ds+Nmr/UZfVaLxp2/PYhw+/dWD/3d62nWrtrWuntqqvZR3bNqieUWuucxuHHc9h\n9ovu7e63zN6xOHXq7oH736T3sWFem3lCxdMyzsQ0gH8t41/L+Ncy/rXMtDjNniQVyGn2JElbGO6S\nVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAlUO94i4LiKejIhH6ihIklRdHUfu9wBP\n17AcVbCxscG5c+fY2NiYdilTN82xGGbdw9Y3qF/3Y3WuUwWpctUx4Fbgc0ALeGSbPpO5ZJquqWMm\nn1JMcyyGWfew9Q3qN2g2pDpmhtJsYJozMQGfAG4D3mG4T0cdM/mUYppjMcy6h61vUL/+M27VNzOU\nZkfVcB97guyI+BngUmaej4gWsO3Vy06fPn3tdqvVotVqjbta9dichefKla0z0Mz7JU9HNc2xGGbd\nw9Y3qB/Q9dg5tpt9aNR1avpWVlZYWVmpb4HjvisA/wF4Dvg68DfAZeBjffpN+P1tb/PI7GUeuXvk\nXhJmYYJsPC0zVXXM5FOKaY7FMOsetr5B/QbNhlTHzFCaDVXDvZbJOiLiHcC/zcw7+zyWdaxDg9Ux\nk08ppjkWw6x72PoG9Rs0G1KVdWp2OBOTJBXImZgkSVsY7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLc\nJalAhrskFchwl6QCGe6SVCDDXZIKZLhrJjktnFSN4a6Zc/bswywuHuP22z/E4uIxzp59eNolSXPH\nq0JqpmxsbLC4eIwrV56gPcPQKo3GSS5efMZL1WpP8aqQKsrmtHD9po6TNDzDXTOl2Wzy4otrwGqn\nZZWrVy9em5RC0nAMd82UhYUFlpbO0Gic5PDhEzQaJ1laOuMpGWlEnnPXTHJaOO11TrMnSQXyA1VJ\n0haGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAY4d7RNwaEY9HxNMR8ZWIuLvO\nwiRJ4xv78gMR8UrglZl5PiIOAf8buCszn+np5+UHJGlEU7v8QGZ+IzPPd25fBi4APzzu8qSd9Jt6\nb7PtwoULU5uWbzenBJzEuoZdplMfzpnMrPwDNIE14FCfx1Kqann5oWw0juaRIyey0Tiay8sPXWtr\nNF6X0MhG483XHptmXfO0rmGXuZvbqbZOdo6fy1We3F4/h4Av0T4l0+/xiQ6Ayre+vp6NxtGEpxIy\n4ak8ePDGTtsTCT/4WKNxNNfX16dS16TWPYl1DbvM3dxOvaxquO+rctQfEfuATwIfz8xPb9fv9OnT\n1263Wi1arVaV1WqP2Zx678qVl6feu/76m4EGcAPtXxy3Tss36evA96trUuuexLqGXeZubudetrKy\nwsrKSn0LrPLOAHwM+E879Jnge5v2Ao/cPXLfi5jWaRng7cBLwHngy8CTwB19+k16DLQHbJ7zPXz4\nrVvOuR882Oycc3/T1M65d9c1T+sadpm7uZ1qqxruzsSkudFv6r3NtkOHDnH58uWpTMu3m1MCTmJd\nwy7TqQ93l9PsSVKBnGZPkrSF4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEu\nSQUy3KVd4kxG2k2Gu7QLzp59mMXFY9x++4dYXDzG2bMPT7skFc4Lh0kTtrGxweLiMa5ceYL2pCKr\nNBonuXjxGa+uqG154TBpxm3OZNRvtihpUgx3acKazSYvvrgGrHZaVrl69SLNZnN6Ral4hrs0YQsL\nCywtnaHROMnhwydoNE6ytHTGUzKaKM+5S7vEmYw0CmdikqQC+YGqJGkLw12SCmS4S1KBDHdJKpDh\nLkkFMtwlqUCGuyQVyHCXpAIZ7pJUoErhHhF3RMQzEfHnEfHrdRUlSapm7HCPiOuA3wHeA7wReH9E\nHKursFmxsrIy7RIqmef657l2sP5pm/f6q6py5P424NnMvJiZV4GHgLvqKWt2zPsOMs/1z3PtYP3T\nNu/1V1Ul3H8YeL7r/gudNknSlFUJ935XK/Pyj5I0A8a+5G9E/CRwOjPv6Ny/F8jM/EhPPwNfksYw\nleu5R8T1wNeAdwF/A/wv4P2ZeWHcYiRJ9dg37hMz86WIOAU8Svv0zpLBLkmzYeIzMUmSdt/EvqEa\nEb8ZERci4nxEfCoiDnc9dl9EPNt5/KcmVUMV8/YFrYi4NSIej4inI+IrEXF3p/2miHg0Ir4WEX8S\nEUemXesgEXFdRDwZEY907jcj4gud+s9GxNi/bU5aRByJiE909uuvRsRPzMv4R8SvRsSfRcRqRPxB\nRByY9bGPiKWIuBQRq11t2453RPyXTu6cj4jbplP1tVr61V5rZk7y8gOPAm/MzNuAZ4H7OkX+OPDP\ngDcAPw2ciYixPzSYhDn9gtbfAf8mM38c+AfAhzs13ws8lpk/BjxO53WYYfcAT3fd/wjwW536vw38\nq6lUNZyPAp/NzDcAbwGeYQ7GPyJeBfwKcCIzj9M+Xft+Zn/sH6D9f7Rb3/GOiJ8G/n5m/ijwy8B/\n381C++hXe62ZObFwz8zHMvP7nbtfAG7t3L4TeCgz/y4z12hvxNsmVceY5u4LWpn5jcw837l9GbhA\ne8zvAh7sdHsQ+LnpVLiziLgVeC/we13N7wQ+1bn9IPDzu13XMCLiFcA/zswHADr793eYn/G/Hrih\nc3TeAP4aOMkMj31mfh74Vk9z73jf1dX+sc7zvggciYhbdqPOfvrVXndm7taFw34J+Gzndu+Xn/6K\n2fvy01x/QSsimsBttHeQWzLzErTfAICF6VW2o/8M/Ds635eIiB8CvtW1w78AvGpKte3kdcDfRsQD\nndNKvxsRf485GP/M/Gvgt4DnaP9//A7wJPDtORn7bjf3jPfNnfZ5yJ1ulTOz6oXDPtc5R7f585XO\nvz/b1ec3gKuZeXazqc+iZu1T3Xmosa+IOAR8ErincwQ/L3X/DHCp89vH5vgHW1+LWd2efcAJ4L9l\n5gngu7RPEcxqvddExI20j2wXaQf4DbR//e8189sywNz8n64rMyt9QJKZtw96PCI+QPvX7Hd2Nb8A\nvLrr/q20fwWcJS8Ar+m6P4s1btH5lfqTwMcz89Od5ksRcUtmXoqIVwLr06twoLcDd0bEe2mfFngF\n8Nu0f32+rnMEOcuvwwvA85n5pc79T9EO93kY/3cDX8/MbwJExB8B/xC4cU7Gvtt24z0PuVNrZk7y\nr2XuAH4NuDMzv9f10CPA+zqfxr8W+BHaX4CaJeeAH4mIxYg4ALyPdt2z7veBpzPzo11tjwC/2Ln9\nAeDTvU+aBZn57zPzNZn5Otrj/Xhm/nPgCeCfdrrNcv2XgOcj4vWdpncBX2U+xv854Ccj4mDng7rN\n2udh7Ht/u+se71/k5ZofAf4lXPt2/bc3T99M0Q/UXntmZuZEfmif9L9I+9zdk8CZrsfuA/6C9od+\nPzWpGirWfwftb+A+C9w77XqGqPftwEvAeeDLnTG/AzgKPNbZls8BN0671iG25R3AI53brwW+CPw5\n8DCwf9r1Daj7LbQPDM4DfwgcmZfxB+7v/H9cpf1B5P5ZH3tgmfYR7Pdov0F9ELhpu/Gm/RdwfwE8\nRfsvg2at9loz0y8xSVKBnGZPkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVKD/D0xb\nXjdKSvHLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110ba6bd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import json\n",
    "result = json.loads(out.stdout)\n",
    "value = result['Call'][0]['Witnesses'][0]['Value']\n",
    "\n",
    "sample = zeros(n)\n",
    "def val(k,v):\n",
    "    sample[k] = v\n",
    "\n",
    "env = {'val': val}\n",
    "    \n",
    "for atom in value:\n",
    "    eval(atom, env)\n",
    "    \n",
    "scatter(arange(n),sample)\n",
    "mean(sample), std(sample), var(sample), min(sample), max(sample)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}