xval.ipynb.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": [
    "from modshogun import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10,) (10,)\n"
     ]
    }
   ],
   "source": [
    "X,Y = sample_gaussian_vs_laplace(n=10)\n",
    "print X.shape, Y.shape\n",
    "samples_p = RealFeatures(X.reshape(1,len(X)))\n",
    "samples_q = RealFeatures(Y.reshape(1,len(Y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "mmd = QuadraticTimeMMD(samples_p, samples_q)\n",
    "alpha = 0.05\n",
    "sigma = 0.25\n",
    "num_runs = 1\n",
    "train_test_ratio = 4\n",
    "mmd.set_statistic_type(ST_BIASED_FULL)\n",
    "mmd.set_null_approximation_method(NAM_PERMUTATION)\n",
    "mmd.set_variance_estimation_method(VEM_DIRECT)\n",
    "\n",
    "mmd.set_kernel(GaussianKernel(10, 2*sigma**2))\n",
    "mmd.set_train_test_mode(True)\n",
    "mmd.set_train_test_ratio(train_test_ratio)\n",
    "# the following computes statistic on the test-dataset\n",
    "statistic = mmd.compute_statistic()\n",
    "# the following performs test on the test-dataset\n",
    "rejected = mmd.perform_test(alpha)\n",
    "print rejected\n",
    "mmd.set_train_test_mode(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "taus = 2**np.linspace(-5, 5, 5)\n",
    "#taus = np.array([1, 10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# shogun cross validation\n",
    "def shogun_xval(num_samples = 125, num_folds = 5, num_runs = 100, alpha = 0.05):\n",
    "    X,Y = sample_gaussian_vs_laplace(n = num_samples)\n",
    "    samples_p = RealFeatures(X.reshape(1, len(X)))\n",
    "    samples_q = RealFeatures(Y.reshape(1, len(Y)))\n",
    "    \n",
    "    train_test_ratio = num_folds-1\n",
    "\n",
    "    mmd = QuadraticTimeMMD(samples_p, samples_q)    \n",
    "    mmd.set_train_test_mode(True)\n",
    "    mmd.set_train_test_ratio(train_test_ratio)   \n",
    "\n",
    "    for tau in taus:\n",
    "        mmd.add_kernel(GaussianKernel(10, tau))\n",
    "\n",
    "    mmd.set_kernel_selection_strategy(KSM_MAXIMIZE_XVALIDATION, num_runs, alpha)\n",
    "    mmd.select_kernel()\n",
    "    \n",
    "    strategy = mmd.get_kernel_selection_strategy()\n",
    "    measure_vec = strategy.get_measure_vector()\n",
    "    measure_mat = strategy.get_measure_matrix()  \n",
    "    \n",
    "    selected_kernel = GaussianKernel.obtain_from_generic(mmd.get_kernel())\n",
    "    print np.log2(selected_kernel.get_width())\n",
    "\n",
    "    # the following computes statistic on the test-dataset\n",
    "    #statistic = mmd.compute_statistic()\n",
    "    # the following performs test on the test-dataset\n",
    "    #rejected = mmd.perform_test(alpha)\n",
    "    #print rejected\n",
    "    mmd.set_train_test_mode(False)\n",
    "    \n",
    "    return taus, measure_mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "taus, measure = shogun_xval(num_samples = 120, num_folds = 3, num_runs = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEPCAYAAACgFqixAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGodJREFUeJzt3Xu0XWV57/HvLwmIKJeCHgqJMYIgiiCoRByiDUUkYAFP\nbxJ6PEI9bSoXpXQgFy8wjlJA7ZCL1RqKSalVbAUkVFRATT14QBANREwgiMEQ7liiYJUkPP1jzm2W\ni7V23p295ppzvuv3GWMN9lx77pnnYSXr2fN51nynIgIzMxttU+oOwMzM6udiYGZmLgZmZuZiYGZm\nuBiYmRkuBmZmRsXFQNKlkh6WdMc4+1wkaaWkpZL2rTIeMzPrreozg4XAof2+KekwYLeI2B2YD/xD\nxfGYmVkPlRaDiLgR+M9xdjkKuKzc97vAdpJ2qjImMzN7trpnBtOB1R3ba8rnzMxsiOouBmZm1gDT\nav7z1wAv6tieUT73LJK8iJKZ2WaICG1qn2EUA5WPXhYDJwBflHQA8EREPNzvQCkJtZWksyPi7Lrj\nqIrza6+cc4ORyC/pF+lKi4GkzwNzgB0l/RQ4C9gSiIhYEBHXSjpc0j3AU8BxVcbTcLPqDqBis+oO\noGKz6g6gQrPqDqBis+oOoAkqLQYRcUzCPidWGYOZmW2aB8jNsajuACq2qO4AKrao7gAqtKjuACq2\nqO4AmkBtubmNpMh5ZmBmVoXU906fGTSEpDl1x1Al59deOecG+eeXysXAzMzcJjIzy5nbRGZmlszF\noCFy71s6v/bKOTfIP79ULgZmZuaZgZlZzjwzMDOzZC4GDZF739L5tVfOuUH++aVyMTAzM88MzMxy\n5pmBmZklczFoiNz7ls6vvXLODfLPL5WLgZmZeWZgZpYzzwzMzCyZi0FD5N63dH7tlXNukH9+qVwM\nzMzMMwMzs5x5ZmBmZslcDBoi976l82uvnHOD/PNL5WJgZmaeGZiZ5cwzAzMzS+Zi0BC59y2dX3vl\nnBvkn18qFwMzM/PMwMwsZ54ZmJlZMheDhsi9b+n82ivn3CD//FK5GJiZmWcGZmY588zAzMySuRg0\nRO59S+fXXjnnBvnnl8rFwMzMqp8ZSJoLXEBReC6NiPO7vr8j8DlgZ2Aq8HcRsajHcTwzMDOboNT3\nzkqLgaQpwN3AwcADwK3A0RGxomOfs4CtIuIMSS8A7gJ2ioj1XcdyMTAzm6CmDJBnAysj4r6IWAdc\nDhzVtc9DwDbl19sAj3cXglGQe9/S+bVXzrlB/vmlmlbx8acDqzu276coEJ0uAb4h6QHg+cDbK47J\nzMy6VF0MUpwB3B4RB0naDbhe0j4R8WT3jpIWAavKzSeApRGxpPzeHIC2bo8915R4nJ/zG9uOiCVN\nisf5jb9dfn0shVUkqnpmcABwdkTMLbdPB6JziCzpWuCciPhOuf0N4LSI+F7XsTwzMDOboKbMDG4F\nXirpxZK2BI4GFnftsxx4M4CknYA9gHsrjqtxcu9bOr/2yjk3yD+/VJW2iSJig6QTgevY+NHS5ZLm\nF9+OBcC5wEJJtwMC3hcRP6syLjMz+21em8jMLGNNaROZmVkLuBg0RO59S+fXXjnnBvnnl8rFwMzM\nPDMwM8uZZwZmZpbMxaAhcu9bOr/2yjk3yD+/VC4GZmbmmYGZWc48MzAzs2QuBg2Re9/S+bVXzrlB\n/vmlcjEwMzPPDMzMcuaZgZmZJXMxaIjc+5bOr71yzg3yzy+Vi4GZmXlmYGaWM88MzMwsmYtBQ+Te\nt3R+7ZVzbpB/fqlcDMzMzDMDM7OceWZgZmbJXAwaIve+pfNrr5xzg/zzS+ViYGZmnhmYmeXMMwMz\nM0vmYtAQufctnV975Zwb5J9fKhcDMzPzzMDMLGeeGZiZWTIXg4bIvW/p/Nor59wg//xSuRiYmZln\nBmZmOfPMwMzMkrkYNETufUvn11455wb555fKxcDMzKqfGUiaC1xAUXgujYjze+wzB/gEsAXwaEQc\n1GMfzwzMzCYo9b2z0mIgaQpwN3Aw8ABwK3B0RKzo2Gc74P8Db4mINZJeEBGP9TiWi4GZ2QQNZIAs\naaqkb00ijtnAyoi4LyLWAZcDR3XtcwxwRUSsAehVCEZB7n1L59deOecG+eeXatxiEBEbgGfK3943\nx3Rgdcf2/eVznfYAdpD0LUm3SnrHZv5ZZma2maYl7PMksEzS9cBTY09GxHsGGMOrgd8HngfcJOmm\niLhnQMdvhYhYUncMVXJ+7ZVzbpB/fqlSisGV5WNzrAFmdmzPKJ/rdD/wWET8CviVpG8DrwKeVQwk\nLQJWlZtPAEvHXsixUz1ve9vb3h7l7fLrYymsIlHSAFnSc4GZEXFX6oHLn5sK3EUxQH4QuAWYFxHL\nO/bZE7gYmAs8B/gu8PaI+FHXsbIeIEuak/NvKM6vvXLODUYiv8FcgSzpCGAp8LVye19Ji1OCKGcO\nJwLXAXcCl0fEcknzJf1luc8K4OvAHcDNwILuQmBmZtXa5JmBpNso+vlLImK/8rkfRsQrhxBfZxxZ\nnxmYmVVhYGcGwLqIWNv13DObF5aZmTVRSjG4U9IxwFRJu0u6mOIiMRug3D/r7PzaK+fcIP/8UqUU\ng5OAvYBfA18A1gInVxmUmZkNV8rMYLeI+PGQ4hkvDs8MzMwmKPW9M6UY/AfF9QG3Av8P+HZELBtI\nlBPgYmBmNnEDGyBHxO8BL6e4FmB74CuSfjb5EK1T7n1L59deOecG+eeXapNXIEs6EHhj+dge+HeK\nMwQzM8tESptoPXAbcC5wbUQ8PYzAesThNpGZ2QQNcmawPfAG4E3A/hTXGNwUER8cRKCpXAzMzCZu\nkDODJ4B7gZ9QrC+0G0VhsAHKvW/p/Nor59wg//xSpcwM7gVWADcCnwaOq6tVZGZm1UhpE02JiNqX\nn3CbyMxs4ga5NtEukq6S9Ej5uELSjAHEaGZmDZFSDBYCi4Fdysc15XM2QLn3LZ1fe+WcG+SfX6qU\nYvDCiFgYEevLxyLghRXHZWZmQ5QyM/gGxZnAF8qn5lEMkQ+uOLbuODwzMDOboEHODP4c+FPgofLx\nx8BxkwvPzMyaJOU6g/si4siIeGH5eFtE/HQYwY2S3PuWzq+9cs4N8s8vVco9kHeVdI2kR8tPE10t\naddhBGdmZsORMjO4Gfh7Ns4MjgZOiojXVRxbdxyeGZiZTdAg1ya6IyL26Xru9oh41SRjnBAXAzOz\niRvkAPmrkk6XNEvSiyW9D7hW0g6Sdph8qAb59y2dX3vlnBvkn1+qTa5NRPFJIoD5Xc8fDQTg+YGZ\nWcttsk3UFG4TmZlN3CDbRGZmljkXg4bIvW/p/Nor59wg//xSuRiYmVna2kTd6xD1eq5qnhmYmU1c\n6ntn308TSdoK2Bp4gaTfAcYOti0wfSBRmplZI4zXJpoP3AbsWf537HE18MnqQxstufctnV975Zwb\n5J9fqr5nBhFxIXChpJMi4uIhxmRmZkOWMjP4E+BrEfELSR8AXg18JCK+P4wAO+LwzMDMbIIGeZ3B\nB8tCcCDwZuBS4NOTDdDMzJojpRhsKP/7VmBBRHwF2LK6kEZT7n1L59deOecG+eeXKqUYrJH0GeDt\nFAvUPSfx58zMrCVSZgZbA3OBZRGxUtLOwN4Rcd0wAuyIwzMDM7MJGtjMICJ+CTwCHFg+tR5YOYFA\n5kpaIeluSaeNs9/+ktZJ+sPUY5uZ2WCk3PbyLOA04IzyqS2Az6UcXNIUimsSDgX2AuZJ2rPPfucB\nX08LOz+59y2dX3vlnBvkn1+qlN7//wSOBJ4CiIgHgG0Sjz8bWBkR90XEOuBy4Kge+50EfIniDMTM\nzIYspRg8HcVgIQAkPW8Cx58OrO7Yvp+upSwk7QK8LSI+zcYlL0ZORCypO4YqOb/2yjk3yD+/VCnF\n4F/LTxNtL+kvgBuAfxxgDBdQtKHGjGxBMDOryyZvexkRH5d0CPBz4GXAhyLi+sTjrwFmdmzPKJ/r\n9FrgckkCXgAcJmldRCzuPpikRcCqcvMJYOlYVR/r+7V4++TM8nF+zYpvs7c7e+pNiMf5JeVzbJnS\nKhKlfLT0/Ig4bVPP9fnZqcBdwMHAg8AtwLyIWN5n/4XANRFxZY/vZf3RUklzcj5ddX7tlXNuMBL5\nDWw5ikN6PHdYShARsQE4EbgOuBO4PCKWS5ov6S97/UjKcXOU819GcH5tlnNukH9+qfqeGUh6N3A8\nsCvw445vbQN8JyL+V/Xh/VY8WZ8ZmJlVYRBnBp8HjgAWl/8de7xm2IVgFOT+WWfn11455wb555dq\nvPsZrAXWAvOGF46ZmdVhkwPkpnCbyMxs4gY5QDYzs8y5GDRE7n1L59deOecG+eeXysXAzMw8MzAz\ny5lnBmZmlszFoCFy71s6v/bKOTfIP79ULgZmZuaZgZlZzjwzMDOzZC4GDZF739L5tVfOuUH++aVy\nMTAzM88MzMxy5pmBmZklczFoiNz7ls6vvXLODfLPL1Xf+xmYWW8SzwF2B15ePmbCqXdL/EfE6N66\n1drNMwOzPiS2ZeMb/p781ps/q4Dl5eMR4F3A08DfAldF8EwNIZs9S+p7p4uBjTQJATux8Y2+87Ed\ncBfFG/4KNr753xPB013HmUJxW9j3A88HzgW+EMH64WRi1puLQctImhMRS+qOoyp15ycxFXgxvd/0\nN7Dxjb7zsTr1N/yx/Mri8maKojATOB/4pwh+NdiMhqfu165qI5Bf0nunZwaWlR79/LHHHsBjbHyj\nvwX4J2B5BI8O6s8vZwbXA9dLvIGiKHxI4u+Az0Tw1KD+LLNB8pmBtdIE+vljjxURPFlTrPsBZwK/\nB1wEfDKCJ+qIxUaP20TWeoPq5zeFxMuB04E/AD4DXBDBI/VGZblzMWiZEehb9s2v6n7+MEzk9ZN4\nCfA+4O3APwMfi+D+CsOblFH+u5kDzwysceru5zdFBD8B3i3xYeAU4A6JK4DzIvhxvdHZqPKZgQ1c\nm/r5TSCxI/Be4Hjg68C5Efyw3qgsF24TWaVy6+c3QVlEjwdOBm4Czonge/VGZW3nYtAyTe9bSuwC\nHAW8ls3o5zc9v8kaZH4SW1Nc0Xwqxf/PcyL49iCOvXnx+LVrM88MbNIkdgP+sHy8DPgKcCMZ9/Ob\nIIJfAhdLfAZ4B3CpxEMUS118zesfWRV8ZmC/UbZ+XsnGArAT8GXgSmCJWzz1kJgG/AnFtQpe/8gm\nxG0iS1KuqbM/GwvAFhRv/lcCN0WwocbwrIPXP7LN4ZvbtMww11SXmCZxkMTFwE+BhcA6is+9vySC\nUyK4cZCFIPc144eRXwTPRHA18DqKTx+9C7hbYr7EVlX9uX7tRoOLwYiQ2EriDyQ+CzwEfBR4AHhz\nBK+I4AMRfN/96OaLICK4PoI5FDOFI4EfS5wi8bx6o7O2cpsoYxLbAIdRtH/mArdTtH++HMF9dcZm\ng+X1j6wfzwxGVHkB0xEUBWAO8B2KArA4godrDM2GwOsfWbfGzAwkzZW0QtLdkk7r8f1jJN1ePm6U\ntHfVMTXRZPqWErtIHC9xA3AvRTH4IjAzgsMiuKTuQpB7X7Yp+UWwPIJ3UlwPsj2wQuJCiRmbe8ym\n5FaV3PNLVWkxkDQF+CRwKLAXME/Snl273Qu8KSJeBXwEuKTKmHIhsZvEqRI3AT8EDgD+Htg5gj+K\n4F/cJhhdEfwkguMp/t2to1j/6BKJl9YcmjVUpW0iSQcAZ0XEYeX26UBExPl99t8eWBYRL+rxvZFu\nE/kaAJsMr380uprSJpoOrO7Yvr98rp//A3y10ohaRGKKxOskzgfuBq6hWPfnBGB6BH8VwXUuBLYp\nETwewYeAXYFlwA0SV0m8tubQrCEasxyFpIOA44ADx9lnEcWqlwBPAEvH1hQZ6/u1ePtkYCnEjcAb\n4bKT4HffCG95FLgS3vlR+NzKiA2/+XmpUfEn5degeEY5v/Ok6T+AUw+Hk6+UWA4nXQufur3z79fY\n/p099YbEP9Dt3PIrvz62TGkViYbRJjo7IuaW2z3bRJL2Aa4A5kZEz/Xcc24TFRcMnfHXcO7uFJ8Z\nX0XR/rkqguW1BjcgI7AYWCvzk9iS4lqF04GHgXPoWv+orbmlGoH86v9oqaSpFEsZHww8SHHTknkR\nsbxjn5nAN4B3RMTN4xwrq2LQ5xqAqygKgK8BsKHy+kf5akQxKAOZC1xIMZ+4NCLOkzSf4gxhgaRL\nKN4Q7wMErIuI2T2O0/pi4GsArOn6rH90eQTrag3MNltjisGgtLUYlPcBeBtFAdgfuIGiAHyl86Of\nI3Cq6vxapPz02sHA++HrL4ND/y+wKIJf1RzawOX22nVryqeJRlKPawBej68BsBYp1z+6IYKD4Kpz\nKM4WvP5RxnxmMAB9rgG4muIM4Fv+6KfloGP9ozcBF+P1j1rBbaKK+T4ANqok9qT49NERwALgE17/\nqLncJqpAlfcByH19FOfXXt25RbAigmOB11BcBDnp9Y/qlPNrNxEuBpvQ4z4AH6P4mOwhvg+AjbII\nVnn9o3y4TdTzz/I1AGYT1WP9o4sorjNa61+W6uOZwaT+LI6huCrT1wCYTZDEtsC7Kf4NzQSmAmso\n1iZb0+frhz1nq4aLQcuMwGednV9LTTa38kx7OjCj/G+vr3cAHqF/wVgDrIngl5udSN/48n3tIP29\nszEL1ZlZniL4BbCifPQksQWwM88uFPt1fD1d4peMf5axBnjcbamJ85mBmbVCeT3PjvQ/uxj7+rnA\nA4xfNB4clSU23CYys5EksTUdZxP0Lhr/A3ic/jOMsbbUL4Yd/6C5GLTMCPQtnV9L5ZhbuUrrTsB0\n+OBb4MNr6V1A1rOJggE82uTVXT0zMDPrI2Ljm7z0ka0jPryke5+yLbUdzz6z2Bd4a8f2thIPMn7R\neCCCX1ec1qT4zMDMbBIkngPsQv8ZxnSK4fhaisJwWQSfGF58bhOZmTVCuZbZCykKw38N8w6GXpuo\nZXJfH8X5tVfOucFw8ovgmQgeLpeuaeStbF0MzMzMbSIzs5y5TWRmZslcDBrCfdl2yzm/nHOD/PNL\n5WJgZmaeGZiZ5cwzAzMzS+Zi0BC59y2dX3vlnBvkn18qFwMzM/PMwMwsZ54ZmJlZMheDhsi9b+n8\n2ivn3CD//FK5GJiZmWcGZmY588zAzMySuRg0RO59S+fXXjnnBvnnl8rFwMzMPDMwM8uZZwZmZpas\n8mIgaa6kFZLulnRan30ukrRS0lJJ+1YdUxPl3rd0fu2Vc26Qf36pKi0GkqYAnwQOBfYC5knas2uf\nw4DdImJ3YD7wD1XG1GC5F0Hn11455wb555ek6jOD2cDKiLgvItYBlwNHde1zFHAZQER8F9hO0k4V\nx9VE29cdQMWcX3vlnBvkn1+SqovBdGB1x/b95XPj7bOmxz5mZlYhD5CbY1bdAVRsVt0BVGxW3QFU\naFbdAVRsVt0BNMG0io+/BpjZsT2jfK57nxdtYh+g+IjUQKNrGEnvrDuGKjm/9so5N8g/vxRVF4Nb\ngZdKejHwIHA0MK9rn8XACcAXJR0APBERD3cfyNcYmJlVp9JiEBEbJJ0IXEfRkro0IpZLml98OxZE\nxLWSDpd0D/AUcFyVMZmZ2bO15gpkMzOrTusGyJJOkrRc0jJJ59UdTxUk/Y2kZyTtUHcsgyTpo+Vr\nt1TSFZK2rTumyUq5qLKtJM2Q9E1Jd5b/3t5Td0yDJmmKpO9LWlx3LFWQtJ2kfyv/3d0p6XX99m1V\nMSivFDwC2Dsi9gY+Xm9EgydpBnAIcF/dsVTgOmCviNgXWAmcUXM8k5JyUWXLrQdOiYi9gNcDJ2SW\nH8B7gR/VHUSFLgSujYiXA68ClvfbsVXFAHg3cF5ErAeIiMdqjqcKnwBOrTuIKkTEDRHxTLl5M8Un\nx9os5aLK1oqIhyJiafn1kxRvJNlcA1T+4nU48I91x1KF8sz7jRGxECAi1kfEz/vt37ZisAfwJkk3\nS/qWpNfWHdAgSToSWB0Ry+qOZQj+HPhq3UFMUspFlVmQNIti2Ybv1hvJQI394pXr4PQlwGOSFpat\nsAWSnttv56o/Wjphkq4HOpejEMWL9QGKeH8nIg6QtD/wr8Cuw49y820ivzMpWkSd32uVcfJ7f0Rc\nU+7zfmBdRHy+hhBtgiQ9H/gS8N7yDKH1JL0VeDgilpbt59b9W0swDXg1cEJEfE/SBcDpwFn9dm6U\niDik3/ck/RVwZbnfreWQdceIeHxoAU5Sv/wkvZLiSsjbJYmihXKbpNkR8cgQQ5yU8V4/AEnHUpya\n//5QAqpWykWVrSZpGkUh+OeIuLrueAboDcCRkg4HngtsI+myiPjfNcc1SPdTdBq+V25/Cej7IYe2\ntYm+TPkmImkPYIs2FYLxRMQPI+J3I2LXiHgJxQu5X5sKwaZImktxWn5kRPy67ngG4DcXVUrakuKi\nytw+lfJZ4EcRcWHdgQxSRJwZETMjYleK1+2bmRUCyot3V5fvlQAHM86wvHFnBpuwEPispGXAr4Gs\nXrwuQX6nrhcDWwLXFyc/3BwRx9cb0ubrd1FlzWENjKQ3AH8GLJP0A4q/k2dGxNfqjcwm4D3Av0ja\nAriXcS7q9UVnZmbWujaRmZlVwMXAzMxcDMzMzMXAzMxwMTAzM1wMzMwMFwOzcZXrueS2UqfZs/g6\nAzMz85mB2RhJW0v6d0k/kHSHpD8tV8d9dfn9d0m6q1w1d4Gki8rnF0r6lKSbJN0jaY6kRZJ+JOmz\nHcf/lKRbyhvF9FwszKwuLgZmG80F1kTEfhGxD/CbZRck7UyxsuxsikXOultH20fE64FTKNYn+mhE\nvALYR9I+5T5nRsRsipuMzCkXJzRrBBcDs42WAYdIOlfSgV03ApkNLImItRGxAfi3rp+9puMYD0bE\n2IJgd1KsRgtwtKTbgB8ArygfZo3QtoXqzCoTESvLltDhwIclfZPfvvHJeAsHjq3C+kzH12Pb08qb\nw/wN8JqI+LmkhcBWg4rdbLJ8ZmBWKltB/1XedOfjFDcGGXMrxV32tivX+P+j8Q7V47ltgSeBX0ja\nCThsQGGbDYTPDMw22hv4mKRngKcp7rn9cYCIeEDS3wK3AD8DVgBry5/r/khedH8dEXdIWkpxH+HV\nwI1VJWG2OfzRUrNEkp4XEU9JmgpcRXH/gpzu/mUjzG0is3Rnlzd5WQbc60JgOfGZgZmZ+czAzMxc\nDMzMDBcDMzPDxcDMzHAxMDMzXAzMzAz4by8Wq1hY8y0YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fae9524e950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_powers = np.mean(measure, 0)\n",
    "plt.plot(np.log2(taus), test_powers)\n",
    "plt.xlabel(\"sigma\")\n",
    "plt.ylabel(\"test power\")\n",
    "plt.grid(True)\n",
    "plt.ylim([0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_test_power_xval(width = 2*0.5**2, num_runs = 10, num_folds = 3, num_samples = 300, alpha = 0.05):\n",
    "    X,Y = sample_gaussian_vs_laplace(num_samples)\n",
    "    test_power = []\n",
    "    for _ in range(num_runs):\n",
    "        inds_x = np.random.permutation(len(X))\n",
    "        inds_y = np.random.permutation(len(Y))    \n",
    "\n",
    "        for k in range(num_folds):\n",
    "            binsize = num_samples/num_folds\n",
    "            train_inds_x = inds_x[(k*binsize):((k+1)*binsize)]\n",
    "            train_inds_y = inds_y[(k*binsize):((k+1)*binsize)]\n",
    "                        \n",
    "            train_inds_x = np.setdiff1d(inds_x, train_inds_x)\n",
    "            train_inds_y = np.setdiff1d(inds_y, train_inds_y)\n",
    "            \n",
    "            X_train = X[train_inds_x]\n",
    "            Y_train = Y[train_inds_y]        \n",
    "\n",
    "            samples_p = RealFeatures(X_train.reshape(1, len(X_train)))\n",
    "            samples_q = RealFeatures(Y_train.reshape(1, len(X_train)))\n",
    "\n",
    "            mmd=QuadraticTimeMMD(samples_p, samples_q)            \n",
    "            mmd.set_kernel(GaussianKernel(10, width))\n",
    "            test_power += [mmd.perform_test(alpha)]\n",
    "\n",
    "    return test_power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.056666666666666664, 0.70499999999999996, 0.040000000000000001, 0.021666666666666667, 0.013333333333333334]\n"
     ]
    }
   ],
   "source": [
    "# compute test power for one kernel\n",
    "\n",
    "test_powers = []\n",
    "for tau in taus:\n",
    "    test_power = get_test_power_xval(width = tau, num_runs=200, num_samples = 150)\n",
    "    test_powers += [np.mean(test_power)]\n",
    "print test_powers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEPCAYAAACgFqixAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHVW1/vHvGwaZZ+UqCCGAoAgiMl4Bg4gEvMJVrkpw\nAkVxAOSCjOEKz0+BMBPQi4CQKIiogEyCMkhUlHmWeRBkEkFvmBRIYP3+qGpoOj1Ud1edXWf3+3me\n8yR1Tp06a9HhrK61q/ZWRGBmZmPbuNQBmJlZei4GZmbmYmBmZi4GZmaGi4GZmeFiYGZmNFwMJJ0q\n6UlJtw2yz/GS7pN0i6S1m4zHzMz61/SZwXRgy4FelLQVsHJErArsAny/4XjMzKwfjRaDiLgK+L9B\ndtkW+FG577XA4pKWbTImMzObW+oxg+WAR3ptP1Y+Z2ZmHZS6GJiZWQvMm/jzHwPe3mt7+fK5uUjy\nJEpmZiMQERpqn04UA5WP/lwAfB34qaQNgVkR8eRAB6qSULeSdHBEHJw6jqY4v+6Vc24wJvKr9It0\no8VA0pnARGBpSX8BDgLmByIiTo6IiyVtLel+4AVgpybjabnxqQNo2PjUATRsfOoAGjQ+dQANG586\ngDZotBhExA4V9tm1yRjMzGxoHkBujxmpA2jYjNQBNGxG6gAaNCN1AA2bkTqANlC3LG4jKXIeMzAz\na0LV706fGbSEpImpY2iS8+teOecG+edXlYuBmZm5TWRmljO3iczMrDIXg5bIvW/p/LpXzrlB/vlV\n5WJgZmYeMzAzy5nHDMzMrDIXg5bIvW/p/LpXzrlB/vlV5WJgZmYeMzAzy5nHDMzMrDIXg5bIvW/p\n/LpXzrlB/vlV5WJgZmYeMzAzy5nHDMzMrDIXg5bIvW/p/LpXzrlB/vlV5WJgZmYeMzAzy5nHDMzM\nrDIXg5bIvW/p/LpXzrlB/vlV5WJgZmYeMzAzy5nHDMzMrDIXg5bIvW/p/LpXzrlB/vlV5WJgZmYe\nMzAzy5nHDMzMrDIXg5bIvW/p/LpXzrlB/vlV5WJgZmYeMzAzy5nHDMzMrDIXg5bIvW/p/LpXzrlB\n/vlV5WJgZmbNjxlImgQcR1F4To2Iw/u8vjRwBvBWYB7g6IiY0c9xPGZgZjZMVb87Gy0GksYB9wKb\nA48D1wPbR8TdvfY5CFggIvaXtAxwD7BsRMzpcywXAzOzYWrLAPL6wH0R8XBEzAbOArbts89fgUXL\nvy8K/L1vIRgLcu9bOr/ulXNukH9+Vc3b8PGXAx7ptf0oRYHo7RTgCkmPA4sAn2o4JjMz66PpYlDF\n/sCtEbGZpJWByyStFRHP991R0gzgoXJzFnBLRMwsX5sI0K3bPc+1JR7n5/x6tiNiZpvicX6Db5d/\n35HCQ1TU9JjBhsDBETGp3N4PiN6DyJIuBg6JiD+U21cA+0bEDX2O5TEDM7NhasuYwfXAKpJWlDQ/\nsD1wQZ997gI+BCBpWeAdwIMNx9U6ufctnV/3yjk3yD+/qhptE0XEK5J2BS7l9UtL75K0S/FynAwc\nBkyXdCsgYJ+I+EeTcZmZ2Rt5biIzs4y1pU1kZmZdwMWgJXLvWzq/7pVzbpB/flW5GJiZmccMzMxy\n5jEDMzOrzMWgJXLvWzq/7pVzbpB/flW5GJiZmccMzMxy5jEDMzOrzMWgJXLvWzq/7pVzbpB/flW5\nGJiZmccMzMxy5jEDMzOrzMWgJXLvWzq/7pVzbpB/flW5GJiZmccMrDqJcRG8mjoOM6vOYwbWhAsl\ndkkdhJnVz8WgJdret5TYENgIOEBivuG/v935jVbO+eWcG+SfX1UuBlbVFOBA4AFgcuJYzKxmHjOw\nIUmsDVwMTAA2AY4D1vT4gVn7eczA6nQAcHQELwKXAy8C/5E2JDOrk4tBS7S1bymxOjAROAkgggCm\nAvtLVD5Ta2t+dck5v5xzg/zzq8rFwIayP3B8BM/3eu5cYGlg0zQhmVndPGZgA5JYCbgBWDmCWX1e\n2xnYLoKtkgRnZpV4zMDqsC9wUt9CUDodWKscXDazLudi0BJt61tKLAd8Eji2v9cjeAk4Btiv2vHa\nlV/dcs4v59wg//yqcjGwgewF/DCCpwbZ52TgQxKrdCgmM2uIxwxsLhJvBu6huJfgsSH2/Tbw5gi+\n0pHgzGxYqn53uhjYXCQOAZau8gXfq3CsEcETjQdnZsPiAeQu05a+pcQSwFeAw6vsX7aRzgD2GPy4\n7civKTnnl3NukH9+VbkYWF+7AhdF8OdhvOdoYOeykJhZF3KbyF4jsQjwILBpBHcP870/Au6O4NBG\ngjOzEfGYgQ2bxF7ABhF8cgTvXQO4Algpgn/VHpyZjYjHDLpM6r6lxAIUl5OO6Df7CO4ArgV26v/4\nefdlc84v59wg//yqcjGwHl8AborgllEcYyqwt8S8NcVkZh3SeJtI0iSK+e/HAadGxFxXqZSV+Vhg\nPuCpiNisn33cJmpIuXLZfcDkCK4e5bFmAqdE8OM6YjOz0WnFmIGkccC9wObA48D1wPYRcXevfRYH\n/gh8OCIek7RMRDzdz7FcDBoisSPw2Qg2r+FYk4AjgbXK6a7NLKFaxgwkzSPpylHEsT5wX0Q8HBGz\ngbOAbfvsswNwTkQ8BtBfIRgLUvUtJeahWLzmOzUd8tfAHGDrN35O3n3ZnPPLOTfIP7+qBi0GEfEK\n8Gr52/tILAc80mv70fK53t4BLCXpSknXS/rsCD/LRua/gKeBmXUcrPfiN3Ucz8w6o8pA3/PA7ZIu\nA17oeTIidq8xhnWADwILA1dLujoi7q/p+F0hImZ2+jMlxlEsdL9/zS2dc4BDJDaO4CpIk18n5Zxf\nzrlB/vlVVaUYnFs+RuIxYIVe28uXz/X2KPB0RLwIvCjpd8B7gLmKgaQZwEPl5izglp4fZM+pnreH\ns73/v8OhrwAX13n8COZIR50PKx0J223Unny97e38t8u/70jhISqqNIAsaUFghYi4p+qBy/fNQzGJ\n2ebAE8B1wOSIuKvXPqsDJwCTgDdRXKv+qYi4s8+xsh5AljSxk7+hlOsXXwMcGcHZDRx/AYq7mSdF\ncFun8+u0nPPLOTcYE/nVc9OZpI8CtwC/KrfXlnRBlSDKMYddgUuBO4CzIuIuSbtI+nK5z90Ug463\nUXw5ndy3EFgjPgQsysjP+gYVwYsUlxTv28TxzaxeQ54ZSLqRop8/MyLeWz73p4h4dwfi6x1H1mcG\nnVbeD3BqBKc3+BmLUZwdrB/Bg019jpkNrLYzA2B2RDzT57lXRxaWtYHExhRjOT9p8nMieJZiNbRv\nNvk5ZjZ6VYrBHZJ2AOaRtKqkEyhuErMadfha5ynA1AjmdOCzpgHbSxt8rAOflUzO16rnnBvkn19V\nVYrBbsAawEsUv0k+wxALmVh7SbwPWBP4YSc+L4IngZ/Al7brxOeZ2chUGTNYOSIe6FA8g8XhMYMa\nSJwL/DaCaR38zJUopiJZOYK+LUcza1DV784qxeC3FPcHXA/8HvhdRNxeS5TD4GIwer3WHJgQwT87\n/NlnALdHVFtO08zqUdsAckR8AHgnxb0ASwC/lPSP0YdovXWob7k/MK3ThaAw+Qpgj/L+g+zk3HfO\nOTfIP7+qhrwDWdLGwCblYwngIoozBOsiEqtQ3Nj3tTQRnPVn+MmNwOeBk9LEYGYDqdImmgPcCBwG\nXBwRL3cisH7icJtoFCROAZ6I4FsJY3g/8CNgtQ5dyWQ25tU5ZrAE8H5gU2A9insMro6I/6kj0Kpc\nDEZO4u3ArcCqEfw9cSy/B74XwVkp4zAbK+ocM5hFcRfpnynmF1qZojBYjRruW+5NcbdxskLQK7+p\nwH7l3EjZyLnvnHNukH9+VVWZm+hB4GhgKeBEYLVyUNm6gMSywGeAY1LHUrqY4t/dpNSBmNnrqrSJ\nxkVE8ukn3CYaGYnDgYUj2DV1LD0kdgB2icC/VJg1rM65id4m6ReS/lY+zpG0fA0xWsMklgJ2Bo5I\nHUsfPwNWkPj31IGYWaFKMZgOXAC8rXxcWD5nNWqob7kbcF4Ef2ng2MPSO7/ySqIjyWh665z7zjnn\nBvnnV1WVYvDmiJgeEXPKxwzgzQ3HZaMksSjFWhJTU8cygOnABhIdnQrdzPpXpRj8XdJnJM1TPj4D\naS9PzFEDKy19Fbg8gvtqPu6I9M0vgn9RzGi6T5KAapbzSlk55wb551dVlQHkFSmmotiofOoPwO4R\n0dHWgweQq5NYkOJS4C0i6Pg8UlVJLAE8AKwTwcOp4zHLUZ33GTwcEdtExJvLx392uhCMBTX3LXcG\nrmlTIegvvwhmAaeQweI3Ofedc84N8s+vqir3GUyQdKGkp8qric6XNKETwdnwScxPcZPZIaljqWga\n8GmJt6QOxGwsq9Imugb4Hq8vkbg9sFtEbNBwbH3jcJuoAomdgU9EsGXqWKqS+D7wVAQdneLEbCyo\nc26i2yJirT7P3RoR7xlljMPiYjA0iXmBu4EvRPC71PFUJbEycC3FOgvPpo7HLCd13nR2iaT9JI2X\ntKKkfYCLJS0laanRh2pQW9/yU8DjbSwEg+UXwQPAZcCXOxZQzXLuO+ecG+SfX1VDrmcAfLL8c5c+\nz28PBODxgxaQGAccAOyZOpYROhz4pcQJEbyUOhizsWbINlFbuE00OImPAVOA9SLojh9qHxKXAOdE\n8IPUsZjlorYxg7ZwMRhYOR30DcC3IzgvdTwjJbEp8APgnRG8kjoesxzUOWZgHTDKvuWWwJso5pBq\npYr5/R54Gvh4s9HUL+e+c865Qf75VeVi0OXKs4IDgUMiSD7V+GiU7a0sF78xa7sql5ZeERGbD/Vc\n09wm6p/EByhaK6vn0FopB8JvA/aM4NLU8Zh1u6rfnQNeTSRpAWAhYBlJS8Jrv6ktBixXS5RWhynA\n1BwKAUAEr5YL8uwHLgZmnTJYm2gX4EZg9fLPnsf5wHebD21sGUnfUmJ9ip/P6bUHVLNh5ncWMEGi\no3e5j0bOfeecc4P886tqwDODiJgGTJO0W0Sc0MGYrLopwBERvJw6kDpFMFviKIqzg4+ljsdsLKgy\nZvAJ4FcR8ZykA4F1gO9ExE2dCLBXHB4z6EViLeDXFFM4/Ct1PHWTWIhiGu6JEdyVOh6zblXnpaX/\nUxaCjYEPAacCJ442QBu1A4BjciwEABH8k2IdjSwWvzFruyrFoGdg8iPAyRHxS2D+5kIam4bTt5RY\nDdgc+H5jAdVshH3Z7wHbSKxQczi1y7nvnHNukH9+VVUpBo9JOoliErSLJb2p4vusOfsBJ0TwXOpA\nmhTB/wGn0b3zLZl1jSpjBgsBk4DbI+I+SW8F1oyIjl725zGDgsSKwE3AKuWXZdYk3gb8CXhHBE+n\njses29S57OU/gb8BG5dPzYHqi6xLmiTpbkn3Stp3kP3WkzRbUtdNRdBh+wAnj4VCABDB48A5wK6p\nYzHLWZVlLw8C9gX2L5+aDzijysEljaO4J2FLYA1gsqTVB9hvKsXVMWNSlb6lxFuBycCxjQdUs1H2\nZY8Avi6xSE3h1C7nvnPOuUH++VVVpff/MWAb4AWAiHgcWLTi8dcH7ouIhyNiNsXNRNv2s99uwNkU\nZyA2sL2A0yPG1n+nCO4DrgS+lDoWs1xVKQYvRzGwEACSFh7G8ZcDHum1/Sh9prKQ9DbgPyPiRBi7\nk5NFxMzBXpdYBvgCcGRHAqrZUPlVMBXYU2rnlWw15NdaOecG+edXVZVi8LPyaqIlJH0JuBxqXXzk\nOIo2VI8xWxCG8A3g7AgeTR1IChHcBNwJfCZ1LGY5GnLZy4g4StIWwLPAasC3IuKyisd/DN5wjfjy\n5XO9rQucJUnAMsBWkmZHxFxz80uaATxUbs4Cbump6j19vy7e3mOgfCQWh8t3h6O/ApcwwuOn3h4w\nv+rbu18C0/aR+CFok/zya+d27556G+JxfpXy2bFM6SEqqnJp6eERse9Qzw3w3nmAeyhukHoCuA6Y\nHBH9Ti8gaTpwYUSc289rWV9aKmniQKerEgdQrP712c5GVZ/B8qt+DARcAxwewVz/RlKqI7+2yjk3\nGBP5VfrurFIMboqIdfo8d1tErFUxkEnANIqW1KkRMVXSLkBExMl99j0NuGgsFoOBSCwMPAhsFsGd\nqeNJrVzr+QBg/W5d69msk0ZdDCR9FfgaMAF4oNdLiwJ/iIiO9m7HcDHYA9g4gv9KHUsblIvf3AHs\nGsEVqeMxa7s6isHiwJLAYRTTH/R4LiL+UUuUw5B7MejvVFXiTRRnBf8Rwc1JAqtJnafiEjsCn45g\nizqOV4ecWw055wZjIr9K352DrWfwDPAMxU1OlsaOwC3dXggacCbw/yTWjeCG1MGY5WDIMYO2yP3M\noC+J+YB7KX4D/mPqeNpG4hvAJm6fmQ2u6nenZx9tr8nAQy4EA/oBsGk5nbeZjZKLQUv0vtZZYh6K\nuaC+kyygmtU9/0sEL1Csd7B3nccdqZznt8k5N8g/v6pcDNrp4xTjNb9JHUjLfRf4uMTyqQMx63Ye\nM2iZ8saqm4EDI7godTxtJ3EMEBHslToWszbymEH3+kj55y+TRtE9jgF2klgqdSBm3czFoCUkTSzP\nCqYAh+Z2d21Tfdly4r7zSLz4Tc5955xzg/zzq8rFoF02o7jR75zUgXSZI4Bdy6k7zGwEPGbQIhK/\nAX4YwQ9Tx9JtJM4BfhvB8aljMWuT2iaqa4vci4HEvwM/plj4fXbqeLqNxHoUq+Wt4v9+Zq/zAHLX\nOfdoiqmZs/wia7ovG8H1wH3ADk1+zkBy7jvnnBvkn19VLgYtIPFeWGQVYEbqWLrcVGDfcmZTMxsG\nt4laQOLnwB8jODZ1LN2svBrreuDbEZyfOh6zNvCYQZeQeCcwE5hQTrFgoyCxHcUUFRvldnmu2Uh4\nzKB77A9MA62XOpAmdbAvex7F5bkf6NDnAXn3nXPODfLPryoXg4QkJgBbU0y4ZjWI4BWK+w72G2pf\nM3ud20QJSZwEPBXBgaljyUm5QtwDwDYR3JQ6HrOUPGbQcuVMm7dR3FfwdOp4ciOxJ7BBBJ9KHYtZ\nSh4zaL9vAtN7CkHufcsE+Z0CfFBi1U58WM4/v5xzg/zzq8rFIAGJtwCfA45OHUuuIngOOJGi6JrZ\nENwmSkDiUGCJCL6WOpacSSxDsY70uyN4PHU8Zil4zKClJJYE7gfeF8FDicPJnsQ04KUI9kkdi1kK\nHjNor12BC/oWgtz7lgnzOxr4YlmEG5Pzzy/n3CD//KpyMeggiUWA3Sjm0LEOiOAvwIXglpzZYNwm\n6iCJbwLr+XLHzpJ4F3AlsFIE/0wdj1knecygZSQWpLgRalIEt6WOZ6yROA+4LMJ3e9vY4jGD9vkC\ncMNAhSD3vmUL8psK7C0xXxMHb0F+jck5N8g/v6pcDDqg/ALaBzgkdSxjVQTXAH8Gt+jM+uM2UQdI\n7ATsEMEWqWMZyyS2pLi6aK0IXk0dj1knuE3UEhLzUExT7bOC9C4FXgY+kjoQs7ZxMWjeJ4C/Ab8d\nbKfc+5ZtyK9c7GYqsH+5Klpt2pBfU3LODfLPryoXgwaVa/FOAQ7xqlutcQ7wFmDj1IGYtYnHDBok\nsS3wLWBdF4P2kPgysG2E20WWP99nkFjZhrgWmBrBuanjsdeVi9/8GdgqgltTx2PWpNYMIEuaJOlu\nSfdK2ref13eQdGv5uErSmk3H1CFbAAtTrMk7pNz7lm3KL4KXgGOBuf49jlSb8qtbzrlB/vlV1Wgx\nkDQO+C6wJbAGMFnS6n12exDYNCLeA3yHYlGSHEwBDvUljK11EvDhch1qszGv0TaRpA2BgyJiq3J7\nPyAi4vAB9l8CuD0i3t7Pa13TJpLYBJgBrBbBnMTh2AAkvgMsHcFXU8di1pS2tImWAx7ptf1o+dxA\ndgYuaTSizpgCHOZC0HrHA5+S+LfUgZilNm/qAHpI2gzYiUEu+ZM0A15bB2AWcEtEzCxfmwiQehvi\neWANWPIoadbEYbx/jzbmU+N2K/ODOBP4hqRf55hfHdu9e+ptiMf5VcpnxzKlh6ioE22igyNiUrnd\nb5tI0loU139PiogHBjhWV7SJJH4BXBnB8cN7n14rHDlqa34S44EbgQkRPDPy47QzvzrknBuMifzS\nX1oqaR7gHmBz4AngOmByRNzVa58VgCuAz0bENYMcq/XFQOLdwOUUXyyeN79LSJwO3BHhRYcsP60o\nBmUgk4BpFOMTp0bEVEm7UJwhnCzpFODjwMOAgNkRsX4/x+mGYnAmcGsE/Q6QWztJrEkxb9GECP6V\nOh6zOrWmGNSl7cVAYhXgamDlCJ4d/vuzP1VtdX4SFwK/jOD7I3t/u/MbjZxzgzGRXyuuJhpL9gO+\nN5JCYK1wGLCP1J6LKsw6yWcGNZBYAbgZWDWCf6SOx0ZG4nfAiRH8JHUsZnXxmUFn7Q38wIWg600F\n9qt7emuzbuBiMErlDUufBo4Z3XHynh+lS/LrueFxq+G+sUvyG5Gcc4P886vKxWD09gTOiODJ1IHY\n6PRa/Ga/1LGYdZrHDEZBYmngXmDtiDdMu2FdqhxAvgf4XAR/SB2P2Wh5zKAzdgd+4UKQj3I+qSPx\n2YGNMS4GIySxGPA1qOeu1dz7ll2W3wxg3fJmtEq6LL9hyTk3yD+/qlwMRu6rwKUR3J86EKtXBC9S\n3DW/T+pYzDrFYwYjILEQxaI8m0dwR+p4rH4SiwMPUKxf/VDicMxGzGMGzdoZ+KMLQb7KGUxPAb6Z\nOhazTnAxGKZyMfW9gUPqPW7efcsuzW8asIPEW4basUvzqyTn3CD//KpyMRi+zwF/iuDG1IFYsyL4\nK/BTiqvGzLLmMYNhxfDaNeifj+CqlLFYZ0hMoFiHY4InIbRu5DGDZmwPPOJCMHZE8CDFWge7pI7F\nrEkuBhVJjAMOoOaxgtePn3ffssvzOxz4b4kFBtqhy/MbVM65Qf75VeViUN3HgOcplrW0MSSCWymm\nKP9c6ljMmuIxg0qfjSgWTT84ggtSxGBpSWwCTAdWi+CV1PGYVeUxg3pNAuYFLkodiCVzFfAksF3q\nQMya4GIwhPKs4EDgkAhebe5z8u5bdnt+5fTWhzHA4jfdnt9gcs4N8s+vKheDoX0AWAY4O3UgltzF\nwHzAh1MHYlY3jxkM+blcBpwZwfROf7a1j8RngC9GsFnqWMyq8JhBDSQ2AFYFzkgdi7XGT4GVJDZM\nHYhZnVwMBjcFOCKC2U1/UO59y1zyK/8tzLX4TS759Sfn3CD//KpyMRiAxHuAdYHTUsdirTMd2Eji\nXakDMauLxwwG/Dx+ClwXwdGd+kzrHhJTgFUj2DF1LGaDqfrd6WLQ72exGvB7isnJnu/EZ1p3kViC\nYvGb90bwl9TxmA3EA8ijszxwUCcLQe59y9zyi2AWcCqwF+SXX2855wb551eVi0E/IrgighNTx2Gt\ndxzwWYllUgdiNlpuE5mNgsRJwJMRfCt1LGb9qfrdOW8ngjHL2JHAtRJLAc9UeDxfTm1h1io+M2gJ\nSRMjYmbqOJqSc34SH4LjPgp7PAUsPsRjAeA5qhWOwQpKY/NkzZ1fvj87GBP5+czArBMiuFz67zkR\ne8wcat9y6dTFGLxgLA1MKP++RD+vLyiNuqA818mCYu3nMwOzLlOxoAz1WIhisaZZuKBkzfcZmNmA\nJOZh9AVlYYqCMlCxeBZ4qXy83OvvLw3x/ID7erxl+FpTDCRNorgEbxxwakQc3s8+xwNbAS8AO0bE\nLf3sk3UxGAN9S+fXpQbKrSwoi/LGAtG7rbUo8KZej/n7bA/03EDPzwfMpmLhqPjcS3DMirDnn0Z6\n3AjmjOS/a6e0YsxA0jjgu8DmwOPA9ZLOj4i7e+2zFbByRKwqaQPg+zAmZ4RcG5iZOogGOb/u1W9u\n5fKfs8pH48pFhaoWlCpFZsniz4fWZe7CVfm4Kr5mh3OGc0kb72NqegB5feC+iHgYQNJZwLbA3b32\n2Rb4EUBEXCtpcUnLRsSTDcfWNkukDqBhzq97tSK3skXU84VaG+mEgyOOP3jk72dehleQHhllyI1o\nuhgsxxsTf5SiQAy2z2Plc2OtGJhZFyrbRHOAf6aOZTQ8HUV7jE8dQMPGpw6gYeNTB9Cg8akDaNj4\n1AG0QdNnBo8BK/TaXr58ru8+bx9iH6AYCKk1upaR9PnUMTTJ+XWvnHOD/POroulicD2wiqQVgSeA\n7YHJffa5APg68FNJGwKz+hsvyPlKIjOz1BotBhHxiqRdgUt5/dLSuyTtUrwcJ0fExZK2lnQ/xaWl\nOzUZk5mZza1rbjozM7PmdN0AsqTdJN0l6XZJU1PH0wRJe0l6VdJSqWOpk6Qjyp/dLZLOkbRY6phG\nS9IkSXdLulfSvqnjqZOk5SX9RtId5f9vu6eOqW6Sxkm6SdIFqWNpQnmp/s/L/+/uKO/l6ldXFYNy\nRaKPAmtGxJrAUWkjqp+k5YEtgIdTx9KAS4E1ImJt4D5g/8TxjEqvmyq3BNYAJktaPW1UtZoD7BkR\nawAbAV/PLD+AbwB3pg6iQdOAiyPincB7gLsG2rGrigHwVWBqRMwBiIinE8fThGOBvVMH0YSIuDwi\neiY2u4biyrFu9tpNlRExG+i5qTILEfHXnqlhIuJ5ii+S5dJGVZ/yF6+tgR+kjqUJ5Zn3JhExHSAi\n5kTEswPt323F4B3AppKukXSlpHVTB1QnSdsAj0TE7alj6YAvAJekDmKU+rupMpsvy94kjaeYluLa\ntJHUqucXr1wHTlcCnpY0vWyFnSxpwYF2bt16BpIuA5bt/RTFD+tAiniXjIgNJa0H/Ixi3veuMUR+\nB1C0iHq/1lUGyW9KRFxY7jMFmB0RZyYI0YZJ0iLA2cA3yjOErifpI8CTEXFL2X7uuv/XKpgXWAf4\nekTcIOk4YD/goIF2bpWI2GKg1yR9BTi33O/6cpB16Yj4e8cCHKWB8pP0boo7IW+VJIoWyo2S1o+I\nv3UwxFEZ7OcHIGlHilPzD3YkoGZVuamyq0mal6IQnB4R56eOp0bvB7aRtDWwILCopB9FxOcSx1Wn\nRyk6DTeU22cDA17k0G1tovMov0QkvQOYr5sKwWAi4k8R8W8RMSEiVqL4Qb63mwrBUMrpzPcGtomI\nWicbS+Qr0co1AAACX0lEQVS1myolzU9xU2VuV6WcBtwZEdNSB1KniDggIlaIiAkUP7ffZFYIKG/e\nfaT8roRi9ugBB8tbd2YwhOnAaZJup5i5MKsfXh9BfqeuJ1DM5HhZcfLDNRHxtbQhjdxAN1UmDqs2\nkt4PfBq4XdLNFP8mD4iIX6WNzIZhd+DHkuYDHmSQm3p905mZmXVdm8jMzBrgYmBmZi4GZmbmYmBm\nZrgYmJkZLgZmZoaLgdmgyvlccpup02wuvs/AzMx8ZmDWQ9JCki6SdLOk2yR9spwdd53y9S9Kuqec\nNfdkSceXz0+X9L+SrpZ0v6SJkmZIulPSab2O/7+SrisXiul3sjCzVFwMzF43CXgsIt4bEWsBr027\nIOmtFDPLrk8xyVnf1tESEbERsCfF/ERHRMS7gLUkrVXuc0BErE+xyMjEcnJCs1ZwMTB73e3AFpIO\nk7Rxn4VA1gdmRsQzEfEK8PM+772w1zGeiIieCcHuoJiNFmB7STcCNwPvKh9mrdBtE9WZNSYi7itb\nQlsD35b0G9648MlgEwf2zML6aq+/92zPWy4Osxfwvoh4VtJ0YIG6YjcbLZ8ZmJXKVtC/ykV3jqJY\nGKTH9RSr7C1ezvG/3WCH6ue5xYDngeckLQtsVVPYZrXwmYHZ69YEjpT0KvAyxZrbRwFExOOSDgWu\nA/4B3A08U76v7yV50ffvEXGbpFso1hF+BLiqqSTMRsKXlppVJGnhiHhB0jzALyjWL8hp9S8bw9wm\nMqvu4HKRl9uBB10ILCc+MzAzM58ZmJmZi4GZmeFiYGZmuBiYmRkuBmZmhouBmZkB/x+N4GqHODke\nUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fae9524e610>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log2(taus), test_powers)\n",
    "plt.xlabel(\"sigma\")\n",
    "plt.ylabel(\"test power\")\n",
    "plt.grid(True)\n",
    "plt.ylim([0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_test_power_new_data(width = 2*0.5**2, num_runs = 100, num_samples = 100, alpha = 0.05):    \n",
    "    test_power = []\n",
    "    for _ in range(num_runs):\n",
    "        #X = norm.rvs(size=num_samples, loc=0.0, scale=1)\n",
    "        #Y = norm.rvs(size=num_samples, loc=1.0, scale=1)\n",
    "                \n",
    "        X,Y = sample_gaussian_vs_laplace(num_samples)\n",
    "        samples_p = RealFeatures(X.reshape(1, len(X)))\n",
    "        samples_q = RealFeatures(Y.reshape(1, len(Y)))\n",
    "\n",
    "        mmd=QuadraticTimeMMD(samples_p, samples_q)            \n",
    "        mmd.set_kernel(GaussianKernel(10, width))\n",
    "        test_power += [mmd.perform_test(alpha)]\n",
    "    return test_power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.31, 0.44, 0.37, 0.059999999999999998, 0.080000000000000002]\n"
     ]
    }
   ],
   "source": [
    "# compute test power for one kernel\n",
    "\n",
    "\n",
    "test_powers = []\n",
    "for tau in taus:\n",
    "    test_power = get_test_power_new_data(width = tau)\n",
    "    test_powers += [np.mean(test_power)]\n",
    "print test_powers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XvcZXPd//HX2/EOlZQfv4hJIkWFkqKacnNPyanj5O5A\nxUiSpB9xK9WvX+mgdLgZOUx30tQtQmkipSKHwRDTDBM5Ralu5BCN8f79sdawXXMd1jXXXnvtva73\n8/HYj9nrtNfn+1hz7c9e389a3yXbREREDGeFpgOIiIj+lSQREREjSpKIiIgRJUlERMSIkiQiImJE\nSRIRETGi2pOEpGmSFkq6QdKho6z3UkmLJb2xY97Nkq6RNE/S5XXHGhERT7RSnR8uaQXg68AOwB3A\nXEln2V44zHqfA3465CMeBabavrvOOCMiYnh1n0lsAyyyfYvtxcBsYLdh1vsgcDpw15D5Il1iERGN\nqfsLeD3gto7p28t5j5H0TGB328dRJIVOBs6XNFfSPrVGGhERy6i1u6mirwCdtYrORLGd7TslrU2R\nLBbYvqi34UVETF51J4k/Aht0TK9fzuv0EmC2JAHPAF4nabHts23fCWD7L5LOpOi+WiZJSMoAVBER\n42R7aO/NMupOEnOBjSVtCNwJTAfe3rmC7Y2Wvpd0CnCO7bMlrQasYPt+SasDOwGfHGlHVRo7iCQd\nZfuopuOoS9o32NK+wVX1x3WtScL2EkkHAOdR1D9Osr1A0oxisU8YuknH+3WAM8uGrAR8x/Z5dcbb\np6Y0HUDNpjQdQM2mNB1AzaY0HUDNpjQdQNNqr0nYngNsOmTezBHWfU/H+z8AL643uoiIGE0uL+1/\ns5oOYCwSq0tstZybz+pmLH1oVtMB1GxW0wHUbFbTATRNbXjokCS3tSbR7yRWA35McQHC+cChNoua\njSoixlL1ezNnEn1O0tSmYxiJxKrAGRT3v6wLXA5cIvEVibWqfUb/tq8b0r7B1vb2VZEkEctFYmXg\ne8D9wN42D9h8Dng+sApwvcSHJVZpMs6ImJh0N8W4SawIfAdYA3ijzT+HWef5wBcoLlo4FDjDZvD/\ns0W0RNXvzSSJGBeJFYCTKYZX2cXmoTHW3xH4EnAvcLDN3PqjjIixpCbREv3UJyoh4BvAc4Ddx0oQ\nADbnA1tSXCVylsSp0uN34fdT++qQ9g22treviiSJqKRMEF8CtgJ2tnmg6rY2S2xOAjYBbgLmSfw/\niafUE21EdEu6m6ISiU8DbwBeazOh53tIrA98hmKolaOAk2wemXCQEVFZahLRNRKHA+8AXm3zly5+\n7lYUZydrA4fYzOnWZ0fE6FKTaImm+0QlDgL2BnboZoIAsLkKVvwkcATwVYmfSmzRzX00renjV7e0\nr/2SJGJEEjOAgygSxJ317OVRbM4CNqe4c/sCiRMk1q1nfxExHuluimFJvIuibjDV5sYe7vdpFGcW\newFfBr5s82Cv9h8xWaS7KZabxFuBzwE79TJBANjcbXMI8DKKUYAXSryzvD8jInosf3h9rtd9ohK7\nAl8FXmezoP79Dd8+mxtt3kLxkKoDgMslXlV3PN3W9j7ttK/9kiTiMRI7AScCb7C5pul4AGwuBl5O\ncRXUf0mcIfHchsOKmDRSkwgAJF4NnE5xJ/XFTcczHIknAR8CDgFOBT5l8z/NRhUxmFKTiMoktgX+\nG3hbvyYIAJt/DBlpdmFGmo2oV5JEn6u7T7S8oe0sYC+bn9e5r+H3P/722dxlsz8wFfhXYL7EG8uh\nQ/pK2/u00772qz1JSJomaaGkGyQdOsp6L5W0WNIbx7ttLB+JF1Dcm/B+m3Objme8bH5nszOwP8Xw\nHr+UeGmzUUW0S601CUkrADcAOwB3AHOB6bYXDrPe+cA/gJNtn1F123L71CTGSWIT4BfAR21Oazqe\niSqfcbEX8Gng58DhNrc2GlREH+uXmsQ2wCLbt9heDMwGdhtmvQ9SFE3vWo5tY5wkng38DDiyDQkC\nMtJsRF3qThLrAbd1TN9eznuMpGcCu9s+Dp7QpzzmtpNBt/tEyxFYLwCOtjm5m5+9fPF0t30299t8\nHHgRxf+X6yVmSKzUzf1U1fY+7bSv/Rr5wxniKxSPt5wQSbOAm8vJe4CrbV9YLpsKMNmnwQuAC+Br\nP4UD51M+TbRf4qth+t1FYf6sk2Dlw6TXv99mTh/Fl+lM9/DvHygu9pjCONRdk9gWOMr2tHL6MMC2\nj+5Y56alb4FnAA8A+1J0PY26bcdnpCYxBomnAxcC/23zqYbD6anyqqddKZ65/QeKYcmvbTaqiGb1\nS01iLrCxpA0lrQJMB87uXMH2RuXr2RR1if1tn11l26hGYk3gPIormT7dcDg9Z+MhI83+LCPNRlRT\na5KwvYRi3J3zgPnAbNsLJM2QtO9wm4y1bZ3x9qOJ9olKrAGcC1wMfMymr26x72Wfr80/bb4KPA/4\nO3CdxBHlndy1aHufdtrXfrXXJGzPATYdMm/mCOu+Z6xtozqJ1YBzKJLsQf2WIJpSPn71EInjKEa7\nvb58+t5pNo82G11Ef8nYTS0lsSrFndR/obibeknDIfUtie2AYyjOrD9i86uGQ4qoXdXvzSSJFpJY\nmWIspkeA6TaPNBxS3yufV/E24LPAVcChNouajSqiPv1SuI4JGm+faHnn8bcpuhL37PcE0S99vjaP\n2nwX2Ay4HLhE4isSa03kc/ulfXVJ+9ovSaJFyl/DJwJPB95s88+GQxo4GWk24onS3dQS5b0A3wC2\nAKbZPNBwSK0g8XyK+ys2objp88xcABBtkJrEJFImiC8C2wM72vy94ZBaR2JHiqfj3UNR3J7bcEgR\nE5KaREtU7BP9FPBaijOIgUoQg9Lna3M+sCXwLeCHEqdKbDDWdoPSvuWV9rVfksSAK6/vfxOwU3n9\nf9SkY6TZTXl8pNnPSDy54dAiapPupgEmcRDwAeBVNnc2Hc9kU46o+xlgJ+ATwMn9fjVZxFKpSbSc\nxAzgYxQJIg/XaVD5CNgvAWtTDB44p+GQIsaUmkRLDNcnKvEu4D+AHQY9QbShz9fmKoqa0OHAsRI/\nldgC2tG+0aR97ZckMWAk3kIx3tBONjc2HU8UypFmz6YYafZHlCPNwgvzZLwYaOluGiASuwDfBP7N\n5pqm44mRSTwN+Cqw2OY9Y60f0WupSbRMeZ3+d4DX21zRdDwxNon/BVwPPNvmnqbjieiUmkRLSJoq\n8WqKBLFH2xJEm/t8be6C788D3tF0LHVp8/GD9reviiSJvrfP8ylGdJ1uc3HT0cR4XXQ2sF95V3zE\nwEl3Ux8rL638CbC3zblNxxPjVyaHBcB7k+Sjn6S7acBJvIDieczvT4IYXOVggCcAM5qOJWJ5JEn0\nIYlNKJ7t/RHQ/zQdT53a3udbtu9bwC4TfTZFP5okx29Sqz1JSJomaaGkGyQdOszyXSVdI2mepCsk\nvbZj2c0dyy6vO9Z+IDEFOB840ua0hsOJLrD5G8VZ4buajiVivGqtSUhaAbgB2AG4A5gLTLe9sGOd\n1Ww/WL7fAjjT9sbl9E3A1rZHHbiuLTUJifWAXwHH2Hyj6XiieyS2p3gg1GZ5HkX0g36pSWwDLLJ9\ni+3FwGxgt84VliaI0hrAXzum1YMY+4LEOsAFwHFJEK10MbAEeFXTgUSMR91fwOsBt3VM317OewJJ\nu0taAJwLHNixyMD5kuZK2qfWSBsk8XTgZ8BpNl984rJ294lOlvaVZw8zaVkBe7Icv8lspaYDALD9\nQ+CHkrYHvk0xXj/AdrbvlLQ2RbJYYPui4T5D0izg5nLyHuBq2xeWy6aW++m7aYk14ce/gT9dCe/9\ndNPxZLq+afC3gU9KW+8GV93bdDyZnlzTpanAFMah7prEtsBRtqeV04cBtn30KNvcCGxj+29D5n8C\nuM/2McNsM5A1CYk1KK5iugL4UPqq20/iFOB3Nl9oOpaY3PqlJjEX2FjShpJWAaYDZ3euIOk5He+3\nArD9N0mrSVqjnL86xYNdrqs53p6ReBJwDjAfOCgJYtKYCewrTY5aWwy+Wv+j2l4CHEDxa3k+MNv2\nAkkzJO1brvYmSddJugo4FnhbOX8d4CJJ84BLgXNsn1dnvL0isSpwBkWNZj+bR0det919opOwfZcB\nDwKv6X003TcJj9+kU3tNwvYcHq8xLJ03s+P954HPD7PdH4AX1x1fr0msTHGV1wMUw20saTik6CEb\nS48VsC9oOp6IsWTsph6SWBE4FXgy8EabfzYcUjRA4qkUF1k8z+bPDYcTk1S/1CSiVPZBnwg8A3hz\nEsTkZXMvcDrkYUTR/5IkeqAcCfTrwMbA7jYPVd+23X2ik7h9M4F9Br2APYmP36Qx0P9BB0GZIL4I\nbA3sbPNAwyFFHygfHnU3sGPTsUSMJjWJmkl8CtgFeK3NqGNQxeQisS8wzeaNTccSk0/V780kiRpJ\nfAx4J/Bqm780HU/0F4knA7cAm9vc0XQ8MbmkcN0wiQ9RFCZ3mEiCaHuf6GRun819wPeA9/YsoC6b\nzMdvskiSqEHZjXAQRYK4s+l4oq/NBN5XXh4d0XfS3dRlEu8CPgNMtbmx6Xii/0lcBnzK5sdNxxKT\nR7qbGiDxFuBzwE5JEDEOM4H9mg4iYjhJEl0isQvwNeB1Ngu697nt7hNN+4CiLvEKiWfVHE7X5fi1\nX5JEF0jsCJwEvMHmmqbjicFS3jtzGvC+pmOJGCo1iQnvm1dRDLGwh83FTcQQg09ic2AOMMXmkabj\nifZLTaIHJF5GkSCmJ0HERNhcR3HPxM5NxxLRKUliOUlsSfEApb1sfl7fftrdJ5r2PcHxDFgBO8ev\n/ZIkloPEC4BzgffbnNt0PNEapwMvkcb3DOKIOqUmMe59sQnwC+CjNqf1Yp8xeUh8GXjQ5oimY4l2\ny9hNteyHpwDXAp+0Obnu/cXkI7EZ8HNgA5vFTccT7ZXCdQ1s/g7s2ssE0fY+0bTvicp7bG4Adqsl\noC7L8Wu/2pOEpGmSFkq6QdKhwyzfVdI1kuZJukLSa6tu24TcBxE9cDzFM7AjGldrd5OkFSh+Fe0A\n3AHMBabbXtixzmq2HyzfbwGcaXvjKtt2fEbfjN0UMVESqwK3Aa+w+X3T8UQ79Ut30zbAItu32F4M\nzGbIafTSBFFaA/hr1W0j2sjmYeBbwD5NxxJRd5JYj+IX0VK3l/OeQNLukhZQXFZ64Hi2bbu294mm\nfSM6AdirPKvoWzl+7bdS0wEA2P4h8ENJrwS+DWw63s+QNAu4uZy8B7ja9oXlsqnlfjKd6YGZBl8L\n7CHpT/0QT6YHe7o0FcZ3H07dNYltgaNsTyunDwNs++hRtrmRoqvpuVW3TU0i2qgcen5/m9c0HUu0\nT7/UJOYCG0vaUNIqwHSKoSweI+k5He+3ArD9tyrbRrTcWcBm0vjPrCO6ZdQkIWlFSb9Y3g+3vQQ4\nADgPmA/Mtr1A0gxJ+5arvUnSdZKuAo6lSAYjbru8sQyqtveJpn0js/kncAqw71jrNiXHr/1GrUnY\nXiLpUUlPtX3v8uzA9hyG1Bhsz+x4/3ng81W3jZhkvglcJnGEzUNNBxOTz5g1CUlnAVsC5wMPLJ1v\n+8ARN+qx1CSizSTmAKfanNp0LNEeVb83q1zddEb5iohmzAQOhiSJ6L1KVzdJehKwge3r6w9p/Np8\nJiFp6tJL2doo7avyGaxMcXn3TjbzuxFXt+T4Da6uXd0kaRfgaopHKyLpxZJylVFEj5SjwZ5MxnOK\nBlSpSVwJvBa40PaW5bzrbG/eg/gqafOZRASAxAbAPOBZNg+OtX7EWLp5n8TiYa5senT5woqI5WFz\nK3AJ8NamY4nJpUqSmC9pT2BFSc+V9DXgNzXHFaW2X6ed9o3LTPqsyynHr/2qJIkPAi8AHga+C9wL\nHFRnUBExrJ8A60u8qOlAYvKoUpN4ju0bexTPcklNIiYLiY8D69rs33QsMdiqfm9WSRK/BNanGEvp\n18CvbF/blSi7JEkiJguJ9Sies76Bzf1NxxODq2uFa9uvBjYDvgasCfxY0v9MPMSoou19omnf+Nj8\nEfgV5RhnTcvxa78x77iWtD3wyvK1JvAjijOKiGjGTOBTwIlNBxLtV6W76RHgSuCzwLm2/9mLwMYj\n3U0xmUisCNwIvMnmyqbjicHUzZrEmsB2wKuAl1LcI3GJ7SO7EWg3JEnEZCNxODDF7t9hxKO/dbMm\ncQ9wE/AH4E7gORQJI3qg7X2iad9yOxl4i8RTavr8SnL82q/K2E03AV8C1gKOAzYti9kR0RCbPwEX\nAHs2HUu0W5XuphVs9/UwHOluislIYkfgC8CWNvU9rD5aqZtjNz1T0pmS7ipfP5C0fhdijIiJuQBY\nA9im6UCivaokiVOAs4Fnlq9zynnRA23vE037lp/No8AJNDieU45f+1VJEmvbPsX2I+VrFrB21R1I\nmiZpoaQbJB06zPI9JV1Tvi6S9MKOZTeX8+dJurzqPiMmkVnAHhJrNh1ItFOVmsQFFGcO3y1nvR3Y\n2/YOY364tAJwA7ADcAfF0B7TbS/sWGdbYIHteyVNA46yvW257CZga9t3j7Gf1CRi0pKYDVxk8/Wm\nY4nB0c2axHsoxrD/U/l6M7B3xTi2ARbZvsX2YmA2sFvnCrYv7XhexaXAeh2LVTHGiMlsJrCfRH4o\nRddVuU/iFtu72l67fO1u+9aKn78ecFvH9O08MQkM9T6K4ZAf2z1wvqS5kvapuM9WaXufaNrXFRcC\nKwOv6MG+niDHr/2qjN20EXAssC3Fl/YlwIdt39TNQCS9huIMZfuO2dvZvlPS2hTJYoHti0bYfhbF\nw+IB7gGuXvoA86UHOtOZbuM06NVw7M/gwBnAxU3Hk+n+nC5NBaYwDlVqEpcC3+DxmsR04IO2Xzbm\nhxf1hqNsTyunDyti99FD1nsh8ANgmkd4doWkTwD32T5mmGWpScSkJvF0ivGcNrLJKM0xpm7WJFaz\n/e2Oq5tOBf6lYhxzgY0lbShpFYoEc/aQQDegSBDv7EwQklaTtEb5fnVgJ+C6ivuNmFRs/kYxQvO7\nm44l2qVKkviJpMMkTSm/7P8PcK6ktSStNdqGtpcABwDnAfOB2bYXSJohaenAZEdSDPnxn0MudV0H\nuEjSPIqC9jm2z1uONg60tveJpn1ddTwwo5cF7By/9qvS3fSHURbb9kbdDWn82tzdJGnq433P7ZP2\ndXNfiOJse3+bX/Zmnzl+g6rq9+aYSWIQtDlJRIyHxIHAtnYG/ovRJUlETEIST6MY2n8Tm780HU/0\nr24WrqNBbe8TTfu6y+Zu4Cxgr17sL8ev/ZIkItrneGBfKX/fMXGVxm7ykHGahpvXpHQ3RTyuLGBf\nDRxsc0HT8UR/qvq9OeId15L+BVgNeIakp8Fjl9U9hdGH1oiIBtlYYibFEOJJEjEho52OzgCuBJ5X\n/rv0dRZktMleaXufaNpXm+8AO0qsW+dOcvzab8QkYftY288GDrG9ke1nl68X2U6SiOhjNvcCp1N9\nxOaIYVWpSbwFmGP7Pkn/AWwF/F/bV/UiwCpSk4hYlsRLgO8DG5dPsYt4TDcvgT2yTBDbA/8KnAQc\nN9EAI6JeNlcAdwM7Nh1LDK4qSWJJ+e/OwAm2fwysUl9I0antfaJpX+2WFrBr0Qftq1Xb21dFlSTx\nR0kzgbdRDOy3asXtIqJ53wWmSjyz6UBiMFWpSawGTAOutb1I0v8GtuinEVlTk4gYmcRxwB02n246\nlugfXatJ2H4QuIvHnxj3CLBoYuFFRA/NBN4nsWLTgcTgGTNJlE+EOxT4WDlrZeDUOoOKx7W9TzTt\nq5/N1cCfKHoEuqof2lentreviiq1hT2AXYEHAGzfATy5zqAioutmAvs1HUQMnio1icttbyPpKttb\nlY8SvcT2C3sT4thSk4gYncTqwK3Ai21uazqeaF4375P4fnl105qS9gF+Bpw40QAjondsHgBOA97X\ndCwxWKoUrr9IcXv/D4BNgY/b/mrdgUWh7X2iaV9PzQTeK408sOd49Vn7uq7t7auiSuH6aNvn2/6o\n7UNsny/p6Ko7kDRN0kJJN0g6dJjle0q6pnxdJOmFVbeNiOpsrgNuAd7QdCwxOKrUJK6yvdWQeb+t\nUpOQtAJwA7ADcAcwF5hue2HHOtsCC2zfK2kacJTtbats2/EZqUlEVCDxTmBPm9c1HUs0a8I1CUnv\nl3QtsKmk33a8/gD8tmIc2wCLbN9iezEwG9itcwXbl9q+t5y8lMefVTHmthExbqcDL5GY0nQgMRhG\n6246DdgFOLv8d+lra9vvqPj568ETrqS4ndEfWPQ+4CfLuW0rtb1PNO3rLZt/UNzntE83Pq/f2tdt\nbW9fFSMWsMpf9/cCb+9FIJJeQzH2/fZjrTvC9rOAm8vJe4CrbV9YLpsKkOlMZxrgTfNgvy9LOx5l\ns7jpeDLdm+nSVBjfWeSYNYmJKOsNR9meVk4fBtj20UPWeyHF1VPTbN84nm3LZalJRIyDxC+Br9mc\n3nQs0YwJ1yS6ZC6wsaQNJa0CTKfovnqMpA0oEsQ7lyaIqttGxHI7nhqHEI/2qDVJ2F4CHACcB8wH\nZtteIGmGpH3L1Y4E1gL+U9I8SZePtm2d8fajtveJpn2NOQN4kcTGE/mQPm5fV7S9fVV07aaakdie\nQ3ETXue8mR3v92GEItpw20bExNk8LPEtir+93IMUI6q1JtErqUlEjJ/Ec4GLgWfZPNx0PNFb/VKT\niIg+ZbOI4p6nPZqOJfpXkkSfa3ufaNrXuAk9A3sA2jchbW9fFUkSEZPbWcBmUmp/MbzUJCImOYnP\nAqvYfKTpWKJ3qn5vJklETHISGwGXURSwH2o6nuiNFK5bou19omlf82xuAq4E3jzebQehfRPR9vZV\nkSQRETDBAna0V7qbIgKJlSkGyNzJZn7D4UQPpLspIiqzWQycTM4mYogkiT7X9j7RtK+vfBP4d4nV\nqm4wYO0bt7a3r4okiYgAwOZW4BLgrU3HEv0jNYmIeIzELsDhNi9vOpaoV2oSEbE8fgKsL/GipgOJ\n/pAk0efa3iea9vUXm0coahOVCtiD1r7xanv7qkiSiIihTgKmS6zRdCDRvNQkImIZEj8EfmRzYtOx\nRD1Sk4iIicgd2AEkSfS9tveJpn196zxgbYmtR1tpgNtXSdvbV0XtSULSNEkLJd0gaZln6UraVNJv\nJD0k6eAhy26WdI2keZIurzvWiCjYLAFOIGcTk16tNQlJKwA3ADsAdwBzgem2F3as8wxgQ2B34G7b\nx3QsuwnY2vbdY+wnNYmILpNYF1gAbGjz96bjie7ql5rENsAi27fYXgzMBnbrXMH2X21fCTwyzPbq\nQYwRMQybPwEXAHs2HUs0p+4v4PWA2zqmby/nVWXgfElzJe3T1cgGRNv7RNO+vjcT2E9i2F+cLWjf\nqNrevipWajqAMWxn+05Ja1MkiwW2LxpuRUmzKIY6BrgHuNr2heWyqQCZznSmxz19AcxZG07fD048\nrg/iyfRyTpemAlMYh7prEtsCR9meVk4fBtj20cOs+wngvs6aRNXlqUlE1Efi/wDPs3lP07FE9/RL\nTWIusLGkDSWtAkwHzh5l/ccClrSapDXK96sDOwHX1RlsRAxrFrCHxJpNBxK9V2uSsL0EOIDimuv5\nwGzbCyTNkLQvgKR1JN0GfBg4QtKtZXJYB7hI0jzgUuAc2+fVGW8/anufaNrX/2zuAn4KvGPosja0\nbzRtb18VtdckbM8BNh0yb2bH+z8Dzxpm0/uBF9cbXURUNBP4msQ3bAZ/LJ+oLGM3RcSYyqubFgLv\nsbm46Xhi4vqlJhERLVCePWQ8p0koSaLPtb1PNO0bKN8CdpVYa+mMlrVvGW1vXxVJEhFRic3fgB8B\n7246luid1CQiojKJ7YETgc1SwB5sqUlERB0uBpYAr2o6kOiNJIk+1/Y+0bRvsAwtYLetfUO1vX1V\nJElExHh9G3i9xNpNBxL1S00iIsZNYhYw3+YLTccSyyc1iYio0/HAvlK+Q9ouB7jPtb1PNO0bWJcB\nD8KHD2o6kDq1+PhVliQREeP2eAH75bs0HUvUKzWJiFguEk+leNDXZuWjTmOApCYREbWyuRc4Hdi7\n6ViiPkkSfa7tfaJp36Db+wpgn7YWsNt//MbWygMbEb0y63rgbmDHpiOJeqQmERETIrEv8DqbPZqO\nJaqr+r2ZJBEREyLxZOAWYHObO5qOJ6pJ4bol2t4nmvYNNklTbe4Dvge8t+l4uq3tx6+K2pOEpGmS\nFkq6QdKhwyzfVNJvJD0k6eDxbBsRfWMm8D6JFZsOJLqr1u4mSSsANwA7AHcAc4Hpthd2rPMMYENg\nd+Bu28dU3bbjM9LdFNEwicuAT9n8uOlYYmz90t20DbDI9i22FwOzgd06V7D9V9tXAo+Md9uI6Csz\ngf2aDiK6q+4ksR5wW8f07eW8urdtjbb3iaZ9g21I+74HvELiWQ2F03VtP35VrNR0AN0iaRbFEAEA\n9wBX276wXDYVINOZznS90xKnwQmflmbM6od4Mr1MopsKTGEc6q5JbAscZXtaOX0YYNtHD7PuJ4D7\nOmoS49k2NYmIPiCxOTAHmGIv04UcfaRfahJzgY0lbShpFWA6cPYo63cGPN5tI6JhNtdR3DPxhqZj\nie6oNUnYXgIcAJwHzAdm214gaYakfQEkrSPpNuDDwBGSbpW0xkjb1hlvP2p7n2jaN9hGaN/xlM/A\nHnRtP35V1F6TsD0H2HTIvJkd7/8Mwxe6hts2Ivre6cAxElPsx+qEMaAyLEdEdJ3El4EHbY5oOpYY\nXtXvzSSJiOg6ic2AnwMb2CxuOp5YVr8UrmOC2t4nmvYNtpHaZ7OAYsSEgb4Btu3Hr4okiYioS2sK\n2JNZupsiohYSq1KMmvAKm983HU88UbqbIqJRNg8D3wL2aTqWfiWxosRGEjtL7N90PMNJkuhzbe8T\nTfsGW4X2nQDsXZ5VDJxuHT+JVSU2l3iLxMclvitxNXAfcCFwIPDcfnxWeGvGboqI/mOzSOK3wB4U\nIzm3msRTgM3K1/M63m9AMbbcgvJ1LvAlYKHN/Y0EW1FqEhFRK4m3APvbvKbpWLpBQsA6PJ4AOl9r\nAtfzeDJVXXJsAAAHMUlEQVRY+vq9zT8bCXgEuU8iIvqCxCrArcBUm2UeGtavyqfsbcjwyeBRlk0E\nC4BbbR5tJOBxSpJoCUlTlw7520Zp32Cr2j6JzwKr2hw81rq9VtZLnsuyiWATOO8+2OlqhiQDm780\nFW+3VP3eTE0iInrhm8BlEofbPNREAEPqBZ2vZzFCvQD+7SVtTvJV5EwiInpCYg5wqs2pNe6jFfWC\nXkh3U0T0FYk9gINtXtmFzxqtXrCEJyaBhQxYvaAXkiRaIn3agy3t61yXlSm6dXaymV9xm1HqBfyV\nYYrH3awXtPn4pSYREX3FZrHEyRTjOR3YuWx56gX9fn9BW+RMIiJ6RmIDYB5wJKkXNCrdTRHRlyQ+\nTVFcXlorSL2gAX0zwJ+kaZIWSrpB0qEjrPNVSYskXS1py475N0u6RtI8SZfXHWs/ytg/gy3tW5bN\nkTb72hxj8xObm/s1QbT9+FVRa5KQtALwdeDfgBcAb5f0vCHrvA54ju3nUvRVHtex+FFgqu0tbW9T\nZ6x97MVNB1CztG+wpX0tV/eZxDbAItu32F5MMcDX0CdV7Qb8F4Dty4CnSlqnXKYexNjv1mw6gJql\nfYMt7Wu5ur+A16N46MhSt5fzRlvnjx3rGDhf0lxJGZM+IqLH+v0S2O1s3ylpbYpkscD2RU0H1WNT\nmg6gZlOaDqBmU5oOoGZTmg6gZlOaDqBpdSeJP1KMo77U+uW8oes8a7h1bN9Z/vsXSWdSdF8NmyQk\nDf5lWiOQ9O6mY6hT2jfY0r52qztJzAU2lrQhcCcwHXj7kHXOBj4AfE/StsA9tv8saTVgBdv3S1od\n2An45HA7yeWvERH1qDVJ2F4i6QDgPIr6x0m2F0iaUSz2CbbPlfR6Sb8HHgD2LjdfBzizPENYCfiO\n7fPqjDciIp6oFTfTRUREPVpzeamkD0paIOlaSZ9rOp46SPqIpEclrdV0LN0k6fPlsbta0g8kPaXp\nmCaqyk2kg0rS+pJ+Lml++fd24NhbDR5JK0i6StLZTcfSbZKeKum/y7+7+ZJeNtK6rUgS5V2RuwBb\n2N4C+GKzEXWfpPWBHYFbmo6lBucBL7D9YmAR8LGG45mQKjeRDrhHgINtvwB4OfCBlrVvqQ8Bv2s6\niJocC5xrezPgRRRDowyrFUkCeD/wOduPANj+a8Px1OHLwEebDqIOtn9me+mwDJdSXOE2yKrcRDqw\nbP/J9tXl+/spvmCG3v800MofZa8HTmw6lm4rz9RfafsUANuP2P77SOu3JUlsArxK0qWSfiHpJU0H\n1E2SdgVus31t07H0wHuAnzQdxARVuYm0FSRNoRi64rJmI+m6pT/K2li0fTbwV0mnlN1pJ0h60kgr\n9/vNdI+RdD7FFU+PzaI4gP9B0Y6n2d5W0kuB7wMb9T7K5TdG+w6n6GrqXDZQRmnfEbbPKdc5Alhs\n+7QGQoxxkrQGcDrwofKMohUk7Qz82fbVZVf2wP29jWElYCvgA7avkPQV4DDgEyOtPBBs7zjSMkn7\nAWeU680ti7tPt/23ngU4QSO1T9LmFHd9XiNJFF0xV0raxvZdPQxxQkY7fgCS9qI4vX9tTwKqV5Wb\nSAeapJUoEsS3bZ/VdDxdth2wq6TXA08Cnizpv2y/q+G4uuV2ip6JK8rp04ERL65oS3fTDym/XCRt\nAqw8SAliNLavs72u7Y1sP5viAG85SAliLJKmUZza72r74abj6YLHbiKVtArFTaRtu0LmZOB3to9t\nOpBus3247Q1sb0Rx7H7eogSB7T8Dt5XflQA7MEqBfmDOJMZwCnCypGuBh4HWHNBhmPad/n4NWIVi\nfC6AS23v32xIy2+km0gbDqtrJG0H/DtwraR5FP8nD7c9p9nIYhwOBL4jaWXgJh6/iXkZuZkuIiJG\n1JbupoiIqEGSREREjChJIiIiRpQkERERI0qSiIiIESVJRETEiJIkIrpE0o+GG+Zc0ickHVy+f7ek\ndTuW/aFtQ79HuyRJRHSJ7TeMNppmaS+eONhfblSKvpYkEVGRpEPKO6mR9GVJF5TvXyPp1M6zAklH\nSLpe0q+ATct5bwJeApxajr75LxR3zx8o6UpJ13QMlRDRF5IkIqr7NfDK8v3WwOqSVizn/ZLyrEDS\nVsBbgRcCOwMvBbD9A+AKYE/bW9l+qPysu2xvDRxPS58ZEoMrSSKiuiuBrSU9mWKMsEsoEsArKRLI\nUq8EzrT9sO37WHZwv6Fjb53Z8fkbdj3qiAloywB/EbWz/YikmynqChcDvwVeAzzH9sJyKPflsXTk\n2yXkbzL6TM4kIsbn18AhwK+Ai4D9gKuGrPMrYHdJq5ZnHbt0LLsPWOYKqIh+lSQRMT6/BtYFLimf\n6fEPHu9qMoDtecD3KM40fgxc3rH9LOD4jsJ1rm6KvpahwiMiYkQ5k4iIiBElSURExIiSJCIiYkRJ\nEhERMaIkiYiIGFGSREREjChJIiIiRpQkERERI/r/VbqCeNX/kj8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fae8f18ad50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log2(taus), test_powers)\n",
    "plt.xlabel(\"width\")\n",
    "plt.ylabel(\"test power\")\n",
    "plt.grid(True)\n",
    "#plt.ylim([0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# TODO\n",
    "\n",
    "try:\n",
    " * max mmd\n",
    " * max test power\n",
    " * weighted version and plain (single kernel)\n",
    " * median\n",
    " * xvalidation (above)"
   ]
  },
  {
   "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
}