MMD tests.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "# this is at git commit: 6ef795ba6e83706b903dc84a43a36e496597c394 \n",
    "from modshogun import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# use scipy for generating samples\n",
    "from scipy.stats import norm, laplace\n",
    "\n",
    "def sample_gaussian_vs_laplace(n=220, mu=0.0, sigma2=1, b=np.sqrt(0.5)):    \n",
    "    # sample from both distributions\n",
    "    X=norm.rvs(size=n, loc=mu, scale=sigma2)\n",
    "    Y=laplace.rvs(size=n, loc=mu, scale=b)\n",
    "    \n",
    "    return X,Y\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(220,) (220,)\n",
      "time taken: 24.8622949123\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEeRJREFUeJzt3X2MXFd9xvHv45g4b8VygdgQh5cqFJKqlUFgWqUSWwEB\nqoKjqg1vUoGUCjWiQIvUxKiVTVUBQQoSVZU/yptMFJS6SJAgaOOEsFRUIqEQQ4jd1Kh1CFa8SUkC\npJDUIb/+MXfbidn1zs7M7sz6fD/SaO+cOffe3x5bz9w5c+/dVBWSpJPfukkXIElaHQa+JDXCwJek\nRhj4ktQIA1+SGmHgS1Ijlgz8JBuS3Jrk9iR3Jnl/174pyb4kdyW5McnGvnV2JjmU5GCSi1byF5Ak\nDSaDnIef5Iyq+kmSU4B/Ad4DvBb4QVV9KMnlwKaquiLJBcC1wIuBrcDNwHPLE/4laaIGmtKpqp90\nixu6dR4EdgB7uvY9wMXd8muB66rqsao6DBwCto+rYEnScAYK/CTrktwOHAVmq+oAsLmq5gCq6ihw\ndtf9HOCevtWPdG2SpAlaP0inqnoceEGSJwM3JpkBjp+iccpGkqbYQIE/r6p+lOSLwIuAuSSbq2ou\nyRbgvq7bEeDcvtW2dm1PkMQ3CEkaQlVlmPWW/NI2yVOBY1X1wySnAzcC7wMuAh6oqisX+dL2JfSm\ncm5igS9tk9SwRU+DJLuravek6xiW9U/WWq5/LdcOJ0X9Q2fnIEf4Twf2JAm9Of9rqupL3Zz+3iSX\nAncDlwBU1YEke4EDwDHgMs/QkaTJWzLwq+oO4IULtD8AvHyRdT4AfGDk6iRJY+OVtsObnXQBI5qd\ndAEjmp10ASOanXQBI5iddAEjmp10AZMy0IVXK7LjNT6HL0mTMEp2eoQvSY0w8CWpEQa+JDXCwJek\nRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqE\ngS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiPWdOAnpx9NUsM9Tj866folaTWlqiaz46Sq\nKqNuA4atP4y6f0labaNk55o+wpckDW7JwE+yNcktSe5MckeSP+nadyX5fpJvdo9X9a2zM8mhJAeT\nXLSSv4AkaTBLTukk2QJsqar9Sc4CvgHsAF4H/LiqPnxc//OBTwMvBrYCNwPPreN25JSOJC3fik7p\nVNXRqtrfLT8MHATOmd/3AqvsAK6rqseq6jBwCNg+THGSpPFZ1hx+kmcD24Bbu6Z3JNmf5GNJNnZt\n5wD39K12hP9/g5AkTcjAgd9N53wGeFd3pH818EtVtQ04Cly1MiVKksZh/SCdkqynF/bXVNX1AFV1\nf1+XjwKf75aPAOf2vba1a1tou7v7ns5W1exAVUtSI5LMADNj2dYg5+En+RTwX1X1Z31tW6rqaLf8\np8CLq+qNSS4ArgVeQm8q5yb80laSxmKU7FzyCD/JhcCbgDuS3E4vYd8LvDHJNuBx4DDwdoCqOpBk\nL3AAOAZcdnzYS5JWn1faStIa4pW2kqQlGfiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXC\nwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8\nSWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqxJKBn2RrkluS3JnkjiTv7No3\nJdmX5K4kNybZ2LfOziSHkhxMctFK/gKSpMGkqk7cIdkCbKmq/UnOAr4B7ADeCvygqj6U5HJgU1Vd\nkeQC4FrgxcBW4GbguXXcjpJUVWWk4pOCE9d/grUZdf+StNpGyc4lj/Cr6mhV7e+WHwYO0gvyHcCe\nrtse4OJu+bXAdVX1WFUdBg4B24cpTpI0Psuaw0/ybGAb8DVgc1XNQe9NATi763YOcE/fake6NknS\nBK0ftGM3nfMZ4F1V9XBvOuUJlj23kmR339PZqppd7jYk6WSWZAaYGce2Bgr8JOvphf01VXV91zyX\nZHNVzXXz/Pd17UeAc/tW39q1/Zyq2j1U1ZLUiO5AeHb+eZJdw25r0CmdTwAHquojfW03AG/plt8M\nXN/X/vokpyZ5DnAecNuwBUqSxmOQs3QuBP4ZuIPetE0B76UX4nvpHc3fDVxSVQ916+wE/hA4Rm8K\naN8C2/UsHUlaplGyc8nAXykGviQt34qelilJOjkY+JLUiIFPyzz5bGCBU0sHdNpc1U+3jLceSVpZ\nTc/hO/8vaa1xDl+StCQDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHg\nS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4k\nNWL9pAtIcgq+8UjSilsyaJN8PMlckm/3te1K8v0k3+wer+p7bWeSQ0kOJrlo6RJ+4Vuw7lFY98jy\nHnl02F9aklo0yJH1J4FXLtD+4ap6Yff4J4Ak5wOXAOcDrwauTpITb/5nW+GewM/WLe/x5SW2K0nq\nt2TgV9VXgQcXeGmhwN0BXFdVj1XVYeAQsH2kCiVJYzHK3Pk7kuxP8rEkG7u2c4B7+voc6dokSRM2\n7Je2VwN/VVWV5K+Bq4C3LXcjSXbD+tN6q78GmBmyHEk6OSWZYUzhOFTgV9X9fU8/Cny+Wz4CnNv3\n2taubbHt7E7OfDe8ZwM8Y5hSJOmkVlWzwOz88yS7ht3WoFM6oW/OPsmWvtd+F/hOt3wD8PokpyZ5\nDnAecNuwxUmSxmfJI/wkn6b3ceIpSb4H7AJ+K8k24HHgMPB2gKo6kGQvcAA4BlxWVbUypUuSliOT\nyuMk3VcAZz4EhzYuf0rnK/Teh4atPyOtW1WeFipp1c1n5zDreoWrJDXCwJekRhj4ktQIA1+SGmHg\nS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWrEsH/isHEbSDLEvZVP\nm6v66Zal+0nS+Bn4Q3mU4e6ln83jrkSSBuWUjiQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9J\njTDwJakRBr4kNcLAl6RGGPiS1AgDX5IasWTgJ/l4krkk3+5r25RkX5K7ktyYZGPfazuTHEpyMMlF\nK1W4JGl5BjnC/yTwyuPargBurqrnAbcAOwGSXABcApwPvBq4OknGV64kaVhLBn5VfRV48LjmHcCe\nbnkPcHG3/Frguqp6rKoOA4eA7eMpVZI0imHn8M+uqjmAqjoKnN21nwPc09fvSNcmSZqwcf0BlGH+\nGghJdsP60+Aq4DXAzJjKkaSTQ5IZxhSOwwb+XJLNVTWXZAtwX9d+BDi3r9/Wrm1BVbU7OfPd8J4N\n8IwhS5Gkk1dVzQKz88+T7Bp2W4NO6aR7zLsBeEu3/Gbg+r721yc5NclzgPOA24YtTpI0Pkse4Sf5\nNL2PE09J8j1gF/BB4B+SXArcTe/MHKrqQJK9wAHgGHBZVQ013SNJGq9MKo+TVFUlOfMhOLRx+VM6\nX6H3PjRs/ZnAuqGqPE1V0tDms3OYdb3SVpIaYeBLUiMMfElqhIEvSY0Y14VXGsgGkgz5TfFpc1U/\n3TLeeiS1xMBfVY8ywplBm8dZiaT2OKUjSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJ\naoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcK/abtm\n+AfQJY3GwF8z/APokkbjlI4kNWKkI/wkh4EfAo8Dx6pqe5JNwN8DzwIOA5dU1Q9HrFOSNKJRj/Af\nB2aq6gVVtb1ruwK4uaqeB9wC7BxxH5KkMRg18LPANnYAe7rlPcDFI+5DkjQGowZ+ATcl+XqSt3Vt\nm6tqDqCqjgJnj7gPSdIYjHqWzoVVdW+SpwH7ktzFz59KMuypJZKkMRop8Kvq3u7n/Uk+B2wH5pJs\nrqq5JFuA+xZbP8luWH8aXAW8BpgZpRwtathz+D1/X5q0JDOMKRxTNdwBeJIzgHVV9XCSM4F9wPuA\nlwEPVNWVSS4HNlXVFQusX1WV5MyH4NBGeMYyK/gKvTEY+tz0Caw7iX2Osm6oqgy5U0krYD47h1l3\nlCP8zcBnuyPH9cC1VbUvyb8Ce5NcCtwNXDLCPiRJYzJ04FfVfwLbFmh/AHj5KEVJksbPK20lqREG\nviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1IhR/wCKTmrD\n3kcfvJe+NH0MfJ3Ao4xwD/7N46xE0uic0pGkRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREG\nviQ1wsCXpEZ4pa1WyLC3ZfCWDNJKMfC1Qoa9LYO3ZJBWilM6Omkkpx9NUsM9Tj866fqlleYRvk4i\nj2z2Zm/S4jzCl6RGrFjgJ3lVkn9L8u9JLl+p/ehk0/uyd5jHpCuXpt2KBH6SdcDfAq8EfgV4Q5Ln\nr8S+Jmd20gWMaHbSBSxi/svepR5fXqBt7UgyM+kahrWWa4e1X/8oVuoIfztwqKrurqpjwHXAjhXa\n14TMTrqAEc1OuoARzY55e8N+shj6y96ZcVa/ymYmXcCIZiZdwKSs1Je25wD39D3/Pr03AWlKDXsa\n6Wmbh51OSk7/S3hkyIOu0x4fbl2vc2jZNJyl8z/w+z+GDY8vb7UH1wNnrkhF0sCGfaPYDbxv3Qhn\nFQ257vBnI/U+zTyyubecXctbe7g3mv59Lt/ib4pL139yvjGmavxzn0l+HdhdVa/qnl8BVFVd2ddn\nbU26StKUqKoMs95KBf4pwF3Ay4B7gduAN1TVwbHvTJI0kBWZ0qmqnyV5B7CP3hfDHzfsJWmyVuQI\nX5I0fVbtStskm5LsS3JXkhuTbFyk3+Ek30pye5LbVqu+xQxyAVmSv0lyKMn+JNtWu8YTWar+JC9N\n8lCSb3aPv5hEnQtJ8vEkc0m+fYI+0zz2J6x/ysd+a5JbktyZ5I4k71yk31SO/yD1T/n4b0hya5eD\ndyZ5/yL9ljf+VbUqD+BK4M+75cuBDy7S7z+ATatV1xI1rwO+CzwLeBKwH3j+cX1eDXyhW34J8LVJ\n173M+l8K3DDpWhep/zeBbcC3F3l9asd+wPqneey3ANu65bPofSe3lv7vD1L/1I5/V98Z3c9TgK8B\nF446/qt5L50dwJ5ueQ9w8SL9wvTc42eQC8h2AJ8CqKpbgY3J1NyIa9AL4Ib6xn+lVdVXgQdP0GWa\nx36Q+mF6x/5oVe3vlh8GDtK7vqbf1I7/gPXDlI4/QFX9pFvcQC8Tj/+/tOzxX81gPbuq5qD3jwGc\nvUi/Am5K8vUkf7Rq1S1soQvIjv9Pc3yfIwv0mZRB6gf4je4j4ReSXLA6pY3FNI/9oKZ+7JM8m94n\nlVuPe2lNjP8J6ocpHv8k65LcDhwFZqvqwHFdlj3+Yz1LJ8lNQP87TOgF+EJzY4t9W3xhVd2b5Gn0\ngv9gd6SklfEN4JlV9ZMkrwY+B/zyhGtqxdSPfZKzgM8A7+qOlNeUJeqf6vGvqseBFyR5MrAvyUur\n6iujbHOsR/hV9Yqq+rW+x692P28A5uY/biTZAty3yDbu7X7eD3yWyd6S4QjwzL7nW7u24/ucu0Sf\nSVmy/qp6eP6jY1X9I/CkJL+4eiWOZJrHfknTPvZJ1tMLy2uq6voFukz1+C9V/7SP/7yq+hHwBeBF\nx7207PFfzSmdG4C3dMtvBn7uHyDJGd07MknOBC4CvrNaBS7g68B5SZ6V5FTg9fR+j343AH8A/3eF\n8UPzU1dTYMn6++f8kmynd6ruA6tb5gmFxedZp3ns5y1a/xoY+08AB6rqI4u8Pu3jf8L6p3n8kzx1\n/kzGJKcDr6B30kW/ZY//at5L50pgb5JLgbuBSwCSPB34aFX9Dr3poM+md9uF9cC1VbVvFWt8glrk\nArIkb++9XH9XVV9M8ttJvgv8N/DWSdV7vEHqB34vyR8Dx4CfAq+bXMVPlOTT9O5s+JQk3wN2Aaey\nBsYelq6f6R77C4E3AXd088gFvJfeGV9TP/6D1M8Ujz/wdGBPkvmTWK6pqi+Nmj1eeCVJjZiW0x8l\nSSvMwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqRH/C4Y9KKv6/GRaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe2e1bcea50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plain performance test\n",
    "\n",
    "X,Y = sample_gaussian_vs_laplace()\n",
    "print X.shape, Y.shape\n",
    "\n",
    "feats_p = RealFeatures(X.reshape(1,len(X)))\n",
    "feats_q = RealFeatures(Y.reshape(1,len(Y)))\n",
    "width=1\n",
    "k = GaussianKernel(10, width) # should have a constructor that does not require to specify cache size\n",
    "\n",
    "mmd = QuadraticTimeMMD() # should take p and q as argument\n",
    "mmd.set_p(feats_p)\n",
    "mmd.set_q(feats_q)\n",
    "mmd.set_kernel(k)\n",
    "\n",
    "# this should be able to receive a CSignal abort, made my Python stall a couple of times\n",
    "#mmd.set_num_null_samples(10000)\n",
    "\n",
    "start = time.time()\n",
    "mmd.set_num_null_samples(1000) \n",
    "null_samples = mmd.sample_null() # seems to take quite long, and not running on all my cpus (only 2/4 it seems). I havent checked properly, but this is only first impression\n",
    "print \"time taken:\", time.time()-start\n",
    "plt.hist(null_samples, bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000,) (1000,)\n",
      "time taken: 0.209638118744\n",
      "time taken: 36.7916610241\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAECZJREFUeJzt3X+sZGV9x/H3B1dgFaRIy24rspY01NRK0YhYtXqDsaCt\nRZPW+qOxamIb26qxjQWpCftXIzaNtT/8g3bZYOvGqFUKVSMSekmwRa38WMQVMbZAKXutCv6oLEr5\n9o+ZlZvrvXtnzpyZufe571cyYc6Z85zz5eTsZ5557pxnUlVIktpw1LwLkCT1x1CXpIYY6pLUEENd\nkhpiqEtSQwx1SWrIuqGeZE+SpST7l607K8lnk9w0/O8zplumJGkUo/TU9wLnrlj3LuAdVfU04GLg\nz/ouTJI0vnVDvaquB+5bsfpe4ITh8x8D7um5LklSBxnljtIku4CrquqM4fKpwKeBAgI8u6runmah\nkqT1df1D6R7gTVV1KvBW4LL+SpIkddW1p/7tqnrcste/VVUnrNHWyWUkqYOqyrhtto24XYaPw+5I\n8vyqui7JC4Av911YV4M3kXHfRzKTGpPsrqrd0z7OZuC5eITn4hGei0d07RCvG+pJ9gELwElJ7mLw\nbZffAd6b5Gjg0HBZkjRn64Z6Vb1qjZfO7rkWSdKEvKN0thbnXcAGsjjvAjaQxXkXsIEszruAzW6k\nP5ROdICkHFOXpPF0zU576pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBD\nXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIeuGepI9SZaS7F+x/k1JDiS5Nck7p1eiJGlUo/xG6V7gr4D3\nHV6RZAF4CfDUqnooyY9PpzxJ0jjW7alX1fXAfStWvxF4Z1U9NNzm61OoTZI0pq5j6qcDz0tyQ5J/\nSfKMPouSJHUzyvDLWu1OrKpnJTkL+CBw2lobJ9m9bHGxqhY7HleSmjQc1l6YeD+j/EZpkl3AVVV1\nxnD548AlVXXdcPkrwNlV9Y1V2vobpZI0pmn/RmmGj8OuAM4ZHvh04NGrBbokabbWHX5Jso/BR4KT\nktwFXAxcBuxNcivwIPCaaRYpSRrNSMMvEx3A4RdJGtu0h18kSZuAoS5JDTHUJakhhrokNcRQl6SG\nGOqS1BBDXZIaYqhLUkMM9Qkk2w8mqfEf2w/Ou3ZJbfKO0kGrTneUdjtW9+NJ2jq8o1SSZKhLUksM\ndUlqiKEuSQ0x1CWpIYa6JDVk3VBPsifJUpL9q7z2R0keTvL46ZQnSRrHKD31vcC5K1cmOQV4IXBn\n30VJkrpZN9Sr6nrgvlVeejfwtt4rkiR11mlMPcmvAXdX1a091yNJmsC2cRsk2Q5cxGDo5Yer12mz\ne9niYlUtjntcSWpZkgVgYeL9jDL3S5JdwFVVdUaSnweuAb7HIMxPAe4BnllVX1ulrXO/9HQ8SVtH\n1+wctaee4YOq+gKwc9mB/wN4elWtNu4uSZqhUb7SuA/4V+D0JHcled2KTYp1hl8kSbPh1LuDVg6/\nSNpQnHpXkmSoS1JLDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUkLEn9GrTMcMbiSRpczPUAXiQ\nrneGStJG4vCLJDXEUJekhhjqktQQQ12SGmKoS1JDDHVJasgov3y0J8lSkv3L1r0ryYEkNyf5xySP\nm26ZkqRRjNJT3wucu2Ld1cBTqupM4A7g7X0XJkka37qhXlXXA/etWHdNVT08XLwBOGUKtUmSxtTH\nmPrrgU/0sB9J0oQmmiYgyZ8AP6iqfetst3vZ4mJVLU5yXElqTZIFYGHi/VStP+dJkl3AVVV1xrJ1\nrwXeAJxTVQ8eoW2nX8TuajAx17jzuITuc790azfLcyJp8+manaP21MOy2auSnAe8DXjekQJdkjRb\n6/bUk+xj8JHgJGAJuBi4CDga+MZwsxuq6vfWaG9PfZV29tQlHUnX7Bxp+GUShvrq7Qx1SUfSNTu9\no1SSGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHU\nJakhhrokNcRQl6SGrBvqSfYkWUqyf9m6E5NcneT2JJ9McsJ0y5QkjWKUnvpe4NwV6y4ErqmqnwWu\nBd7ed2GSpPGtG+pVdT1w34rV5wOXD59fDry057okSR10HVM/uaqWAKrqIHByfyVJkrrq6w+l0/2h\nU0nSSLZ1bLeUZEdVLSXZCXztSBsn2b1scbGqFjsed0tLth+EQzvGa3XsUtUDO6dTkaS+JFkAFibe\nT9X6newkTwKuqqqnDpcvAb5ZVZckuQA4saouXKNtp1/E7ipJjf/BIXT7sNG9XZdz0vX/bZbnX1I/\numbnuqGeZB+Dd4+TgCXgYuAK4EPAE4E7gZdX1f19FtaVod7PsSTN19RCfVKG+urtDHVJR9I1O72j\nVJIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQl\nqSGGuiQ1xFCXpIYY6pLUkIlCPcnbk9yWZH+S9yc5uq/CJEnj6xzqSXYBbwCeVlVnMPgR61f0VZgk\naXzbJmj7beD7wGOTPAw8BvjvXqqSJHXSuadeVfcBfw7cBdwD3F9V1/RVmCRpfJMMv5wGvBXYBfwU\ncFySV/VVmCRpfJMMvzwD+HRVfRMgyUeAZwP7Vm6YZPeyxcWqWlxv58n2g3BoxwT1SdKmkWQBWJh4\nP1XVtYBfAP4BOAt4ENgLfK6q/mbFdlVV6bD/gi61hfHbdWkzWbvZnZNux5I0X12zc5Ix9VuA9wGf\nB25hkHCXdt2fJGlynXvqIx/Anvqq7eypSzqSmffUJUkbj6EuSQ0x1CWpIYa6JDXEUJekhhjqktQQ\nQ12SGjLJNAEjWzFNwCgenkYdktS6mdx8BBeP2eoD34fbj2735qNjGcys0IU3H0lbQecbN2cT6uMe\n42XfgSuObzfUZ9nOUJc2I+8olSQZ6pLUEkNdkhpiqEtSQwx1SWqIoS5JDZko1JOckORDSQ4kuS3J\n2X0VJkka36R3lL4H+HhV/UaSbcBjeqhJktRR51BP8jjgl6rqtQBV9RDw7Z7qkiR1MMnwy08DX0+y\nN8mNSS5Nsr2vwiRJ45tk+GUb8HTg96vq35P8BXAhq070snvZ84XhQ7NxzHCqhnEd+zAc6vCmf+xS\n1QM7x28nbW1JFughHDvP/ZJkB/BvVXXacPm5wAVV9ZIV2zn3y1zbzb5G55qRJjfzuV+qagm4O8np\nw1UvAL7YdX+SpMlN+u2XNwPvT/Jo4KvA6yYvSZLU1UShXlW3AGf1VIskaULeUSpJDTHUJakhhrok\nNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JD\nDHVJasjEoZ7kqCQ3Jrmyj4IkSd310VN/C/42qSRtCBOFepJTgBcDf9dPOZKkSUzaU3838DageqhF\nkjShzj88neRXgKWqujnJApC1t9697PnC8CFJOmyYowsT76eqWyc7yZ8CvwU8BGwHjgc+UlWvWbFd\njd+Rf9l34Irju30ACOO369Jms7SbfY1VdYQ3eEmjSFJd/i11Hn6pqouq6tSqOg14BXDtykCXJM2W\n31OXpIZ0HlNfrqquA67rY1+SpO7sqUtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIa\nYqhr00q2H0xS4z+2H5x37dK09HJHqTQfh3Z0nHRsR++lSBuEPXVJaoihLkkNMdQlqSGGuiQ1xFCX\npIYY6pLUkM6hnuSUJNcmuS3JrUne3GdhkqTxTfI99YeAP6yqm5McB3w+ydVV9aWeapMkjWmSH54+\nWFU3D59/FzgAPKGvwiRJ4+tlTD3Jk4Azgc/0sT9JUjcTTxMwHHr5MPCWYY99FbuXPV8YPqRHDOZj\nOeTt+9qykizQQzimqsvcGT8sYhvwz8Anquo9a2xT48/P8bLvwBXHd5zXg/HbdWmzWdrNvsaqytit\nOl0ns61RmqUk1eU6nXT45TLgi2sFuiRptib5SuNzgFcD5yS5KcmNSc7rrzRJ0rg6j6lX1aeBR/VY\niyRpQt5RKkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQyae+0XafI4ZTkswrmMfhkMdOkJd\n2h27VPXAzvGP1U33uXfGr3OWx9qKDHVtQQ/Scc6Yo2bXLjOe3OzQjo7/bx3qnOWxth6HXySpIYa6\nJDXEUJekhhjqktQQQ12SGmKoS1JDJgr1JOcl+VKSLye5oK+iJEndTPLLR0cBfw2cCzwFeGWSJ/dV\nWJsW513AhjH8kV0BXhfLLc67gE1vkp76M4E7qurOqvoB8AHg/H7KatXivAvYSBbmXcDGsTjvAjaQ\nxXkXsOlNEupPAO5etvxfw3WSpDmZ0TQB53xrvO1vOWY6dUhS21LVZQ4GSPIsYHdVnTdcvhCoqrpk\nxXbdDiBJW1xVZdw2k4T6o4DbgRcA9wKfBV5ZVQc67VCSNLHOwy9V9X9J/gC4msHY/B4DXZLmq3NP\nXZK08fR+R2mSE5NcneT2JJ9McsIa2/1nkluS3JTks33XMU+j3JSV5C+T3JHk5iRnzrrGWVjvPCR5\nfpL7k9w4fLxjHnXOQpI9SZaS7D/CNs1fE7D+udhi18UpSa5NcluSW5O8eY3tRr82qqrXB3AJ8MfD\n5xcA71xju68CJ/Z9/Hk/GLxRfgXYBTwauBl48optXgR8bPj8bOCGedc9p/PwfODKedc6o/PxXOBM\nYP8arzd/TYxxLrbSdbETOHP4/DgGf6ecKC+mMffL+cDlw+eXAy9dY7vQ5twzo9yUdT7wPoCq+gxw\nQtLcr7qMenPa2H/d34yq6nrgviNsshWuCWCkcwFb57o4WFU3D59/FzjAj97vM9a1MY1QPbmqlg4X\nDJy8xnYFfCrJ55K8YQp1zMsoN2Wt3OaeVbbZ7Ea9Oe0Xhx8pP5bk52ZT2oa0Fa6JcWy56yLJkxh8\ngvnMipfGujY6ffslyaeA5e8UYRDSq419rfWX2OdU1b1JfoJBuB8YvoNr6/g8cGpVfS/Ji4ArgNPn\nXJPmb8tdF0mOAz4MvGXYY++sU6hX1QvXem34B5AdVbWUZCfwtTX2ce/wv/+T5KMMPq63EOr3AKcu\nWz5luG7lNk9cZ5vNbt3zsPzirapPJHlvksdX1TdnVONGshWuiZFstesiyTYGgf73VfVPq2wy1rUx\njeGXK4HXDp//NvAjRSZ5zPCdiSSPBX4Z+MIUapmHzwE/k2RXkqOBVzA4J8tdCbwGfnhn7v2Hh6wa\nsu55WD4umOSZDL5i2+Q/3KGw9ljxVrgmllvzXGzB6+Iy4ItV9Z41Xh/r2pjG3C+XAB9M8nrgTuDl\nw2J+EvjbqvpVBkM3Hx1OIbANeH9VXT2FWmau1rgpK8nvDl6uS6vq40lenOQrwP8Cr5tnzdMwynkA\nfj3JG4EfAA8Avzm/iqcryT4GM1OelOQu4GLgaLbQNXHYeueCrXVdPAd4NXBrkpsYDFdfxOBbY52u\nDW8+kqSGtPiVQknasgx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIa8v/NJlQ8uxRW2AAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe2e1aae150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# same as above but with more data\n",
    "\n",
    "X,Y = sample_gaussian_vs_laplace(n=1000)\n",
    "print X.shape, Y.shape\n",
    "\n",
    "feats_p = RealFeatures(X.reshape(1,len(X)))\n",
    "feats_q = RealFeatures(Y.reshape(1,len(Y)))\n",
    "width=1\n",
    "k = GaussianKernel(10, width)\n",
    "\n",
    "joint_features = RealFeatures(np.hstack((X, Y)).reshape(1, len(X)+len(Y)))\n",
    "k.init(joint_features,joint_features)\n",
    "start = time.time()\n",
    "k.get_kernel_matrix()\n",
    "print \"time taken:\", time.time()-start\n",
    "\n",
    "mmd = QuadraticTimeMMD()\n",
    "mmd.set_p(feats_p)\n",
    "mmd.set_q(feats_q)\n",
    "mmd.set_kernel(k)\n",
    "\n",
    "start = time.time()\n",
    "mmd.set_num_null_samples(100) # Is the kernel matrix precomputed in here? It should be. Why isnt this faster?\n",
    "# when doing permutation test, is the kernel matrix only computed once and then summed over differently? Or is it copied around?\n",
    "# what takes the time here?\n",
    "null_samples = mmd.sample_null()\n",
    "print \"time taken:\", time.time()-start\n",
    "plt.hist(null_samples, bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5000,) (5000,)\n",
      "time taken: 12.797287941\n",
      "time taken: 157.862884998\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADidJREFUeJzt3X2sZHV9x/H3Z11R6LZUsawNK9jiUzQaqgYxmLC1sQJN\nwD9MfEptSdoSU6KpacUYE9c0aeo/bWq1pbTUiI3RRFPcFqzY2FtLG5ACC1RBoQ+Kq2x9oo0CFvXb\nP+Zgxsu9O+feO3fm7pf3KznZc2Z+c84nZ9nPnPubOZdUFZKknnYtO4AkaftY8pLUmCUvSY1Z8pLU\nmCUvSY1Z8pLU2MyST7IvyaeSfDbJ7UneuM64dye5K8mhJGfMP6okaaN2jxjzPeDNVXUoyR7gpiTX\nVtWdDw9Ich5welU9PcmLgMuAs7YnsiRprJlX8lV1b1UdGta/DdwBnLJq2IXAlcOYG4ATk+ydc1ZJ\n0gZtaE4+yVOBM4AbVj11CnDP1PZhHvlGIElasNElP0zVfAR403BFL0na4cbMyZNkN5OC/0BVfWyN\nIYeBp0xt7xseW70ff1GOJG1CVWWzL5y5MJlv/4OjPH8+cPWwfhZw/Trjaszxlr0AB5adwZzmPFYz\nmnPN4xTUJhdqK90580o+ydnA64Dbk9wyHPBtwGnDgS+vqmuSnJ/kbuA7wEWbeseRJM3VzJKvqn8G\nHjNi3CVzSSRJmhvveF3byrIDjLSy7AAjrSw7wEgryw4wwsqyA4y0suwAI60sO8B2yzBftJiDJVWb\n/fBAko5Rky+dbLZrJ5W52e70Sl6SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16S\nGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPk\nJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakx\nS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGrPkJakxS16SGptZ8kmuSHIk\nyW3rPH9OkvuS3Dwsb59/TEnSZuweMeZ9wB8DVx5lzKer6oL5RJIkzcvMK/mqug741oxhmU8cSdI8\nzWtO/sVJDiW5Osmz57RPSdIWjZmumeUm4NSquj/JecBVwDPWG5zkwNTmSlWtzCGDJDWyMixbl6qa\nPSg5DfibqnreiLH/Cbygqr65xnNVVU7tSHpUSVIwu2vXeTUAm+3OsdM1YZ159yR7p9bPZPLG8YiC\nlyQt3szpmiQfBPYDJyX5EvAO4Digqupy4JVJ3gA8BDwAvGr74kqSNmLUdM3cDuZ0jaRHoWNhukaS\ndAyy5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz\n5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWp\nMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUte\nkhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqbWfJJrkhyJMltRxnz7iR3JTmU\n5Iz5RpQkbdaYK/n3AS9f78kk5wGnV9XTgYuBy+aUTZK0RTNLvqquA751lCEXAlcOY28ATkyydz7x\nJElbsXsO+zgFuGdq+/Dw2JE57HtDkjyWyU8d2eQuvlZV188xkiQt1TxKfkOSHJjaXKmqlTnu/nVw\n8p/Cc767uZd/ek9y/DfgwZM39/rH/wAe3OSH2Vt5LcDjj1Q98OTNvDI5/l54cAs/fR2bx96qrWVf\nXu5l8pyNtTIsWzePkj8MPGVqe9/w2Jqq6sAcjrmex8LLfwBXnri5l5/wXXjgZKhNHj67lvNagK1M\nkT2499F57K3aSvZH65Sm52yc/cMC8M4t7WnslWNYfwrkIPB6gCRnAfdV1cKnaiRJjzTzSj7JB5m8\npZyU5EvAO4DjgKqqy6vqmiTnJ7kb+A5w0XYGliSNN7Pkq+q1I8ZcMp84kqR58o5XSWrMkpekxix5\nSWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrM\nkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpek\nxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWrMkpekxix5\nSWrMkpekxix5SWrMkpekxix5SWrMkpekxix5SWpsVMknOTfJnUm+kOTSNZ4/J8l9SW4elrfPP6ok\naaN2zxqQZBfwHuAXgK8ANyb5WFXduWrop6vqgm3IKEnapDFX8mcCd1XVF6vqIeBDwIVrjMtck0mS\ntmxMyZ8C3DO1/eXhsdVenORQkquTPHsu6SRJWzJzumakm4BTq+r+JOcBVwHPWGtgkgNTmytVtTKn\nDJLUxMqwbN2Ykj8MnDq1vW947Ieq6ttT6x9P8idJnlhV31y9s6o6sMmskvQosX9YAN65pT2Nma65\nEXhaktOSHAe8Gjg4PSDJ3qn1M4GsVfCSpMWaeSVfVd9PcglwLZM3hSuq6o4kF0+ersuBVyZ5A/AQ\n8ADwqu0MLUkaZ9ScfFX9HfDMVY/92dT6e4H3zjeaJGmrvONVkhqz5CWpMUtekhqz5CWpMUtekhqz\n5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWp\nMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUte\nkhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz\n5CWpMUtekhqz5CWpsVEln+TcJHcm+UKSS9cZ8+4kdyU5lOSM+caUJG3GzJJPsgt4D/By4DnAa5I8\na9WY84DTq+rpwMXAZduQdYFWlh1gpJVlBxglyf5lZxjjWMh5LGQEc+4kY67kzwTuqqovVtVDwIeA\nC1eNuRC4EqCqbgBOTLJ3rkkXamXZAUZaWXaAsfYvO8BI+5cdYIT9yw4w0v5lBxhp/7IDbLcxJX8K\ncM/U9peHx4425vAaYyRJC7Z72QHm7CH4xC546f9s7uX/t2e+cSRpuVJVRx+QnAUcqKpzh+23AlVV\n75oacxnwD1X14WH7TuCcqjqyal9HP5gkaU1Vlc28bsyV/I3A05KcBnwVeDXwmlVjDgK/CXx4eFO4\nb3XBbyWkJGlzZpZ8VX0/ySXAtUzm8K+oqjuSXDx5ui6vqmuSnJ/kbuA7wEXbG1uSNMbM6RpJ0rFr\nW+94TfKEJNcm+XySTyQ5cZ1x/5Xk1iS3JPnMdmbaSs5h7K4kNyc5uMiMw7Fn5kzyuCQ3DOfys0l+\nb4fm3JfkU0PG25O8cSfmHMZdkeRIktsWmO2YuAFxVs4kz0zyL0keTPLmZWQccszK+dqhg25Ncl2S\n5+7QnBdMdeW/JnnpzJ1W1bYtwLuAtwzrlwK/v864/wCesJ1Z5pFzeP63gL8CDu7UnMAJw5+PAa4H\nzt5pOYEnA2cM63uAzwPP2mk5h+deApwB3LagXLuAu4HTgMcCh1afG+A84Oph/UXA9Ys8dxvI+STg\nBcDvAm9edMYN5DwLOHFYP3cHn88TptafC9w9a7/b/btrLgTeP6y/H3jFOuPCcn+PzqicSfYB5wN/\nsaBcq43KWVX3D6uPY3Jev7X90X7EzJxVdW9VHRrWvw3cweLvrRh7Pq9jsefwWLkBcWbOqvp6Vd0E\nfG/B2aaNyXl9VT381evrWc59PmNy3j+1uQf4+qydbnexnlzDt2yq6l7g5HXGFfDJJDcm+fVtzrSW\nsTn/EPgdJnmXYVTOYUrpFuBeYKWqPrfAjDD+fAKQ5KlMrpRv2PZkP2pDORfoWLkBcUzOnWCjOX8N\n+Pi2JlrbqJxJXpHkDuAaYOY055ZvhkrySWD6CiJMSvDtawxfrxzPrqqvJvkpJmV/x3D1NDdbzZnk\nl4AjVXVo+H0X2/J10Hmcz6r6AfBzSX4CuDbJOVX1jzst57CfPcBHgDcNV/RzNa+cenRI8vNMvh34\nkmVnWU9VXQVcleQlwAeAZx5t/JZLvqpett5zw4dVe6vqSJInA/+9zj6+Ovz5tSR/zeTHlrmW/Bxy\nng1ckOR84Hjgx5NcWVWv32E5p/f1v0muBl4IzLXk55EzyW4mBf+BqvrYPPPNM+cSHAZOndreNzy2\nesxTZozZbmNy7gSjciZ5HnA5cG5VLXqKEzZ4PqvquiS7k5xUVd9Yb9x2T9ccBH51WP8V4BH/kJOc\nMFzNkeTHgF8E/m2bc602M2dVva2qTq2qn2VyQ9in5l3wI4w5n096+FsiSY4HXsbkA5xFmplz8JfA\n56rqjxYRag1jc8LkJ4BF3cz3wxsQkxzH5L+31d/mOgi8Hn54V/qaNyBuszE5py3rZsiZOZOcCnwU\n+OWq+vclZIRxOU+fWn8+wNEKnmHAdn5a/ETg75l8c+Ja4CeHx38a+Nth/WeYlNAtwO3AWxf1afZG\ncq4afw7L+XbNmPP5XODm4XzeCvz2Ds15NvD9qb/7m5lcQe2onMP2B4GvAN8FvgRctIBs5w657nr4\n3wSTX+P9G1Nj3sPk2xi3As9f9N/zmJxMpsruAe4Dvjmcvz07MOefA9+Y+rfzmR16Pt/C5CL4ZuCf\ngBfO2qc3Q0lSY/7v/ySpMUtekhqz5CWpMUtekhqz5CWpMUtekhqz5CWpMUtekhr7fyZE1CSyV5fC\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe2e1aaee50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# same as above, how far can be go in terms of data?\n",
    "\n",
    "X,Y = sample_gaussian_vs_laplace(n=5000)\n",
    "print X.shape, Y.shape\n",
    "\n",
    "feats_p = RealFeatures(X.reshape(1,len(X)))\n",
    "feats_q = RealFeatures(X.reshape(1,len(Y)))\n",
    "width=1\n",
    "k = GaussianKernel(10, width)\n",
    "\n",
    "joint_features = RealFeatures(np.hstack((X, Y)).reshape(1, len(X)+len(Y)))\n",
    "k.init(joint_features,joint_features)\n",
    "start = time.time()\n",
    "k.get_kernel_matrix()\n",
    "print \"time taken:\", time.time()-start\n",
    "\n",
    "mmd = QuadraticTimeMMD()\n",
    "mmd.set_p(feats_p)\n",
    "mmd.set_q(feats_q)\n",
    "mmd.set_kernel(k)\n",
    "\n",
    "start = time.time()\n",
    "mmd.set_num_null_samples(10)\n",
    "null_samples = mmd.sample_null()\n",
    "print \"time taken:\", time.time()-start\n",
    "plt.hist(null_samples, bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100,)\n",
      "0 0.000113010406494\n",
      "1 3.46533489227\n",
      "2 2.56397795677\n",
      "3 3.12962388992\n",
      "4 2.77851605415\n",
      "5 2.36319494247\n",
      "6 2.64674282074\n",
      "7 2.62075686455\n",
      "8 2.64666509628\n",
      "9 2.66596198082\n",
      "10 2.69098114967\n",
      "11 2.74386906624\n",
      "12 9.1876540184\n",
      "13 7.1437599659\n",
      "14 3.66940402985\n",
      "15 2.86414289474\n",
      "16 2.65622401237\n",
      "17 2.89886879921\n",
      "18 2.60687017441\n",
      "19 2.91620898247\n"
     ]
    }
   ],
   "source": [
    "# are p-values uniformly distributed if H0 is true?\n",
    "# Initial try, too slow. Made smaller below, but this case also is interesting\n",
    "\n",
    "X,_ = sample_gaussian_vs_laplace(n=100)\n",
    "X2,_ = sample_gaussian_vs_laplace(n=100)\n",
    "\n",
    "print X.shape\n",
    "\n",
    "feats_p = RealFeatures(X.reshape(1,len(X)))\n",
    "feats_q = RealFeatures(X2.reshape(1,len(X2)))\n",
    "width=1\n",
    "k = GaussianKernel(10, width)\n",
    "\n",
    "mmd = QuadraticTimeMMD()\n",
    "mmd.set_p(feats_p)\n",
    "mmd.set_q(feats_q)\n",
    "mmd.set_kernel(k)\n",
    "\n",
    "\n",
    "mmd.set_num_null_samples(500) # if this is not called, there is no default and p-values are all nan\n",
    "\n",
    "#mmd.set_null_approximation_method()\n",
    "\n",
    "num_runs = 20\n",
    "p_values = np.zeros(num_runs)\n",
    "\n",
    "last = time.time()\n",
    "for i in range(num_runs):\n",
    "    print i, time.time()-last\n",
    "    last = time.time()\n",
    "    \n",
    "    # this is a usual use-case: repeating the same experiment multiple times.\n",
    "    # I guess we can save computation here, e.g. the precomputed kernel matrix\n",
    "    # .... maybe use some kind of lazy evaluation?\n",
    "    # also, notice the highly inhomogeneous times. Why is that?\n",
    "    p_values[i] = mmd.compute_p_value(0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADTVJREFUeJzt3W+sZPVdx/H3B5bypyQbxLIbS1uCD9rG0GwaS5pQ41aT\ngk+EYEItJmoTCU9oGx+BTQzrwxpDbDSY1FJCTJtGMdBCpKWKN7pNkNp2LaUrUhUUs7tFpI0o4Cpf\nH9whnV3u7p07c2bO3e99v5LJnXPmnPl97y+/+5lzf3POTKoKSVJPZ41dgCRpeQx5SWrMkJekxgx5\nSWrMkJekxgx5SWps05BPcmmSR5I8keTxJB+ZrL89ybNJvjG5XbP8ciVJW5HNzpNPshfYW1WHklwI\nfB24Fvgg8J9Vdcfyy5QkzWPXZhtU1VHg6OT+i0kOA2+ePJwl1iZJWtCW5uSTXAbsA/5msuqWJIeS\nfDrJ7oFrkyQtaOaQn0zV3At8rKpeBO4ELq+qfawf6TttI0nbzKZz8gBJdgEPAg9V1Sc3ePxtwANV\n9a4NHvPDcSRpDlW18JT4pnPyE58BvjMd8En2TubrAa4Hvn2qnYcoVOuSHKiqAyO1XTDWa3YGH0dj\n9mVH9uewhjpA3jTkk1wF/BLweJJvsv5X/nHgxiT7gFeBp4GbhyhIkjScWc6u+Spw9gYPfWn4ciRJ\nQ/KK1zPP2tgFNLI2dgHNrI1dgF5vpjdeF2ogKefke+g2Jy9tZ0Nlp0fyktSYIS9JjRnyktSYIS9J\njRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRny\nktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSY\nIS9JjRnyktSYIS9JjRnyktSYIS9JjW0a8kkuTfJIkieSPJ7ko5P1FyV5OMmTSb6cZPfyy5UkbUWq\n6vQbJHuBvVV1KMmFwNeBa4EPA89X1W8nuRW4qKpu22D/qqosoXatWJKC04+XJbaO40g7yVDZuemR\nfFUdrapDk/svAoeBS1kP+nsmm90DXLdoMZKkYW1pTj7JZcA+4FFgT1Udg/UXAuCSoYuTJC1m16wb\nTqZq7gU+VlUvrv/rfoJT/h+f5MDU4lpVrW2lSEnqLsl+YP/gz7vZnPyk8V3Ag8BDVfXJybrDwP6q\nOjaZt//LqnrnBvs6J9+Ec/LS6qxsTn7iM8B3Xgv4iS8Cvzq5/yvAFxYtRpI0rFnOrrkK+CvgcdYP\n4wr4OPAY8MfAW4BngBuq6vsb7O+RfBMeyUurM1R2zjRds1ADhnwbhry0OquerpEknYEMeUlqzJCX\npMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYM\neUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlqzJCXpMYMeUlq\nzJCXpMYMeUlqzJCXpMZ2jV2Ati45/yi8vGfsOlbrXJLU6ts971jVS3tX3640jFQt9+8mSVVVltrI\nDrMediPkHWGcdsdsOzh+NYahstPpGklqzJCXpMY2DfkkdyU5luRbU+tuT/Jskm9Mbtcst0xJ0jxm\nOZK/G7h6g/V3VNW7J7cvDVyXJGkAm4Z8VR0EXtjgId+MkqRtbpE5+VuSHEry6SS7B6tIkjSYeUP+\nTuDyqtoHHAXuGK4kSdJQ5roYqqqem1r8Q+CB022f5MDU4lpVrc3TriR1lWQ/sH/w553lYqgklwEP\nVNUVk+W9VXV0cv/XgfdU1Y2n2NeLoQbmxVCrbdfxqzEMlZ2bHskn+Rzrry4XJ/kX4Hbg/Un2Aa8C\nTwM3L1qIJGl4fqzBGcgj+dW26/jVGPxYA0nSpgx5SWrMkJekxgx5SWrMkJekxvxmKOm0/EYqndkM\neem0XmGkUzd32Nc7almcrpGkxgx5SWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5\nSWrMkJekxgx5SWrMkJekxgx5SWrMkJekxgx5SWrMLw1ZQHL+UXjZL3eQtG0Z8gt5ec9I3xo0QpuS\nzkRO10hSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY5uGfJK7\nkhxL8q2pdRcleTjJk0m+nGT3csuUJM1jliP5u4GrT1p3G/DnVfV24BHgN4YuTJK0uE1DvqoOAi+c\ntPpa4J7J/XuA6wauS5I0gHnn5C+pqmMAVXUUuGS4kiRJQxnqjdcxPlRdkrSJeb805FiSPVV1LMle\n4Hun2zjJganFtapam7NdSWopyX5g/+DPW7X5QXiSy4AHquqKyfIngP+oqk8kuRW4qKpuO8W+VVUt\nv8ooSY33zVA7qd0x2x6v3a5/N5rNUNm5acgn+Rzrry4XA8eA24H7gT8B3gI8A9xQVd9fZqHbkSG/\nE9o25DWOlYX8wg0Y8stomZ3V7phtG/Iax1DZ6RWvktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRny\nktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSYIS9JjRnyktSY\nIS9JjRnyktSYIS9Jje0auwBJGzmXJLX6ds87VvXS3tW3q2Ux5KVt6RVghIwne0ZoVEvkdI0kNWbI\nS1JjhrwkNWbIS1JjhrwkNWbIS1JjZ/wplEneAPz02HVI0nZ0xoc8cBNc8jtw+Surbfb5c+Cp1TYp\nSVvUIeTfADecBb+3e7XNfgG4brVNStIWOScvSY0Z8pLUmCEvSY0tNCef5GngB8CrwPGqunKIoiRJ\nw1j0jddXgf1V9cIQxUiShrXodE0GeA5J0pIsGtAFfCXJ15LcNERBkqThLDpdc1VVHUnyJtbD/nBV\nHTx5oyQHphbXqmptwXYlqZUk+4H9Qz/vQiFfVUcmP59Lch9wJfC6kK+qA4u0I0ndTQ5+115bTnL7\nEM8793RNkguSXDi5/0bgA8C3hyhKkjSMRY7k9wD3Tb5seBfw2ap6eJiyJElDmDvkq+qfgX0D1iJJ\nGpinP0pSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVm\nyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDVmyEtSY4a8JDW2a+wCJCk5/yi8vGf1\nLZ93rOqlvatvd3UMeUnbwMt7oEZoNyO8sKyW0zWS1JghL0mNGfKS1JghL0mNGfKS1JghL0mNGfKS\n1JjnyUuaci5JxjhhfSRj/b6ruwjLkJc05RVGuihphDZhxN93ZRdhOV0jSY0Z8pLU2EIhn+SaJH+f\n5B+S3DpUUZKkYcwd8knOAn4fuBr4CeBDSd4xVGE6lbWxC2hkbewCmlkbuwBtYJEj+SuBp6rqmao6\nDnweuHaYsnRqa2MX0Mja2AU0szZ2AdrAIiH/ZuBfp5afnayTJG0THU6hPA5/+io88YPVNvvcOcAF\nq21TkrZmkZD/N+CtU8uXTta9zvIvNjgCHDlvuW2cyljn9+60dpfV9m+N1O4szsR2Z+nPZbS7iHHa\nXdVFWKmar50kZwNPAj/Leso+Bnyoqg4PV54kaRFzH8lX1f8luQV4mPW5/bsMeEnaXuY+kpckbX+L\nnCd/2guhktyY5O8mt4NJ3jX12NOT9d9M8ti8NXQyQ3/+/FSf/W2Sn5l1351owf50fE6ZdXwleU+S\n40mu3+q+O8mC/bn1sVlVW76x/uLwXeBtwDnAIeAdJ23zXmD35P41wKNTj/0TcNE8bXe8zdifF0zd\nvwL47qz77rTbIv05WXZ8bqEvp7b7C+BB4Pqt7LuTbov052T9lsfmvEfym14IVVWPVtVrpzU+yonn\n0Ac/N2faLP3531OLFwL/Puu+O9Ai/QmOz2mzjq+PAPcC35tj351kkf6EOcbmvAN5qxdC/Rrw0NRy\nAV9J8rUkN81ZQycz9WeS65IcBv4M+OhW9t1hFulPcHxO27Qvk/wYcF1V/QEnno/o2Hy9RfoT5hib\nS78YKsn7gQ8D75tafVVVHUnyJtYLPlxVB5ddy5muqu4H7k/yU8AfAW8fuaQz2lR/vo8T+9PxuTW/\nCzjfPpyT+3M66Lc8NucN+ZkuhJq82fop4JqqeuG19VV1ZPLzuST3sf4vzE7+I5r5wjKAqvrrJLuS\nXLzVfXeIrfbnwdf6s6qed3yeYJa+/Eng80kC/Cjwc0n+d8Z9d5p5+/N4VX1xrrE555sHZ/PDNw/e\nwPqbB+88aZu3Ak8B7z1p/QXAhZP7bwS+Cnxg7DdExrzN2J8/PnX/3cA/zrrvTrst2J+Ozy325Unb\n380P33h1bA7bn3ONzbmO5OsUF0IluXn94foU8JvAjwB3Tl6RjlfVlcAe4L7JJb27gM9W1cPz1NHF\njP35C0l+Gfgf4L+AXzzdvqP8ItvEnP35wcnujs8pM/blCbtstu+qat+OFulP5hybXgwlSY15mpgk\nNWbIS1JjhrwkNWbIS1JjhrwkNWbIS1JjhrwkNWbIS1Jj/w+n05nZNBEmvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe2e1941a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# are p-values uniformly distributed if H0 is true?\n",
    "# permutation test version\n",
    "num_runs = 100\n",
    "p_values = np.zeros(num_runs)\n",
    "\n",
    "last = time.time()\n",
    "for i in range(num_runs):\n",
    "    X,_ = sample_gaussian_vs_laplace(n=100)\n",
    "    X2,_ = sample_gaussian_vs_laplace(n=100)\n",
    "\n",
    "    feats_p = RealFeatures(X.reshape(1,len(X)))\n",
    "    feats_q = RealFeatures(X2.reshape(1,len(X2)))\n",
    "    width=1\n",
    "    k = GaussianKernel(10, width)\n",
    "\n",
    "    mmd = QuadraticTimeMMD()\n",
    "    mmd.set_p(feats_p)\n",
    "    mmd.set_q(feats_q)\n",
    "    mmd.set_kernel(k)\n",
    "\n",
    "    mmd.set_num_null_samples(200)\n",
    "    p_values[i] = mmd.compute_p_value(0.05)\n",
    "\n",
    "# ouch: this doesnt look uniform (it has to be)\n",
    "plt.hist(p_values);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.46565\n"
     ]
    }
   ],
   "source": [
    "# does the type1 error (wron alarm) calibration work.\n",
    "# if so, it should be equal to the alpha on average\n",
    "\n",
    "num_runs = 100\n",
    "rejections = np.zeros(num_runs)\n",
    "\n",
    "last = time.time()\n",
    "for i in range(num_runs):\n",
    "    X,_ = sample_gaussian_vs_laplace(n=100)\n",
    "    X2,_ = sample_gaussian_vs_laplace(n=100)\n",
    "\n",
    "    feats_p = RealFeatures(X.reshape(1,len(X)))\n",
    "    feats_q = RealFeatures(X2.reshape(1,len(X2)))\n",
    "    width=1\n",
    "    k = GaussianKernel(10, width)\n",
    "\n",
    "    mmd = QuadraticTimeMMD()\n",
    "    mmd.set_p(feats_p)\n",
    "    mmd.set_q(feats_q)\n",
    "    mmd.set_kernel(k)\n",
    "\n",
    "    mmd.set_num_null_samples(200)\n",
    "    alpha=0.05\n",
    "    rejections[i] = mmd.perform_test()\n",
    "\n",
    "# we expect 0.05 (false) rejection rate here\n",
    "print np.mean(rejections)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Other\n",
    "\n",
    " * I just realised that `CHypothesisTest::compute_p_value` uses this `find_position_to_insert`. This method should go away, and we should just do `np.mean(stat<null_samples)` which automatically gives a quantile\n",
    " \n",
    " * The print out of `KernelManager::kernel_at() : getting the kernel 0` should be removed.\n",
    "\n",
    " \n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}