{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Intuition for the Maximum Mean Discrepancy two-sample test\n", "\n", "_Thomas Viehmann_\n", "\n", "This note sketches the intuition behind [A. Gretton et al.: A Kernel Two-Sample Test. JMLR 2012](http://www.gatsby.ucl.ac.uk/~gretton/mmd/mmd.htm).\n", "\n", "Given a (high-dimensional) space $\\mathbb{R}^d$ and iid samples $X_i \\in \\mathbb{R}^d, i=1,...N_X$ sampled from $X ~ P_X$ and $Y_i \\in \\mathbb{R}^d, i=1,...N_Y$ sampled from $Y ~ Y_X$, we may wish to make a statistical test of the hypothesis that $P_X$ is different from $P_Y$.\n", "\n", "For simplicity, we stay in the classical, frequentist scheme.\n", "\n", "So let us import the world:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import torch\n", "%matplotlib inline\n", "from matplotlib import pyplot\n", "import tqdm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then we can sample a few points (to facilitate plotting, we'll operate in 1-d, so not that high dimensional, but hey.)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyUklEQVR4nO3dd3wc13Xo8d9Z9EYQRG8EO0GCXWBXI2k1UoqbbEt2bDmRoyfbcpyXOI7t5DmK0/3ec5GdWG6KSxTpxZaLLFGWZIoyRZGiCPYC9oZG9N6Bve+P2YVAECAWwM7OAHu+nw8+C+zOzD1YLg8vz9wixhiUUkpNfh6nA1BKKRUcmtCVUmqK0ISulFJThCZ0pZSaIjShK6XUFBHpVMNpaWlm1qxZTjWvlFKT0oEDB+qMMenDveZYQp81axYlJSVONa+UUpOSiFwe6TUtuSil1BShCV0ppaaIURO6iOSLyE4RKRWREyLy2WGOuV1EmkXksO/ry/aEq5RSaiSB1ND7gL8wxhwUkSTggIi8aow5OeS4N4wx9wY/RKWUUoEYtYdujKkyxhz0fd8KlAK5dgemlFJqbMZUQxeRWcBKYN8wL68XkSMi8pKIFI1w/iMiUiIiJbW1tWOPViml1IgCTugikgg8B/yZMaZlyMsHgQJjzHLgW8CvhruGMeZ7xphiY0xxevqwwyiVUkqNU0Dj0EUkCiuZP22M+cXQ1wcneGPMdhH5dxFJM8bUBS9UNVkYY6hr6+FKQwfljR0UpCawIn+602GN7vIeOL/znZ8zCqHofSDiXExKjcGoCV1EBPghUGqM+doIx2QB1cYYIyJrsHr+9UGNVE0KzR29PPQfb3O4rOma5z+8diZfuKeQabFRzgR2I8bAvu/Cy18E4wUE8O0TcPEN2Pq/IcKFcSs1RCA99I3AR4FjInLY99yXgJkAxpgngfuBT4pIH9AJPGB054yw09rVy8f+421KK1v4q7sLWZiVSM70OJ47UM4Pd1/ktdIa/vX+Zdy2wEXltv5e2P45OPAjKLwX3vtdiEkErxde+3vY/TVoOA8f+DHEz3A6WqVuSJzKu8XFxUan/k8d7d19PPSU1TP/94+s4s6irGteP1zWxOd/foRLdR0898kNLM1LdijSQYyBnz0EJ38NN/85bP5f4BlyW+nIs/D8Z2DGHHjkdYiKcyRUpfxE5IAxpni413SmqJowYwyP/ddBDl5p5JsPrLwumQOsyJ/Os4+sJzUxmk/91wGaO3odiHSI09utZL75b+Bdf3t9MgdY/gA88AzUnoLdXw99jEqNgSZ0NWH7Ljaw83Qtf3V3IduWZY943IyEaL794VVUNXXxuZ8fwdGqXE87vPQFSF8EG//sxsfOfxcs/QDs/gbUnw9FdEqNiyZ0NWHffu0caYkxPLRh1qjH3lSQwpe2LuLVk9V8b9cF+4Mbya7/A81X4N6vBXbD885/gMgY2P6XVqlGKRfShK4m5HBZE7vP1fEnt8wmNioioHP+aOMs7lmSxf995QxlDR02RziM2jOw51uw/EEo2BDYOUlZsOlLcH4HlD5vb3xKjZMmdDUh337tHMlxUXxkXUHA54gIX75vMSLw9VfP2BjdCF7+IkTHwx1fGdt5q/8EMpfCb78Ifd32xKbUBGhCV+NWWtXC70qr+eONs0mMGdteKdnJcXx84yx+ebiC0qqhE49tVHsazv0ONn4WEjPGdm5EJNzxOLRUWDdTlXIZTehq3P5t5zkSYyL5eAC18+F86rZ5JMVE8tXfngpuYDey/4cQEQ2rHhrf+XM2w4y5sP8HwY1LqSDQhK7Gpamjh5eOX+XBNfkkx49vFmVyfBSfvH0eO0/Xsu9CCCYWd7fBkWeg6L2QkDa+a3g8sPphKNsHVUeDG59SE6QJXY3LztM19HsN25blTOg6f7RxFlnTYvmX356yfxjjsf+G7hZY/YmJXWfFhyEyTnvpynU0oatxefVkNRlJMSzLndiMz9ioCD69eR6HrjSx/1JjkKIbhjFWuSVrKeStnti14lJg6fvh2M+gsyko4SkVDJrQ1Zh19/Xz+9O1bFmUiccz8ZUI71+VR3JcFD/aczEI0Y2gbB9UH7d658FYPXH1J6C3w1oaQCmX0ISuxmzv+Xrae/q5c3FmUK4XFx3BA2vyeflENRVNnUG55nX2/wBikq0Zn8GQsxJyi63r6kQj5RKa0NWYvXqymvjoCNbPTQ3aNT+6rgBjDD/dezlo1xzQ3Qonn4flH4LohOBdt/iPof4slOsic8odNKGrMfF6Db8rrebW+ekBzwwNRF5KPHcuzuLZ/Vfo7OkP2nUBOPsK9Hdbm1UEU+E28ETqzFHlGprQ1Zgcr2ymuqWbO4JUbhns4xtn0dTRy68PVwT3wqW/gYQMyF8T3OvGTYfZt1nX17KLcgFN6GpMXj1ZjUdgU+EYZ1kGYO3sGSzKnsZ/vHkpeEMYezvhzCuw6F7wBO9/FAMW/wE0XoTqE8G/tlJjpAldjcmrJ6spnjWDGQnRQb+2iPDQ+gJOV7dy8EpTcC56fif0tsOi+4JzvaEWbgPxaNlFuYImdBWwmtYuTl1tZYsNvXO/bcuyiY3y8NzB8uBcsPR5iJ0Os24JzvWGSkyHmRussotSDtOErgJ22NdrvqkgxbY2kmKjuGdJNr85UklX7wRvjvb3WrsSLdxq7ybPi+6DmpNQd86+NpQKgCZ0FbDDZU1EeoQlE5wdOpr3r8qjtauPV09WT+xCF3dBV7N95Ra/Rfdaj1p2UQ7ThK4CduhKE4tzpgV1uOJw1s9NJSc5duJll9LfQFQCzN0cnMBGkpwHuTdp2UU5ThO6Cki/13C0vIkV+dNtbyvCI7x3VS67ztRS09I1vot4vXDqRVhwJ0TFBjfA4Sy6DyoPQnOQav9KjYMmdBWQszWttPf0s3Lm9JC0975VeXgN/PLQOMekXz0K7TWw4O7gBjaSBfdYj+d2hKY9pYahCV0F5JDvhujKfPtuiA42Nz2RVTOn89zB8vGNSb+w03qcsym4gY0kfSEkZb/TrlIO0ISuAnLoSiMp8VEUpMaHrM3335THmeo2TlSOY4u68zshowiSgj+jdVgi1j8eF35vlXuUcoAmdBWQw2VW/VyCsfRsgLYuySbCI7x4rGpsJ/Z2wpW3YG6Ieud+czdBZwNcPRLadpXy0YSuRtXa1cvZmjZWzgxNucUvJSGajfPSePFo1djKLpf3WItxharc4jfnduvxvJZdlDM0oatRHS1vxhhCMsJlqHuXZnOloYPjFWMou1zYaW0EXbDBvsCGk5gBmUu0jq4cowldjerQFWtruOUOJPQ7izKJ9AgvHK0M/KTzr0P+WogOXb1/wJzbrXJPT0fo21ZhTxO6GtWhK03My0gkOc7G6fMjmB4fzS3z03gh0LJLWw1UHwt9/dxv7ibo74Ere5xpX4U1TejqhowxAzdEnbJtWQ4VTZ0cKW8e/eALv7ceQ10/95u5wSr3aB1dOWDUhC4i+SKyU0RKReSEiHx2mGNERJ4QkXMiclREVtkTrgq1iqZO6tt7HCm3+N2xOJOoCOHFQMouF3ZCXApkL7c/sOFEx8PMdXDhdWfaV2EtkB56H/AXxphFwDrg0yKyeMgx9wDzfV+PAN8JapTKMWdr2gAozEpyLIbkuChunZ/Oi0er8HpvUHYxxuoZz77Nns0sAjVnE1Qfh9YJLi6m1BiNmtCNMVXGmIO+71uBUiB3yGHvBn5iLG8B00UkO+jRqpA7V20l9HnpiY7GsW1ZNpXNXRwubxr5oPrz0FoJc24LWVzD8g9fvPSGo2Go8DOmGrqIzAJWAvuGvJQLlA36uZzrkz4i8oiIlIhISW1t7RhDVU44W9NKWmIMKTbsUDQWWxZZo11ePnF15IP8NyILNoYmqJFkLYPoRLiy19k4VNgJOKGLSCLwHPBnxpihg4KHmz543f+NjTHfM8YUG2OK09PTxxapcsTZmjbmZzjbOwer7LJ+biqvnKgeebTL5b0QnwppC0Ib3FARkZC32opHqRAKKKGLSBRWMn/aGPOLYQ4pB/IH/ZwHjGHgsHIjYwznqtuYn+l8Qge4syiLi3XtnPPV9a9zZS/MXG+tq+K0gg3WLkadjU5HosJIIKNcBPghUGqM+doIhz0PfMw32mUd0GyMGeMCHMptqlu6ae3uc0UPHeDOxdZCW8OWXVqvQuNFK6G7wcz1gIErQ6uTStknkB76RuCjwGYROez72ioij4rIo75jtgMXgHPA94FP2ROuCqWzNa0AzMtwboTLYJnTYlmRP51Xhtua7rK/fu6ShJ5XDJ4onWCkQipytAOMMbsZvkY++BgDfDpYQSl3OOsb4eKWkgvAXUVZ/OtvT1HZ1EnO9Lh3Xriy19puLsuh8edDRcVBzkqto6uQ0pmiakRna9pIiY8i1eERLoPdWWSVXV4ZWna5vBfyV1s3JN2iYD1UHrKW81UqBDShqxGdq2llfkZSSNdAH83c9ETmZSReW3bpbLIm8swM8eqKo5m5Aby9UF7idCQqTGhCV8MyxnCmuo15Liq3+N1VlMm+iw00tvdYT5S9DRj31M/9Zq61HnU8ugoRTehqWHVtPTR39rpmhMtgdy7Oot9reO1UjfXElT3giYTcYmcDGyouBTIWv3PDVimbaUJXw/KPcJnvkhEugy3NTSYjKYYdp3xll8t7IXuFM+ufj2bmeijfD/19TkeiwoAmdDUs/+QdN41w8fN4hC2LMth1po6erg6oPOi+cotfwQboabPWaFfKZprQ1bDOVreRFBtJRlKM06EMa0thJm3dfZw8+Ia1oYRbJhQNNXOd9Vj2trNxqLCgCV0N62xNK/MzEl01wmWwjfPSiIn0UH3Ct6Jh3mpnAxpJch4k5VhlF6VspgldDetcTZsr6+d+cdER3DI/jairBzDTZ1obNLtVXrEmdBUSmtDVdRrae6hr63Fl/XywLYsyWdh3mta0FU6HcmN5q6HxErTpktHKXprQ1XUu1rUDMCc9weFIbuyO3H5ypZ4jZr7TodyYvxxUoROMlL00oavrlDd2AJCf4sJhgIOkNVsjR37TcN1eKu6SvRwkQssuynaa0NV1yhqshJ7n8oRO+X76JYpfX02ltrXb6WhGFh0PWUs0oSvbaUJX1ylr6CQtMYa4aAc3Wg5EeQnd6UvpNlHs9M8adau81VBxCLz9TkeipjBN6Oo6ZY0d5M+IG/1AJ/X3QuUh4mavJTs59p1lANwqbzX0tELtaacjUVOYJnR1nbLGDtfXz6k+AX2dSP5qbl+Ywe5zdfT0eZ2OamT+G6NadlE20oSurtHX76Wyqcv9PXR/YswtZkthBm3dfey/1OBsTDcyY461WJcmdGUjTejqGlXNXfR7zSS4IVoCCRkwfSYb5qUSHelhR6mLyy4iVi9d10ZXNtKErq5R3mjtruP6kkv5fitBihAfHcn6OansPO3ihA5WvLWnoKvZ6UjUFKUJXV2jzD8G3c0ll44GaDhvTan32bIog4t17VyobXMwsFHkFQMGKg46HYmaojShq2uUN3TgEa7dgNltKg5Yj4MS+qaF1lourh7tknuT9agzRpVNNKGra5Q1dpKdHEdUhIs/GhUHAbE2tfDJnxHPgsxEd5ddYpMhdZ41Hl0pG7j4b61yQllDB3kpLu6dg7WhRdp8iJ12zdObCjN4+2IDrV29DgUWgJxVVvxK2UATurqGNanIxTdEja8GnbPqupc2L8ygt9+w+2ydA4EFKHcVtFZBS5XTkagpSBO6GtDV2091S7e7e+gtFdBeYyXGIW4qSGFabKS76+j+f4i0l65soAldDahomgRDFv0jRIbpoUdGeLh1QTqvn6nF6zUhDixAWUutlRd1pIuygSZ0NWBgDLqbSy6VB8ETaSXGYWwuzKC2tZsTlS0hDixA0fGQsVh76MoWmtDVAP+yua4eg15xEDKLICp22JdvW5COiNuHL66EykPW/QClgkgTuhpQ1thBdISHzKThk6XjvF6oPDxsucUvNTGG5XnTec3NwxdzVkFnIzRedDoSNcVoQlcDyhs6yU2Jw+MRp0MZXsMF6G4e9oboYJsLMzha3kRdm0s3vfDHr3V0FWSa0NWAskaXj0GvHPmG6GCbCzMwBn5/2qWbMmcshshYq+yiVBCNmtBF5CkRqRGR4yO8fruINIvIYd/Xl4MfpgoFa1KRi2+IVhyEyDhIL7zhYUU508hIinFv2SUiyrqp61/CQKkgCaSH/iPg7lGOecMYs8L39ZWJh6VCra27j8aOXnffEK08aG24HBF5w8NEhE0LM9h1ppbefpduepGzCqqOQH+f05GoKWTUhG6M2QW4eOcAFQwDI1zc2kPv74Oqo6PWz/02FabT2tXHgcuNNgc2TrmroLcD6nRLOhU8waqhrxeRIyLykogUjXSQiDwiIiUiUlJb69L6Zpiq9E0qynVrDb22FPo6R62f+908P52oCHHv5tE5emNUBV8wEvpBoMAYsxz4FvCrkQ40xnzPGFNsjClOT08PQtMqWKpbrBEhWdNcOmTRfwMxZ2VAhyfGRLJm9gz3jkdPnQfRSXpjVAXVhBO6MabFGNPm+347ECUiaROOTIXU1ZYuRCA9KcbpUIZXeQhikq29OQO0aWEGZ2vaBspJruLxQM4KTegqqCac0EUkS0TE9/0a3zXrJ3pdFVo1LV2kJsS4dx30ysOQs9xKhAHasigTwL1rpOesgOoT0NfjdCRqighk2OIzwF5goYiUi8jDIvKoiDzqO+R+4LiIHAGeAB4wRuc0TzZXW7rISnZp77yvB6qPX7OhRSBmpyUwOy3BvWWX7BXQ323dH1AqCG48/gswxjw4yuvfBr4dtIiUI6pbuslJdmn9vOYk9PcEXD8fbNPCDP5z32U6evqIjx714x5a/t+n8pA1HFOpCXLp/69VqNW0dJHp1oQ+xhuig20uzKCnz8uecy6sAs6YY90X0Dq6ChJN6Iruvn7q23vcuyhX5SGInQ4ps8Z86prZM0iIjnDnrFERvTGqgkoTuqK21Tdk0a019KrDVu9cxr5oWHSkh1vmp7PzVA2uvLWTsxKqT0KfSxcSU5OKJnRFdUsXABluHIPe22UlvHGUW/w2F2ZQ1dzFqautQQwsSHJWgrfXGu2i1ARpQlfunlRUc8JKeDkrxn2J2wutSWyuHO3i/7207KKCQBO64mqz1UPPdGNCn8ANUb+MpFiW5SWzo7Q6SEEF0fQCiEvRhK6CQhO6orq1i+gIDynxUU6Hcr3KQxCfCsn5E7rM5sIMDpU1Ue+2TS9ErH+sKg87HYmaAjShK6qbu8iYFoOM46aj7SqPjPuG6GDvWpSJMbDTjZte5Ky0Jhf1djodiZrkNKErqlu63Vk/7+20JhVNoNziV5QzjcxpMe4su+SsBG+f3hhVE6YJXVHd0uXO+vnV42D6xzzlfzgiwubCTHadqaWnz2WbXvh/P62jqwnShK7cm9AHboiuCMrlthRm0N7Tz76LLps1mpwH8Wm6NrqaME3oYa61q5f2nn4yp7lwUlHlIUjIgGm5QbncxnlpxER62FHqsuGL/hujVYedjkRNcprQw9zAGHQ3ruNSeSgoN0T94qIj2DgvjR2nqt03azRnJdSegp52pyNRk5gm9DA3MEvUbeu4dLdZ+20G4YboYFsWZVDW0MnZmragXnfCclaC8cLVY05HoiYxTehhzp/QXddDv3rMSnBBTuibCzMA3Fd2GbyUrlLjpAk9zF1t8c8SdVkNPcg3RP2yk+MoypnG79w2fHFaNiRmaUJXE6IJPczVtHSTFBvpvs0fKg9BUg4kZQX90ncszuTglcaBVSZdI2elJnQ1IZrQw9zVZhcPWQxyucVvYNao2xbrylkJdWehq8XpSNQkpQk9zFW3drlvlmhXC9SftS2hF+VMI3d6HK+cdFnZJWclYODqUacjUZOUJvQw51/HxVWqjliPNiV0EeFdizLYfa6Wzp5+W9oYF11KV02QJvQw5vUaalpduI6LTTdEB7tjcRZdvV7eOOuixboSM2BaniZ0NW6a0MNYfXsPfV7jviGLlYcgeSYkpNnWxNo5M0iKjXTfaBfdY1RNgCb0MObaSUWVh2ztnQNERXjYtDCDHaU19HtdNGs0ZyU0XIDORqcjUZOQJvQwVu3GMeidjdB40bb6+WB3LM6kvr2HQ1dclDz9v7f/PoJSY6AJPYzV+XbvSU9yUUL379wTgoR+28J0oiKEV9002kVnjKoJ0IQexuraegBITXBTQvctIWtzyQVgWmwU6+ak8vKJq+5ZrCt+hrXPqC6lq8ZBE3oYa2jvISE6grjoCKdDeUfFQZgx19o4OQTuXpLFpfoOzlS7aLGu3Ju0h67GRRN6GKtv6yY10UW9c7ASeu6qkDV3x+JMROC3x6+GrM1R5a6C5jJoc9lMVuV6mtDDWH17D6mJ0U6H8Y6WKmittHqoIZKRFEtxQQq/PeGmhO77/bXsosZIE3oYq2vrITXBRQndXz8PYUIHuKsoi9KqFi7Xu2RziezlIB6oOOB0JGqS0YQexurbut11Q7TiAHgiIWtpSJu9q8ha0fFlt/TSoxMgfZEmdDVmoyZ0EXlKRGpE5PgIr4uIPCEi50TkqIiErgCqxs0YQ4PbSi4VByFjMUTFhbTZ/BnxLMmd5r46euVBcMvoGzUpBNJD/xFw9w1evweY7/t6BPjOxMNSdmvp7KPPa9xzU9TrtRJYiMstfncXZXHwStPAZCvH5d70ziQrpQI0akI3xuwCGm5wyLuBnxjLW8B0EckOVoDKHnXt1qSiNLf00BsuQFdzSEe4DHb3Eqvs8opbyi7+90FvjKoxCEYNPRcoG/Rzue+564jIIyJSIiIltbUuWuUuDNX7JhXNcMtNUYduiPrNy0hibnoCL7ml7JKxGCJjNaGrMQlGQpdhnhu28GeM+Z4xptgYU5yenh6EptV41fum/bvmpmjFAYhKgPRCx0LYujSbty7UDyyJ4KiIKGu0S6UmdBW4YCT0ciB/0M95QGUQrqtsVN9u9dBdU3KpOGAlMI9zs1a3LcvGa1w0yShnlbW2TX+f05GoSSIYCf154GO+0S7rgGZjTFUQrqts5C+5pLih5NLfC1VHHauf+y3MtMouLx51ycc39ybo64TaU05HoiaJQIYtPgPsBRaKSLmIPCwij4rIo75DtgMXgHPA94FP2RatCpr69m6mx0cRFeGCqQjVJ6C/2/GELiJsW5bDvov11LS6YLTLwI1RHY+uAhPIKJcHjTHZxpgoY0yeMeaHxpgnjTFP+l43xphPG2PmGmOWGmNK7A9bTVR9W4/eEB3Gvb6yy8tuKLvMmAOxyZrQVcBc0D1TTqhr6ybNLTdEy0sgPtVaNtZhCzKTmJ+RyAtuKLuIQG6xJnQVME3oYcpVC3OVvQ15a6wE5gLblmXz9qUGatwwyShvNdSchO5WpyNRk4Am9DDlmmn/HQ1QfxbyVzsdyYB7l2VjDGw/5oJeev5qMF7tpauAaEIPQ339Xho7etwxBt2fqPLck9DnZSRRmJXkjrJLbrH1WL7f2TjUpKAJPQw1dvRiDO7ooZe9bS0Vm+OuNd3uXZZNyeVGyhs7nA0kbjqkLYQyTehqdJrQw1B9u4tmiZa/DZlFEJPodCTXePcKa/WK54+4YI5c/mqrh64rL6pRaEIPQ/5JRY730L39UH7AVeUWv/wZ8dxUkMKvDlU4v4F03hrobID6887GoVxPE3oYcs20/9rT0NNqJSwXes+KHM5Ut1Fa5fAIk3zf+1P+trNxKNfThB6GXLMwlz9B5bszoW9blkOkR/j14QpnA0lbCDHT9MaoGpUm9DBU39ZDhEdIjotyNpDy/RA3w5oR6UIzEqK5bUE6zx+pxOt1sOzi8VizaPXGqBqFJvQwVN/eTUp8NB6PwxN5yvZb9XOXTCgazrtX5lLV3MW+izfa4yUE8tdAzQmdYKRuSBN6GKpr63G+ft7ZCHWnXTWhaDh3LMokITrC+bJL3hrfBCNdH12NTBN6GKpv63Z+hMvAhCJ31s/94qIjuGtJFi8eq6Krt9+5QPJ8C5dpHV3dgCb0MNTQ7oJZomX7rQlFDi+ZG4j7V+XR2tXHy07uNxqXAmkLNKGrG9KEHobq21ywjsuVvZBRBDFJzsYRgHVzUslLieNnJeXOBpK/Fq68BV6vs3Eo19KEHma6evtp7e4j1cm10Pt6rCn/szY6F8MYeDzCB27K583zdc4uBVCwEbqaoLbUuRiUq2lCDzMN7f5Zog6WXKqOWFurzVzvXAxj9P6brKUAnjvg4M3RAt/7dXmPczEoV9OEHmYGpv072UO/4ktIBRuci2GM8lLi2Tg3jZ8dKHNuTPr0ApiWqwldjUgTepgZWJjLyR765T2QOh8SM5yLYRw+UJxHeWMnb12odyYAEesfwct7dKEuNSxN6GHG8R6612vdEJ1EvXO/u4qymBYbyX+XlDkXRMEGaLsKDReci0G5lib0MNPYYSX0FKcSes1J6GqelAk9NiqCP1iRw0vHr9Lc0etMEDN979uVvc60r1xNE3qYaeroJcIjTIuNdCaAy5Ovfj7Yg2tm0t3n5ecHHRrCmL7Q2lBb6+hqGJrQw0xTZw/JcVGIU+unXH4TkvNh+kxn2p+gopxkVs2cztNvXXZmnXQRa3TQ5TdD37ZyPU3oYaaxo5fpTq2yaIzVs5xEwxWH84frCrhQ186e8w7dHC3YAI2XoMUFuykpV9GEHmaaO3pJjncooTdcgPaaSVtu8du6NJuU+Ch+uveyMwH43z8tu6ghNKGHmabOHlLiHboh6i8TFEyOGaIjiY2K4IPF+bxaWs3V5q7QB5C5FKKTNKGr62hCDzNNTpZcLr0J8WmQNt+Z9oPow2tn4jWGZ96+EvrGIyJh5lqto6vraEIPM01OlVyMgQuvw+xbXb2hRaAKUhO4dX46z+6/Qm+/A4tlzboFak9BS1Xo21aupQk9jPT2e2nr7nOm5FJ7ypoQM3dT6Nu2ycfWF1Dd0s32Yw4k1bmbrccLr4e+beVamtDDSHOnNRlmuhM99POvWY9zpk5C37QwgznpCXz/jQuhH8KYucQqX13YGdp2latpQg8jTb5Zoo5sDn1+J6TOg+n5oW/bJh6P8Imb53C8ooW3LoR4z1GPB+bcbr2vuq6L8tGEHkaaOvw99BCXXPq6rRt4U6h37ve+VbmkJkTzgzccWFtl7iZrGGj1idC3rVwpoIQuIneLyGkROSciXxjm9dtFpFlEDvu+vhz8UNVE+RN6SqhLLmVvQ2/HO3XfKSQ2KoI/XFfAjlM1nKtpC23j/n8gteyifEZN6CISAfwbcA+wGHhQRBYPc+gbxpgVvq+vBDlOFQRN/hp6XIh76Bd2gkTArJtD226IfHR9AdGRHn64+2JoG07OtfYZPa8JXVkC6aGvAc4ZYy4YY3qAZ4F32xuWssNADT3UPfTzr0HeaoidFtp2QyQtMYb3r8rlFwfLqWvrDm3jczZZE4x6HZjgpFwnkISeCwxeALrc99xQ60XkiIi8JCJFw11IRB4RkRIRKamtrR1HuGoimjp68QgkxYRwpcWOBqg8PKWGKw7nT26ZQ2+/l++HupY+d5O1nV/ZvtC2q1wpkIQ+3CyQobfVDwIFxpjlwLeAXw13IWPM94wxxcaY4vT09DEFqiauqbOH6fHReDwhnNhz8feAmZI3RAebk57Ifctz+OneywP7tobErJvBE6l1dAUEltDLgcFjzfKAa5Z5M8a0GGPafN9vB6JEJC1oUaqgcGTa//mdEDMNcm8KbbsO+MzmeXT29od2xEtMEuSteWecvwprgST0/cB8EZktItHAA8Dzgw8QkSzxLbAtImt813VobVE1kpBP+/d64ewrMOc2a/2RKW5eRhLblmbz4z2XaAxlL33eFqg6ossAqNETujGmD3gMeBkoBf7bGHNCRB4VkUd9h90PHBeRI8ATwAPGkdX/1Y2EfKXFykPQWgWF94WuTYd9ZvN82nv6QzvipfBe6/H09tC1qVwpoHHoxpjtxpgFxpi5xph/9D33pDHmSd/33zbGFBljlhtj1hljdF1PFwp5yeXUb6zhigvuDF2bDluYlcTWpVn8aM+lgVFFtktfCDPmwqkXQtOeci2dKRpGQl5yOfWiddMuLiV0bbrAn26ZT3tPH//++vnQNCgChdvg4i7obApNm8qVNKGHCf9KiyGbVFR7BurOvFMOCCOFWdN438o8fvTmJcoaOkLU6L3g7YNzvwtNe8qVNKGHCf9KiykJIeqhn37ReizcGpr2XOZzdy1ABP7PK6dD02BeMSRkaNklzGlCDxP+dVxCttLiqRchewUk54WmPZfJTo7j4Ztn8+vDlRwtb7K/QU8ELLwHzr5qLYamwpIm9DDhv0EXkpUWW69C+X5YFH7llsE+eftcUhOi+aftpaFZL33RfdDTZtXSVVjShB4mBpbODUUP3T98Lgzr54MlxUbx2XfN560LDbx6str+BmffCtGJWnYJY5rQw4R/pcWQjEMvfQFmzIH0QvvbcrkH18xkYWYSf/ebk3T09NnbWGQMzL8DTm0Hb7+9bSlX0oQeJkK20mLrVWtdkcXvmRKbQU9UVISHf3zvEiqaOvnG787a32DRe61NL3RJ3bCkCT1MhGylxSPPgvHCio/Y284kUjxrBg+szueHuy9SWtVib2ML7rbG/R9+2t52lCtpQg8TIVlp0Rg4/F+Qvw7S5tnXziT0hXsKSY6L4ku/PIbXa+MN0sgYWPpBa5RRZ6N97ShX0oQeJkIy7b/iANSdhhUftredSWh6fDR/vXURh6408dO3Ltvb2IoPQ383HH/O3naU62hCDxMhmfZ/+GmIjLPquOo671uVy+0L0/mn7aWcq2m1r6Hs5ZC5BA5p2SXcaEIPE02dPfb20Hs74dhzsPgPpuxWcxMlInz1/mUkxETy2WcP09Pntashq5deeRBqSu1pQ7mSJvQw0dTRa++QxVMvQnez3gwdRUZSLP/yvqWcqGzha6+esa+hpR+0djLSm6NhRRN6mGi2u+Ry+GlIngmzbrGvjSnizqIsHlyTz3d3nWfPuTp7GklMh/l3wZH/B/299rShXEcTehjo7ffSaudKi9UnrC3QVn0MPPqRCsTfbFvM7LQEHnvmEOWNNq3IeNPHrTHpx35uz/WV6+jfvjDgX2lxul099N1ft6acr/mEPdefghJiIvn+x4rp7fPyyE8O2DOLdP4d1s3R3V+3tgNUU54m9DAwsI6LHQm94aI1PK74j8JuI4uJmpueyBMPrqT0agt/+fOjwV/ASwRu/p/WUFL/csZqStOEHgaaO21cafHNb1o339Y/Fvxrh4FNhRl8/q5CXjxaZc/SAIvfAymz4Y2vWRO/1JSmCT0MNLbbtNJi61XrZuiKj0BSVnCvHUYevW0O71+Vxzd3nA3+5tIRkbDxs9YQxguvB/faynU0oYcB21Za3Ptta9uzjX8a3OuGGRHhX9+/lHuWZPH3L5zk2bevBLeBFR+GxCzY/bXgXle5jib0MGDLSotNZbD/KSh6n7VUrpqQyAgP33xgJbctSOeLvzzGLw6WB/HiMbDhMWvji/OvBe+6ynU0oYeBE5UtpCZEMy02SCstGgPbPwcY2PLl4FxTER3p4ck/vIl1s1P58/8+wvd2nQ/ejdLVfwIz5sILf27N6lVTkib0Kc4Yw97z9ayfm4oEa33yk7+GM7+FTX8NKQXBuaYCIC46gh/98Wq2Lcvmn7af4u9+c5L+YKzOGBUL930DGi/C77868espV9KEPsVdrGvnaksXG+amBeeCnU3w0uetBaDWPhqca6prxERG8K0HVvLwzbP50Z5L/I+flgyUzSZk9q3WDew9T1iTwdSUowl9ittzvh6ADXNTg3PBHX8H7bVw3zetERTKFh6P8L/uXczf/UERvz9Ty9ZvvsH+Sw0Tv/Cd/wCxyfD8n+o2dVOQJvQpbu/5enKSYylIjZ/4xQ49DSVPwbpPQc7KiV9PjeqhDbP4xSc3EhXp4UPf3cvXXjlNV+8EEnH8DLj7X6CiBF7+ko5Nn2I0oU9hXq9h74V61s9Nm3j9/Mwr8PxnYM4m2PK3wQlQBWRpXjIvfOZm3r0ilydeO8edX9/FjtLqCVzwA7Du07DvSXjzG0GLUzlPE/oUdrq6lYb2HtZPtNxSXgI/ewiylsCHfgqRNi7Dq4aVFBvF1z+0gqc/sZboSA8P/7iEjz31NiXjKcOIWKWXJffD7x63tg1UU4Im9CnMXz+fUEI/8zI8fT8kZsJHfg4xSUGKTo3HxnlpbP/TW/jS1kKOVzRz/5N7+eB397KjtJq+/jEswOXxwHu+A3Nuh18/Bm89qeWXKUAT+hS293wds1LjyZ0eN/aTezrgxb+A//ogTMuFj/4SEjOCH6Qas+hID4/cOpfdf7WJL9+7mLKGDh7+cQnr/vk1/uGFkxwrbw5sI+rIaPjQf8K8d8Fv/wr+8/3Wcg5q0pJAJi6IyN3AN4EI4AfGmH8Z8rr4Xt8KdAAfN8YcvNE1i4uLTUlJyXjjVqPo6/ey8iuvcu/yHP75fUsDPm/t47/iHt7k8dTXoP6ctejWli9bsw0D8KHv7gXg//2P9QPPLX38ZQCOPX7XNd+PdK7f4GsE0t5wbY8W242eD9SNfqdQ6OnzsvN0Db84WM5rp2ro7TekJUZz87w0NsxLY3nedOamJxAZMUL/zRgo+SG8/DcQHQ/rP63r87iYiBwwxhQP99qo485EJAL4N+AOoBzYLyLPG2NODjrsHmC+72st8B3fo3LI8coWWrv7Rh+u2N0KNaeg5iSU7WMnPyeebvAsgo/+CuZuCkm8avyiIz3cVZTFXUVZNLb3sONUDbvP1rL7XB2/OlwJQEykh8KsJGalJTBzRjz5KfGkJ8WQmhjNjIRokpY8REL+RiJf/jzs+Aq89o+w4G6YcxtkLIKMImuETLAmpylbjNpDF5H1wOPGmLt8P38RwBjzz4OO+S7wujHmGd/Pp4HbjTFVI113vD30o68/x7RdOspiNF5j6Os3zElPJNIjVi/MeAFjbUnW0wbdbdDfPXBOG3G82LeWZ/s3cS6mEJCAe53+Xu6+i9ZNurWzZ1ByqYH4mEhau67fvCHJtwzBscfvuu5cv7WzZwDD95yHnpMU+047Q88bLrbBhj4faE/d3zP3tzv4d3IDr9dwoa6N4xUtHK9opvRqC5frO6hs6mSkikxMpIcFkVe5X3ay1ewincaB1/qIoJ14OiWOPonEi2DwoJX3sbs69wOs+8j48tiEeuhALlA26Odyru99D3dMLnBNQheRR4BHAGbOnBlA09eLTkimIX72uM4NN0kxUURmJr7zhHh8XxHWzc2YRIidDukLuefZOirIoMW3E30S2hOb7DweYV5GEvMyknjPytyB53v7vVQ1dVHX3k19Ww8N7d20dvXR3t1Pe08fPX0zOd1fzLHePyexp47MrvNkdl8kob+Z2P4OYr3tRJh+wOBBJyeNR2RSpi3XDaSH/gHgLmPMJ3w/fxRYY4z5zKBjXgT+2Riz2/fzDuDzxpgDI11Xa+juNNF6sNbQlbLXjXrogYxyKQfyB/2cB1SO4xillFI2CqSHHgmcAbYAFcB+4MPGmBODjtkGPIY1ymUt8IQxZs2Nrqs9dKWUGrsJ1dCNMX0i8hjwMtawxaeMMSdE5FHf608C27GS+TmsYYt/FKzglVJKBSag5fKMMduxkvbg554c9L0BPh3c0JRSSo2FzhRVSqkpQhO6UkpNEZrQlVJqitCErpRSU0RAi3PZ0rBILXB5nKenAXVBDCdY3BoXuDc2jWtsNK6xmYpxFRhj0od7wbGEPhEiUjLSOEwnuTUucG9sGtfYaFxjE25xaclFKaWmCE3oSik1RUzWhP49pwMYgVvjAvfGpnGNjcY1NmEV16SsoSullLreZO2hK6WUGkITulJKTRGuTegiMkNEXhWRs77HlBGOuyQix0TksIiUjPV8O+ISkXwR2SkipSJyQkQ+O+i1x0WkwhfvYRHZOsF47haR0yJyTkS+MMzrIiJP+F4/KiKrAj3X5rg+4ovnqIjsEZHlg14b9s80RHHdLiLNg/58vhzouTbH9ZeDYjouIv0iMsP3mp3v11MiUiMix0d43anP12hxOfX5Gi0uez9fxhhXfgFfBb7g+/4LwL+OcNwlIG2859sRF5ANrPJ9n4S1nvxi38+PA58LUiwRwHlgDhANHPG3M+iYrcBLgADrgH2BnmtzXBuAFN/39/jjutGfaYjiuh14YTzn2hnXkOPvA16z+/3yXftWYBVwfITXQ/75CjCukH++AozL1s+Xa3vowLuBH/u+/zHwnhCfP+7rGmOqjDEHfd+3AqVYe6wG2xrgnDHmgjGmB3jWF9/QeH9iLG8B00UkO8BzbYvLGLPHGOPfgfgtrF2u7DaR39nR92uIB4FngtT2DRljdgENNzjEic/XqHE59PkK5P0aSVDeLzcn9ExjTBVYCRLIGOE4A7wiIgfE2oR6rOfbFRcAIjILWAnsG/T0Y77/Cj41wVLQSJtzB3JMIOfaGddgD2P18vxG+jMNVVzrReSIiLwkIkVjPNfOuBCReOBu4LlBT9v1fgXCic/XWIXq8xUo2z5fAW1wYRcR+R2QNcxLfz2Gy2w0xlSKSAbwqoic8v0r6XRciEgi1l+8PzPGtPie/g7w91gfqr8H/i/wx+MNdZjnho5DHemYQM4dr4CvLSKbsP7C3Tzo6aD/mY4hroNYa2W0iXV/41fA/ADPtTMuv/uAN40xg3uBdr1fgXDi8xWwEH++AmHr58vRhG6MeddIr4lItYhkG2OqfP+FqxnhGpW+xxoR+SXWf112AQGdb1dcIhKFlcyfNsb8YtC1qwcd833ghUDjGsZENvCODuBcO+NCRJYBPwDuMcbU+5+/wZ+p7XEN+ocXY8x2Efl3EUkL5Fw74xrkAYaUW2x8vwLhxOcrIA58vkZl++fLjhsDwfgC/jfX3nz86jDHJABJg77fA9wd6Pk2xiXAT4BvDPNa9qDv/yfw7ARiiQQuALN550ZK0ZBjtnHtTau3Az3X5rhmYu1BuyHQP9MQxZXFOxPu1gBXfO+do++X77hkrPpsQijer0FtzGLkm3wh/3wFGFfIP18BxmXr5ytov0Swv4BUYAdw1vc4w/d8DrDd9/0c3y9+BDgB/PVo54corpux/rt0FDjs+9rqe+2nwDHfa88zKMGPM56tWKNozvt/f+BR4FHf9wL8m+/1Y0Dxjc4N4p/faHH9AGgc9P6UjPZnGqK4HvO1ewTrZtqGG50bqrh8P3+cIR2AELxfzwBVQC9WL/Jhl3y+RovLqc/XaHHZ+vnSqf9KKTVFuHmUi1JKqTHQhK6UUlOEJnSllJoiNKErpdQUoQldKaWmCE3oSik1RWhCV0qpKeL/A3hjvjSFhPTCAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dist_x = torch.distributions.Beta(2, 5)\n", "dist_y = torch.distributions.Beta(5, 5)\n", "x = dist_x.sample((15,))\n", "y = dist_y.sample((15,))\n", "z = torch.linspace(-0.5, 1.5, 100)\n", "pyplot.scatter(x, torch.zeros_like(x), marker='+')\n", "raw_density_x = dist_x.log_prob(z).exp()\n", "raw_density_y = dist_y.log_prob(z).exp()\n", "density_x = torch.where(raw_density_x.isnan(), torch.tensor(0.0), raw_density_x)\n", "density_y = torch.where(raw_density_y.isnan(), torch.tensor(0.0), raw_density_y)\n", "pyplot.plot(z, density_x)\n", "pyplot.plot(z, density_y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One technique to estimating (the density of) $P_X$ (with $P_Y$ being analogous) is via a kernel (Kernel Density Estimation). Given a kernel, for example the Guassian or Radial Basis Function (RBF - but note that sometimes RBF is considered a general class) kernel $K(x,y) = C(\\sigma, d) \\exp\\left(- \\frac{|x-y|^2}{2\\sigma^2}\\right)$, for some fixed standard deviation $\\sigma > 0$. We will not worry about the choice of $\\sigma$ for now.\n", "\n", "We can estimate $P_X$ as $\\hat P_X(x) = \\frac{1}{N_x} \\sum_{i=1}^N K(X_i, x)$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwFUlEQVR4nO3dd3yU153v8c9vRhWhBhIIkIQoondEN8axccM9xTYmceI4ccimOHnt7l3v3ddmsyW5d3M3u9nESby4xHbc0tyNE+zENgZMEaIjOkhICCGh3qWZc/94RiDQjDQjTdHM/N6vlzyjec4z80Mavj6cOc85YoxBKaVU+LOFugCllFL+oYGulFIRQgNdKaUihAa6UkpFCA10pZSKEDGheuGMjAyTl5cXqpdXSqmwtHv37mpjTKa7YyEL9Ly8PAoLC0P18kopFZZEpMTTMR1yUUqpCKGBrpRSEaLfQBeRHBH5QESKReSQiDzqps11IlIvIntdX98LTLlKKaU88WYMvQv4a2NMkYgkA7tF5D1jzOGr2n1sjLnd/yUqpZTyRr89dGNMhTGmyHW/ESgGxgW6MKWUUr7xaQxdRPKA+cAON4eXicg+EXlXRGZ6OP8RESkUkcKqqirfq1VKKeWR14EuIsOBPwDfMcY0XHW4CBhvjJkL/Ax43d1zGGM2GGMKjDEFmZlup1EqpZQaIK/moYtILFaYv2iMefXq4z0D3hizUUR+ISIZxphq/5WqwkZXO1QdhdozUFcKIybC1FtBJNSVKRXR+g10ERHgaaDYGPOfHtpkAZXGGCMii7F6/hf9WqkKDw3n4NnboebklY9PXQO3/RhSxoamLqWigDc99BXAF4ADIrLX9dj/BnIBjDFPAJ8Fvi4iXUArcL/RnTOiT9MFeO5OaKqEu5+AUdMhNQf2vQR/+QH8fIkV6nPuDXWlSkUkCVXuFhQUGL30P4I0X4TnbreGWT7/Bxi//MrjF0/CG9+AszvgwTdhwsqQlKlUuBOR3caYAnfH9EpRNXjGwCtrrdBe+3LvMAcYOQnW/Q5GTobffxkazwe/TqUinAa6Grzjm6ye95ofwcTrPLeLT4Z7n4eOJvj9w+DoClqJSkUDDXQ1OMbA5v9njZXPW9d/+1HT4fb/gpIt8MG/Bb4+paKIBroanNOboWwXrHgU7LHenTP3fpj/Bdj639b0RqWUX2igq8H5+D9g+GgroH2x+p8hNgn+/C+BqUupKKSBrgbu7E6rh778WxCb4Nu5SSOtXv2Rt63nUUoNmga6GrjN/wGJ6bDwoYGdv+yvIGkUvP99ayxeKTUoGuhqYJouWLNbFn0F4ocP7DnikmDV/4KSrXD8Pf/Wp1QU0kBXA3P0XcDAjLsG9zwLvwTpE6xeutPph8KUil4a6Gpgjr4LqbkwetbgnsceC9c9BhcOwakP/FObUlFKA135rqPZCl9/raA489PWTJkdTwz+uZSKYhroynenPoSuNpi2xj/PFxMHBV+2xuQvnuy/vVLKLQ105bsjGyE+Fcav8N9zLnwIbLGwc4P/nlOpKKOBrnzjdMCxP0L+jd5fGeqN5NEw69Ow50Vou3pDLKWUNzTQlW/KdkFLtf+GW3pa8jXoaIS9L/n/uZWKAhroyjdH3rGGRiav9v9zj1sI2Ytg5//oFEalBkADXfnm6LuQdw0kpAbm+Rc/AjWnrNUYlVI+0UBX3qsvh4vHIf+mwL3GtNshLhn2vhy411AqQmmgK++V7bJuc5YE7jXihsGse+DwG9DeFLjXUSoCaaAr75UXgj0esmYH9nXmPgCdzVD8ZmBfR6kIo4GuvFdWCGPmWBcCDUKXw8krO0v5wTuHOVhe37tB7lJrfRed7aKUT2JCXYAKE45OOLcXCga4VC5gjOHDo1X8cGMxxy80YRN48uPTzMlO5eFrJnDXvHFWQxGY9wB88AOoLYH08f75MygV4bSHrrxTeQi6WiG7YMBP8c9vHeahZ3fR6XDyxOcXUvSPN/L9O2bQ1ung0Vf28mpR2eXGc+6zbvf/ZpCFKxU9NNCVd7o/EB03sEDfeqKaZ7edYd2SXDZ9dxW3zMoibVgcX1oxgXe+vZKlE0fw2KsH2He2zjohfTzkrbSGXXTzC6W8ooGuvFO+29pdKC3X51NbOrp47NX9TMhI4h9vn0FczJVvu1i7jV+sW0jm8Hge+XUhFxrarAPzHoDa05f/Z6KU6pMGuvJO2S7rKs4BLJf7403HOFvTyv/59GwSYu1u24xIiuPJBwtoaO1i/Qu76XI4YdptYI+DQ68PsnilooMGuupfSw1cPAHZC30+dU9pLc9sPc26JbksnTiyz7Yzxqbww0/Poqi0jtf2lFtXo05eDYdf16UAlPKCBrrqX3mRdZu9yOdT/+2dYkYnJ/DYrdO8an/3vHHMzU7lJ+8fp63TATPvgYZyKNvp82srFW000FX/ynaB2GDsfJ9OK65oYHdJLV9ZOYHkBO+W2hUR/u6WaZTXtfLijlKYcot1MdOh1wZSuVJRRQNd9a+8EDKnQ3yyT6e9tKOUuBgbn1mQ7dN5yydnsDI/g59/cIJGEq211w+9rsMuSvVDA131zem0rhD1cf55c3sXr+0p57bZY0hP8v3K0r+9eSo1zR089fFpa9il6Tyc3e7z8ygVTfoNdBHJEZEPRKRYRA6JyKNu2oiI/FRETojIfhFZEJhyVdDVlUBbHYzz7Vf61r5zNLV38cAS36c5AszJTmPN7Cye+vgUNdnXQ0yCDrso1Q9veuhdwF8bY6YDS4FviMiMq9rcCuS7vh4BfunXKlXoVB21bkfN9Om0l3aWMmX0cArGpw/4pR+9YQrNHQ5e2VdjLdl7+A1rCzyllFv9BroxpsIYU+S63wgUA+OuanYX8LyxbAfSRGSM36tVwVd1xLrNnOL1KQfK6tlfVs8Di3ORAcxb7zY1K5kVk0fy609KcEy/G5oqoVSHXZTyxKcxdBHJA+YDO646NA442+P7MnqHPiLyiIgUikhhVVWVj6WqkKg6CsljfNqh6KWdJSTE2rjHxw9D3fnS8glU1LfxvmOuNdvlyNuDfk6lIpXXgS4iw4E/AN8xxly9Lbu7blivBTiMMRuMMQXGmILMzEzfKlWhUXUEMqd63bzT4eSd/RXcOmsMqYneTVXsy/XTRpEzIpGnd1TBpE9B8du6totSHngV6CISixXmLxpjXnXTpAzI6fF9NnBu8OWpkDLG6qFnendREEDhmVoa2rq4eeZov5RgtwlfXJbHzjM1lI++HupL4fwBvzy3UpHGm1kuAjwNFBtj/tNDszeBB12zXZYC9caYCj/WqUKhvszaOciHHvr7xZXExdhYme+/f4F9riCHxFg7T16YZl3gpMMuSrnlTQ99BfAF4HoR2ev6WiMi60VkvavNRuAUcAJ4EvirwJSrgqp7houXPXRjDO8drmTFpJEkxftv75TUxFg+s3AcLx1qoXPcEjjyjt+eW6lI0u/fOmPMFtyPkfdsY4Bv+KsoNURcmuHiXaAfv9BEaU0LX1s10e+lPLgsjxe2l7I7cTlLy34MNadhxAS/v45S4UyvFFWeVR2BpEwYNsKr5u8drgRg9XT/jJ/3NGV0Mgty03i8wjX8o8MuSvWiga488/ED0feLK5mbncrolISAlHP/oly2VA+nZcR0HXZRyg0NdOXepRku3n0geqGxjb1n6wLSO+9225wxJMXZ2RKz1LrAqOlCwF5LqXCkga7cazwP7fVe99A/OHIBY2D1jMAFelJ8DHfMHcsvz08DDBz7Y8BeS6lwpIGu3Lv0gah3PfT3DlcyLi2RaVm+LbHrq/sW5bCnI5vmhCw4qoGuVE8a6Mo9H6YsdnQ52XriIjdMHzWotVu8MS8njSmjk/mQAjj1AXS2BvT1lAonGujKvaojkJhuzXLpx76yOlo7HayYnBHwskSE+xbl8puGmdDZAqc3B/w1lQoXGujKve4ZLl70uLeduIgILJ3Q9ybQ/nLP/HHslpm024bB0XeD8ppKhQMNdNWbMVBV7PX4+baT1cwam0rqsMEvxuWNEUlxrJyWzVYzG3PsT7pYl1IuGuiqt+ZqaK31avy8tcPBntI6lk8KTu+82+cKstnYMQ9pPAcV+4L62koNVRroqreLJ6zbkZP7bbq7pJYOh5NlQQ70VVMy2ZuwBCei0xeVctFAV73VlVq36Xn9Nt12spoYm7Aoz7vlAfwlxm7j+oUz2OPMp6tYrxpVCjTQlTt1JdZtak7f7YBtJy8yLyfNr6sreuuzC7N537GAmMr90KDL7yulga56qy2B4VkQ2/eaLA1tnewvC/74ebcpo5M5m3ktAEYvMlJKA125UVcC6eP7bbbzVA1OA8smBX7+uSdLlqzgrDOTxgM67KKUBrrqrbYE0voP9G0nLxIfY2N+blrga/Lgzrnj+Ij5JJRt0atGVdTTQFdXcnRCQ5lXPfRtJ6spyEsnIdYehMLcSx0WS2PuDcQ52+g48VHI6lBqKNBAV1eqLwPj7LeHXtfSwZHzjUG7OrQvc1feQYuJp2znG6EuRamQ0kBXV+qe4ZKW22ez3SW1ACyaENzpiu4szR9LkX0Ow0v/rFeNqqimga6udGkOet899F1naom1C3Oz0wJfUz9sNqFt4mpGOSqpOLE31OUoFTIa6OpKtSUgdkjJ7rPZ7pIaZo1LJTEudOPnPc289rMAnNjyhxBXolToaKCrK9WVQOo4sHu+UKit08G+s/VBvzq0L2NyJ3MmdhLDS/+Mw6nDLio6aaCrK3kxZfFgeT0dDicF49ODVJR3OiauZrbzCDsOnwx1KUqFhAa6upIXFxXtOmN9ILpwiAX6+KX3ECNOjmx5PdSlKBUSGujqss5WaKqEtLw+mxWeqWFSZhIjh8cHpy4vxY9fTHNMGiPKP6CqsT3U5SgVdBro6rLuGS59TFl0Og2FJbVDavz8Epsdx8TrWWnbx+8Lz4S6GqWCTgNdXVbrmoPex5DLiaom6ls7KRiKgQ6kzL6NkdLIvu0f4NQPR1WU0UBXl126qMhzoO86UwPAoryhNX5+yeQbMNiY0fwJW09Wh7oapYJKA11dVlcC9ngYPtpjk8IztWQmx5M7YlgQC/NBYjomZzE3xuzjpR2loa5GqaDSQFeX1ZZY4+c2z2+LXWdqWJSXjogEsTDf2KbczHROs+9wMRca20JdjlJBo4GuLutnyuL5+jbKaltZOH5ojp9fMuVmAK6Rvfx219kQF6NU8PQb6CLyjIhcEJGDHo5fJyL1IrLX9fU9/5epgqK7h+5BUenQnH/ey6gZkJLN55IP8+KOUroczlBXpFRQeNNDfxa4pZ82Hxtj5rm+/mXwZamga6uHtro+PxAtKqklPsbGjDEpwatrIEQg/0bmd+3hYn0jmw5XhroipYKi30A3xmwGaoJQiwolL6Ys7i6tZU52KnExYTBSN+VmYrpauDXlDM9uOxPqapQKCn/9zVwmIvtE5F0RmempkYg8IiKFIlJYVVXlp5dWftFQbt2muh9yae9ycKi8gQW5Q3y4pduEa8Eez8Ojj7HzdA3FFQ2hrkipgPNHoBcB440xc4GfAa97amiM2WCMKTDGFGRmZvrhpZXfNJyzblPGuD18sLyBDoeT+eES6HFJMGElM5u2kxBr4zntpasoMOhAN8Y0GGOaXPc3ArEiErpt4NXANJ4HsUHSKLeHi1w7FC0YnxbEogYp/2bstSf5ygwnr+8tp66lI9QVKRVQgw50EckS16RkEVnses6Lg31eFWSN56ww97AOelFpLTkjEhmVnBDkwgZhyk0AfGHkUdo6nbyiUxhVhPNm2uLLwCfAVBEpE5GHRWS9iKx3NfkscFBE9gE/Be43Rjd2DDuN5yE5y+0hYwy7S2rDZ/y8W3oeZExl9PmPWDZxJM9uPUNHl05hVJHL87Y0LsaYtf0cfxx43G8VqdBoPO9xDnp5XSsXGtuH/vxzd6bcBNuf4Ov3jOLBF4t558A57pnf9/Z6SoWrMJh/poKi4ZzHHnpRaR1A+PXQAfJvBmcn19gPkj9qOBs2n0b/AakilQa6gq52aK2B5LFuDxeV1JIYa2daVnKQC/OD3KUQn4rt+Ca+unIixRUNbD2hH/GoyKSBrqCxwrr12EOvZW5OKjH2MHy72GNh8vVwfBN3zcsiY3g8Gz4+FeqqlAqIMPwbqvyu8bx162YOemuHg8PnwuiCInfyb4amSuKrDvLQijw2H6vSC41URNJAV5cvKkruHej7y+rocpowD/QbAYFjm1i3JJfEWDtPai9dRSANdHW5h+4m0Ls/EJ2fmxa8evwtKQOyC+DYH0kbFsfaxbm8sfccZ2taQl2ZUn6lga6si4rs8ZDYuxdeVFpL3shhjBweH4LC/GjKzXCuCBrP89VrJ2AT2LBZe+kqsmigK6uHnjLGWna2B2MMe0rD8IIid6bcat0e38SY1EQ+syCb3xSe1R2NVETRQFfQUOF2uOVsTSvVTR3MD8cLiq42eiakZMPRPwKwftUkuhxOnv74dIgLU8p/NNCVNW3R7fi5a0GucB4/7yYCU2+BUx9AZxt5GUncPmcsL2wv0UW7VMTQQI92xvQZ6MPi7EwdHYYXFLkz5RbobIEzHwPw9esm0dzh0A0wVMTQQI927Q1WyLmZg15UWsu8nLTwvKDInbyVEDsMjlnDLtPHpLB6+mh+tfUMjW2dIS5OqcGLkL+pasAauq8SvTLQWzq6KK5ojIwPRLvFJsDET1nj6K71XB69IZ/61k7dAENFBA30aNfoPtD3l9XjcJrw2tDCG1NvgYYyqDwEwOzsVK6fNoqntpymqb0rxMUpNTga6NHOwzou3R+Izs+JoB46QL616QXH3r300KM35FPXor10Ff400KOdhx56UUkdEzOSSE+KC0FRAZScBWMXwNHLgT43J43rpmby1MenaNZeugpjGujRrqECElIhbtilh4wxFJXWhs+G0L6atgbKd19e8gCrl17b0snzn5SEsDClBkcDPdo1VvRaB73kYgs1zR2RN37ebeoa67ZHL31+bjrXTsnkSe2lqzCmgR7tGit6jZ/vLrHGz8NyyzlvjJoBaePh6MYrHn70hnxqmjv49XbtpavwpIEe7RrPQ8qVPfTdpbUkx8cwZVSEXFB0NRGYdhuc+gjamy49vHC81UvfsFl76So8aaBHM6fDCvSre+hnapk/Ph2bTTycGAGmrgFHO5z8yxUPf2e11UvXsXQVjjTQo1lzNRjHFTNc6ls7OXahkYJIHW7plrvMWi74qmGXBbnprJqSyYbNJ7WXrsKOBno0czNlcU9pLcZE8Ph5N3uMtTXdsT+C48rg/s5qnfGiwpMGejRzE+hFJbXYBOblpIWmpmCaeiu01sLZ7Vc8PD83neumai9dhR8N9GjWdMG6HT7q0kOFJbVMH5NCUnxMiIoKosk3gD0OjmzsdejbrnnpL+iMFxVGNNCjWUu1dZuUAUCXw8nes3WRP9zSLT4ZJl4HR9+5tFhXtwW56azMz2DD5lO0djhCU59SPtJAj2bN1RCXDLGJABw530hLhyN6Ah1g+h1QewYqD/Y69OgN+Vxs7uDFHdpLV+FBAz2aNVdd6p1DFFxQ5M7UNSA2KH6r16GCvBEsnzSS/9l8irZO7aWroU8DPZo1V0FS5qVvC0tqyUpJYFxaYgiLCrKkDBi/wm2ggzWWXtXYzis7S4NcmFK+00CPZs3VV/TQi0pqWTg+HZEIvqDInWm3w4XDUH2i16GlE0eyeMIIfvnRSdq7tJeuhjYN9GjWI9Ar6lspr2uNruGWbtNvt26PuO+lf+v6yVQ2tPNaUXkQi1LKd/0Guog8IyIXRKT3p0bWcRGRn4rICRHZLyIL/F+m8jun05rl4hpyKTwThePn3VKzrTXSPQy7XDM5g9njUvmfzadwOI3bNkoNBd700J8Fbunj+K1AvuvrEeCXgy9LBVxbHTi7LgX6ztM1JMXZmTk2JbR1hcr0O6w10uvLeh0SEb5+3SROVzfz7sGKEBSnlHf6DXRjzGagpo8mdwHPG8t2IE1Eem8hr4aWZtcc9GHWkMuuMzUsGJ9OjD1KR+Gm32ndHnnH7eGbZ2YxMSOJX3xwEmO0l66GJn/87R0HnO3xfZnrsV5E5BERKRSRwqqqKj+8tBqwZtfPPymD+pZOjlY2sihvRGhrCqWMyZA5HQ6/6faw3SasXzWJwxUNfHRM37tqaPJHoLubEuG2C2OM2WCMKTDGFGRmZrprooLlUqBnUlhSgzGweEIUBzrAjLugZCs0Vro9fPf8cYxJTeAXH54McmFKeccfgV4G5PT4Phs454fnVYF06bL/THaeriHWLtGxIFdfZt4NGCh230uPi7HxlZUT2Xm6hr1n64JZmVJe8Uegvwk86JrtshSoN8boJ0dD3aUx9JHsPFPDnOw0EmLtoa0p1EZNt4ZdDr3mscl9i3JIjo/h6S2ng1iYUt7xZtriy8AnwFQRKRORh0VkvYisdzXZCJwCTgBPAn8VsGqV/zRXQWI6rQ7hQFl9dI+f9zTzHijZBg3u+yTD42NYuySXjQcqKK9rDXJxSvXNm1kua40xY4wxscaYbGPM08aYJ4wxT7iOG2PMN4wxk4wxs40xhYEvWw2a67L/PWdr6XIaFk+Iwvnn7vQz7ALwxeV5ADy37UwwKlLKa1E6R01ZV4lmsut0LSKwcLz20AHInAqjZsLBVz02GZeWyK2zsnh5RylNugGGGkI00KOV67L/nWcuMi0rhdTE2FBXNHTMvMfaxaje86X+X1k5kcb2Ln6766zHNkoFmwZ6tGquwjEsg6KSOhbn6XDLFWbeY90efsNjk3k5aRSMT+eZrad1OQA1ZGigRyNHF7TWUOVIprXTwaJon39+tYzJkDUbDnkedgF4+JoJlNW28udi9/PWlQo2DfRo1HIRgOPNCYBeUOTWzE9D2S5rNyMPbpwxmqyUBH6t+46qIUIDPRq5rhLdWxPL5FHDGZWcEOKChqDZn7VuD/zOY5MYu40HluTy8fFqTlU1BakwpTzTQI9GrqtEd16wsXzSyBAXM0Sl5ULuctj/214bSPd0/+IcYu3CC9t1RyMVehro0ch1lei5zuEa6H2Z8zmoPgYV+zw2GZWcwC2zxvC73Wdp6dApjCq0NNCjkWvIpYYUlkzQQPdoxt1gi+1z2AXgwWXjaWzr4o29uoSRCi0N9GjUXIUDG2NHjyE9KS7U1Qxdw0ZA/k1w4Pfg9LyfaMH4dKZlJfP8JyW6VroKKQ30KNTVeIGLJoWlk3UJ437N+Rw0nYczH3tsIiI8uCyP4ooGikprg1icUlfSQI9C9dUVXDQpOn7ujSm3QFwy7O972OWueWNJirPz8k69clSFjgZ6FGqvr+QiKXpBkTdiE62NLw6/AR0tHpslxcdw1/xxvL3/HPWtnUEsUKnLNNCjkLRU40gcSUqCrt/ilXlroaMRit/qs9naRbm0dTp5Y6/nNWCUCiQN9CjT0tFFclcdSelZoS4lfIxfAekTYM+v+2w2OzuVWeNSeGlHqX44qkJCAz3KFJ48z3BpZeRot/t4K3dEYN4664PRmr53Klq7OJcj5xvZV1YfpOKUukwDPcrsLT4OwNhxuSGuJMzMWwsI7Hu5z2Z3zh3LsDg7L+/QK0dV8GmgR5mjp04BEJ8yOsSVhJnUbJh0Pex5sc856ckJsdw5dyxv7jtHY5t+OKqCSwM9ilTUt9Jcc976JikjtMWEo/nroKEMTn/UZ7P7F+fS2ungrX26V7oKLg30KLL5WBUjabC+GaZz0H029TZISIM9L/TZbG52KlNHJ/ObQp2TroJLAz2KfHSsitzENusb7aH7LjYB5twLxW9DS43HZiLCvYty2He2jqPnG4NYoIp2GuhRosvhZMvxamaNcIDYIT4l1CWFpwVfBEc77H2xz2b3zB9HrF34je45qoJIAz1K7Curo6Gti8nDOyExzZqKp3yXNQtylkLhM+B0emw2IimOG2eM5rU9ZXR0eW6nlD9poEeJj45VYxMYE9cKiXrJ/6As+grUnIJTH/TZ7N6CHGpbOnlf9xxVQaKBHiU+OlbFvJw04jrrITE91OWEtxl3wrAMq5feh5X5mYxNTdBhFxU0GuhRoKa5g/1ldayaMsr6ME8DfXBi4mHBF+DoRqgv89jMbhM+uzCbzcerOFfXGsQCVbTSQI8CW05UYwxcOyUDWus00P1h4UPWXqO7n+uz2ecKcjAGfr/bc/Ar5S8a6FHgz8WVjEiKY052GrTWWjvxqMFJH2/tZlT0HHR1eGyWM2IYyyeN5LeFZ3E6dcEuFVga6BGuo8vJX45cYPX0UdhNl7UMrPbQ/WPxV6GpEg691mez+xblUFbbyienLgapMBWtNNAj3I7TF2ls6+KmGVlW7xw00P1l0g2QMRU++Zk1/OLBzTOzSE2M1Q9HVcBpoEe4TYcqSYy1c01+hga6v9lssOwbcP4AnN7ssVlCrJ27543lj4fOU9+iC3apwNFAj2DGGN47XMm1UzJIiLVroAfCnPsgKRM+ebzPZvcuyqGjy8nrupuRCiCvAl1EbhGRoyJyQkQec3P8OhGpF5G9rq/v+b9U5asD5fWcb2izhltAAz0QYhNg0Vfh+CaoOuqx2cyx1m5Gr+w6q7sZqYDpN9BFxA78HLgVmAGsFZEZbpp+bIyZ5/r6Fz/XqQZg06FK7Dbh+mmjrAe6F5TSQPevRQ9DTEK/vfT7CnIormjgYHlDkApT0cabHvpi4IQx5pQxpgN4BbgrsGUpf9h0+DyL80aQnhRnPdDdQ9dpi/6VlAFz18K+30Cj58v875w3joRYGy/v0t2MVGB4E+jjgJ4fz5e5HrvaMhHZJyLvishMd08kIo+ISKGIFFZVVQ2gXOWt09XNHKts4qaZPXYmaq3VlRYDZfm3wNlpzXjxIDUxlttmj+WNPeU0tXcFsTgVLbwJdHfL8l09CFgEjDfGzAV+Brzu7omMMRuMMQXGmILMzEyfClW+2XTI2pnoxhlXBbqutBgYIyfB7M/BrqehyXNn5YEluTR3OHhDPxxVAeBNoJcBOT2+zwbO9WxgjGkwxjS57m8EYkVEd1AIobf2n2NOdirZ6cMuP9iq67gE1LV/C52tffbSF+SmMS0rmZd2lOqHo8rvvAn0XUC+iEwQkTjgfuDNng1EJEvE6vaJyGLX8+plcSFy4kIjB8sbuHveVSNjrbW6dG4gZeTDrM/Azqeg2f3bX0RYtySXQ+ca2F9WH+QCVaTrN9CNMV3AN4E/AcXAb40xh0RkvYisdzX7LHBQRPYBPwXuN9r9CJnX95zDJnD73DFXHmit1R56oF37t9DZ0ueMl7vmjyMx1s5LO/TDUeVfXs1DN8ZsNMZMMcZMMsb8wPXYE8aYJ1z3HzfGzDTGzDXGLDXGbAtk0cozYwyv7y3nmvxMRiUnXHmwRQM94EZNg5l3w84NHvcdTUmI5c65Y3lz3zka2vTKUeU/eqVohNldUktZbSt3zxvb+6D20INj1d9BRzN8/GOPTR5Ykktrp4PXivTDUeU/GugR5rU95STE2rhpZtaVBxyd1kqLOgc98EZNh3nrrF567Rm3TeZkpzI3O5Xntp3RZXWV32igR5COLifvHKjgphlZDI+PufJga511qz304Lj+H6w5/392f9G0iPDQigmcqm7mo2N6TYbyDw30CLL5WBV1LZ3cPd/dcIte9h9UKWNh+Tfh4B+gfLfbJmtmj2FUcjzPbD0d5OJUpNJAjyC/232WEUlxrMx3c9HWpYW50oJaU1Rb8ai1EuOm77ldLz0uxsaDy8bz8fFqjlU2hqBAFWk00CPEubpW3jtcyb0FOcTa3fxaLwW6jqEHTXwyXPcYlGyBI2+7bbJ2cS7xMTZ+tfVMcGtTEUkDPUJ0z2letyTXfQNdOjc0FnwJRs+Cdx+D9qZeh0cOj+fueeN4taiM2mbPe5Mq5Q0N9AjQ3uXg5Z2l3DB9NDkjhrlvpEvnhoY9Bm7/L2gohw//j9smD12TR3uXk5d26oVGanA00CPAxgMVXGzu4MFl4z03aq0FselKi6GQsxgWfgm2/xIq9vc6PC0rhVVTMnl6y2madRVGNQga6BHguW0lTMxMYsWkPtZD676oyKa/8pBY/U/Wz//t74LT0evwt2/Ip6a5gxe2l4SgOBUp9G93mNtfVsfes3U8uHQ8Nlsfy+LqVaKhlZgON/8Qygth55O9Di8cn87K/Aw2bD5Fa0fvwFfKGxroYe7ZrWdIirPzmYXZfTfUpXNDb869kH8TvP9PcKG41+Fv35DPxeYOXtyhvXQ1MBroYexkVROv7y1n7eJckhNi+26sS+eGngjc9XOIGw6vfhW62q84vChvBMsnjeSJj7SXrgZGAz2M/ff7x0mItbP+ukn9N9Yhl6Fh+Ci463E4fwD+8m+9Dn/7hnyqm9q1l64GRAM9TB0938hb+8/xxeV5ZAyP7/8EXTp36Jh6Kyx8CLb9DE59eMWhpRNHsnzSSH72lxM6L135TAM9TP3k/WMMj4vha9dO7L9x90qLGuhDx80/sHY4+v2XofbK3vj37phBY1sn//nesRAVp8KVBnoYOlhez7sHz/PlayaQNiyu/xO6V1rUpXOHjrgkuP9lcHTBK+us9dNdpmWl8IWl43lxRwmHzzWEsEgVbjTQw4wxhh/96SipibE8vHKCdyfpZf9DU8Zk+OwzcOEQvP5XVyzg9d0bp5CaGMv33zqkm0krr2mgh5k39p5j87EqvrM6n5T+ZrZ0u7R0blrA6lIDlL8aVn8fDr8OH/zw0sNpw+L4m5unsvN0DW/vrwhZeSq8aKCHkYtN7fzzW4dYkJvGg8vyvD9Re+hD2/Jvw7zPw+YfwbbLm0vfvyiXmWNT+P6bh7jQ2BbCAlW40EAPI//81mGa2x38+2fmYO/rqtCr6dK5Q5sI3PHfMOMu2PQPUPgrAOw24b/um0dTexd//dt9ulWd6pcGeph4/3Alb+47xzevn0z+6GTfTtaVFoc+ewx8+inrStK3vwt7XwZgyuhk/vH2GXx8vJqnt+jORqpvGuhhoKy2hcdePcC0rGTWr/LiIqKrnd8PwzIgIdX/xSn/iYmDe5+HCdfC6+thy0/AGNYtyeXmmaP50Z+OcKCsPtRVqiFMA32Ia2zr5OFnC2nvcvD4A/OJi/HxV2YMnN5shYT4MEyjQiM2Edb9DmZ+2lrz5d2/Q4yTf//MHDKGx7P+hd1U1LeGuko1RGmgD2FdDifffGkPJ6uaeOLzC5k8ysehFoCLJ6Cxwgp0FR5i4uEzT8Oyb8LO/4GX15JGE08+WEBDayeff2oHF5va+38eFXU00Icop9PwT28e4qNjVfzr3bNYMbmPtc77cnqzdauBHl5sNutq0jX/ASf/Ak9cw6yuQzz9pUWU1bbyxV/tpKGtM9RVqiFGA30Iaut08K1X9vDijlLWr5rE2sUe9gn1xunNkJINI7xYIkANPYu/Cl95z+q1P3sbi0//kg1rZ3KkopEHn95JZYNOZ1SXaaAPMTXNHax7agfv7K/g72+dxt/dMnXgT+Z0wpmPdfw83I2dD1/bDHPug80/YtX7d/DbT9Vz9Hwjt/9sC7vO1IS6QjVEaKAPIR8cvcCdj2/hYHk9v1i3gK+tmoQMJogvHIaWizrcEgnik+GeJ+DBN8Eez4Kt6ynM+wVL7MdYu2E7T24+RZfDGeoqVYhpoA8B5XWtfO3XhTz0q13Exdh45ZGlrJk9ZvBPfGn8fOXgn0sNDRNXwde3wk3/RtLFgzze9vdsTP4h2//4Amt+8iF/Lq7UtV+imITql19QUGAKCwtD8tpDgTGGotJaXtheyjsHKrAJfOv6fL66cqLvUxM9eel+qD4K397jn+dTQ0tHMxQ9j9n2M6ShnBpJ49XOZRwfvYZlKz7FzbPGkhhnD3WVys9EZLcxpsDtMW8CXURuAf4bsANPGWP+71XHxXV8DdACfMkYU9TXc0ZjoLd1Oth1poYtJ6r58EgVRysbGR4fwz3zx/G1VRPJTh/mvxdzdMGPJsCsT1uXlavI1dUBxzfh3Psy5tgfsZsuqkwKO5hD87hrGDl1KdNnFzBuxACmvaohp69Aj/HiZDvwc+BGoAzYJSJvGmMO92h2K5Dv+loC/NJ1GxWMMbR3OWlq76KprYv61k6qGtu50NjO+fpWTlY1c/xCI6erm+l0GGLtwvycdH54z2zumjeWpPh+fw2+q9gH7Q06fh4NYuJg+u3Ypt8OLTU4j74L+//EqtLNJJ/bAueg7S+xHLaNpy4xh87kXGwjJhCfnkViWhbDR2QxLHkEickpJCUk+LZOkBpSvEmSxcAJY8wpABF5BbgL6BnodwHPG6u7v11E0kRkjDHG7+t+7v/wD6Rs/id/P61HV/8DxrgeNK77xtBrzDIRyHV9AcTabcTF2IhPsZEYaycxzo7NIVCI9RUIba5LxPN0/DyqDBuBbf46MuevA6cTR9UxKo7soO7kLmKrDzGp7TAZzR9ir3T/L/M2E0sTcXQSQ6fE4MSOExtGBHPpIzfh0tly6T/KB+cnfY6l6/yfY94E+jjgbI/vy+jd+3bXZhxwRaCLyCPAIwC5uQObWx2XlErNMC83driC72866XFHAHHdsYl1XwRsNsEugk0gxm4jxibE2m3Ex9iIj7URH2MnZB2eUTOtTYlVdLLZsI+eRvboaWSv+uKlh01XB7XnT9NwsYKWmvO0N1zA2daAaWvEdDQhjg5wdFi3TgdinGAcrpOtroy5dH8g9EPbmOTRgXleL9q4i6OrfyPetMEYswHYANYYuhev3cu0Rath0eqBnKqUAiQmjvTsqaRnD+IaBzUkeTOdogzI6fF9NnBuAG2UUkoFkDeBvgvIF5EJIhIH3A+8eVWbN4EHxbIUqA/E+LlSSinP+h1yMcZ0icg3gT9hTVt8xhhzSETWu44/AWzEmrJ4Amva4kOBK1kppZQ7Xs2XM8ZsxArtno890eO+Ab7h39KUUkr5Qi/9V0qpCKGBrpRSEUIDXSmlIoQGulJKRYiQrbYoIlVAyQBPzwCq/ViOvwzVumDo1qZ1+Ubr8k0k1jXeGJPp7kDIAn0wRKTQ02pjoTRU64KhW5vW5RutyzfRVpcOuSilVITQQFdKqQgRroG+IdQFeDBU64KhW5vW5RutyzdRVVdYjqErpZTqLVx76Eoppa6iga6UUhFiyAa6iIwQkfdE5LjrNt1DuzMickBE9opIoa/nB6IuEckRkQ9EpFhEDonIoz2OfV9Eyl317hWRNYOs5xYROSoiJ0TkMTfHRUR+6jq+X0QWeHtugOta56pnv4hsE5G5PY65/Z0Gqa7rRKS+x+/ne96eG+C6/rZHTQdFxCEiI1zHAvnzekZELojIQQ/HQ/X+6q+uUL2/+qsrsO8vY8yQ/AJ+BDzmuv8Y8O8e2p0BMgZ6fiDqAsYAC1z3k4FjwAzX998H/sZPtdiBk8BEIA7Y1/06PdqsAd7F2lVqKbDD23MDXNdyIN11/9buuvr6nQapruuAtwdybiDruqr9HcBfAv3zcj33tcAC4KCH40F/f3lZV9DfX17WFdD315DtoWNtPP2c6/5zwN1BPn/Az2uMqTDGFLnuNwLFWHus+tulDbyNMR1A9wbeV9f7vLFsB9JEZIyX5wasLmPMNmNMrevb7Vi7XAXaYP7MIf15XWUt8LKfXrtPxpjNQE0fTULx/uq3rhC9v7z5eXnil5/XUA700ca165Hr1tNuxwbYJCK7xdqE2tfzA1UXACKSB8wHdvR4+Juufwo+M8ihIE+bc3vTxptzA1lXTw9j9fK6efqdBquuZSKyT0TeFZGZPp4byLoQkWHALcAfejwcqJ+XN0Lx/vJVsN5f3grY+8urDS4CRUTeB7LcHPoHH55mhTHmnIiMAt4TkSOu/0uGui5EZDjWX7zvGGMaXA//EvhXrDfVvwI/Br480FLdPObtBt5ebew9QF4/t4h8Cusv3DU9Hvb779SHuoqw1spoEuvzjdeBfC/PDWRd3e4AthpjevYCA/Xz8kYo3l9eC/L7yxsBfX+FNNCNMas9HRORShEZY4ypcP0T7oKH5zjnur0gIq9h/dNlM+DV+YGqS0RiscL8RWPMqz2eu7JHmyeBt72ty43BbOAd58W5gawLEZkDPAXcaoy52P14H7/TgNfV43+8GGM2isgvRCTDm3MDWVcP93PVcEsAf17eCMX7yysheH/1K+Dvr0B8MOCPL+D/ceWHjz9y0yYJSO5xfxtwi7fnB7AuAZ4HfuLm2Jge978LvDKIWmKAU8AELn+QMvOqNrdx5YdWO709N8B15WLtQbvc299pkOrK4vIFd4uBUtfPLqQ/L1e7VKzx2aRg/Lx6vEYenj/kC/r7y8u6gv7+8rKugL6//PaH8PcXMBL4M3DcdTvC9fhYYKPr/kTXH3wfcAj4h/7OD1Jd12D9c2k/sNf1tcZ17NfAAdexN+kR8AOsZw3WLJqT3X9+YD2w3nVfgJ+7jh8ACvo614+/v/7qegqo7fHzKezvdxqkur7pet19WB+mLe/r3GDV5fr+S1zVAQjCz+tloALoxOpFPjxE3l/91RWq91d/dQX0/aWX/iulVIQYyrNclFJK+UADXSmlIoQGulJKRQgNdKWUihAa6EopFSE00JVSKkJooCulVIT4/x2qpE7UA4f4AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sigma = 0.1\n", "raw_hat_p_x = torch.exp(-((x[None] - z[:, None])**2)/(2*sigma**2)).sum(1)\n", "hat_p_x = (raw_hat_p_x / raw_hat_p_x.sum() / (z[1]-z[0]))\n", "pyplot.plot(z, hat_p_x)\n", "pyplot.plot(z, density_x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given our two samples, $X_i$ and $Y_i$, we may do this for both, so we have two probability densities $\\hat P_X$ and $\\hat P_Y$. Now, even if $P_X = P_Y$, we would not expect exact equality $\\hat P_X = \\hat P_Y$ due to the sampling and estimation error. But we would expect that they are somehow _close_.\n", "\n", "To measure this closeness, we can look at our kernel again. It turns out that this is a so-call [reproducing kernel](https://en.wikipedia.org/wiki/Reproducing_kernel_Hilbert_space), that is if we define the two functions $K_x, K_y : \\mathbb{R}^d \\rightarrow \\mathbb{R}$, $K_x(y) := K(x,y) =: K_y(x)$ we can define a scalar product on the Hilbert space spanned by such functions by setting $\\langle K_x, K_y \\rangle_{\\mathcal{H}} := K(x,y)$, this is called a reproducing kernel Hilbert space in the literature. It turns out that for any function $f \\in \\mathcal{H}$, we have $f(x) = \\langle f, K_x \\rangle_\\mathcal{H}$.\n", "\n", "Clearly, our estimated distributions $\\hat P_X$ and $\\hat P_Y$ are in $\\mathcal{H}$, so we can measure closeness in the norm of $\\mathcal{H}$.\n", "\n", "This leads to the maximum mean discrepancy test statistic $\\widehat{MMD} = \\| \\hat P_X - \\hat P_Y \\|_{\\mathcal{H}}^2$. By the the above representation, we can write\n", "\n", "$$\n", "\\begin{aligned}\n", "\\widehat{MMD} &= \\| \\hat P_X - \\hat P_Y \\|_{\\mathcal{H}} \\\\\n", " &= \\langle \\hat P_X - \\hat P_Y, \\hat P_X - \\hat P_Y \\rangle_{\\mathcal{H}} \\\\\n", " &= \\langle \\hat P_X , \\hat P_X \\rangle_{\\mathcal{H}} + \\langle \\hat P_Y, \\hat P_Y \\rangle_{\\mathcal{H}}\n", " - 2 \\langle \\hat P_X , \\hat P_Y \\rangle_{\\mathcal{H}} \\\\\n", " &= \\sum_{i=1}^{N_X} \\sum_{j=1}^{N_X} K(X_i, X_j) + \\sum_{i=1}^{N_Y} \\sum_{j=1}^{N_Y} K(Y_i, Y_j)\n", " - 2 \\sum_{i=1}^{N_X} \\sum_{j=1}^{N_Y} K(X_i, Y_j).\n", "\\end{aligned}\n", "$$\n", "This is relatively straightforward to compute." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def mmd(x, y, sigma):\n", " # compare kernel MMD paper and code:\n", " # A. Gretton et al.: A kernel two-sample test, JMLR 13 (2012)\n", " # http://www.gatsby.ucl.ac.uk/~gretton/mmd/mmd.htm\n", " # x shape [n, d] y shape [m, d]\n", " # n_perm number of bootstrap permutations to get p-value, pass none to not get p-value\n", " n, d = x.shape\n", " m, d2 = y.shape\n", " assert d == d2\n", " xy = torch.cat([x.detach(), y.detach()], dim=0)\n", " dists = torch.cdist(xy, xy, p=2.0)\n", " # we are a bit sloppy here as we just keep the diagonal and everything twice\n", " # note that sigma should be squared in the RBF to match the Gretton et al heuristic\n", " k = torch.exp((-1/(2*sigma**2)) * dists**2) + torch.eye(n+m)*1e-5\n", " k_x = k[:n, :n]\n", " k_y = k[n:, n:]\n", " k_xy = k[:n, n:]\n", " # The diagonals are always 1 (up to numerical error, this is (3) in Gretton et al.)\n", " # note that their code uses the biased (and differently scaled mmd)\n", " mmd = k_x.sum() / (n * (n - 1)) + k_y.sum() / (m * (m - 1)) - 2 * k_xy.sum() / (n * m)\n", " return mmd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gretton et. al. recommend to set the parameter $\\sigma$ to the median distance between points $\\sigma = \\mathrm{Median(|Z_i - Z_j|)/2}$ where $Z_\\cdot$ is the combined sample of $X_\\cdot$ and $Y_\\cdot$. We stick with this recommendation without further questioning. In practice we might take the median of a reasonable subset (e.g. 100 samples)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor(0.3407)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dists = torch.pdist(torch.cat([x, y], dim=0)[:,None])\n", "sigma = dists[:100].median()/2\n", "our_mmd = mmd(x[:, None], y[:, None], sigma)\n", "our_mmd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Detecting different distributions\n", "\n", "But now we need to make precise when $\\hat P_X$ and $\\hat P_Y$ are fare enough apart to say with some confidence that they $P_X$ and $P_Y$ are distinct.\n", "\n", "Clearly, one option, in particular for evaluation of a detector, is to take the Receiver Operating Characteristic point of view and consider what any threshold for $\\widehat{MMD}$ means in terms of true positives and false positives. This is what we do in our experiments.\n", "\n", "## Bootstrapping\n", "\n", "The other route is to try to give a threshold for a given confidence level or equivalently to convert $\\widehat{MMD}$ values (given the sample sizes) into $p$-values. We could try to derive expressions for this (and earlier papers of the same group of authors do, see the link above), but the conceptually easiest way is to sample from the distribution $\\widehat{MMD}$ under the null-hypothesis using bootstrapping. In this technique, we approximate sampling from the null-hypothesis by shuffling between the $X_\\cdot$ and $Y_\\cdot$, so that both the $x$ and $y$ argument come from the same distribution. If one sample is sufficiently large (e.g. the training sample), we might compute the thresholds just on that, too.\n", "\n", "*Note*: To do this efficiently, it is recommended to use custom CPU or GPU kernels.\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████| 1000/1000 [00:00<00:00, 4415.20it/s]\n" ] } ], "source": [ "N_X = len(x)\n", "N_Y = len(y)\n", "xy = torch.cat([x, y], dim=0)[:, None].double()\n", "\n", "mmds = []\n", "for i in tqdm.tqdm(range(1000)):\n", " xy = xy[torch.randperm(len(xy))]\n", " mmds.append(mmd(xy[:N_X], xy[N_X:], sigma).item())\n", "mmds = torch.tensor(mmds)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([ 73., 199., 198., 155., 124., 62., 66., 37., 26., 18., 11.,\n", " 9., 4., 6., 6., 3., 1., 0., 1., 1.]),\n", " array([0.05134737, 0.07508753, 0.09882769, 0.12256785, 0.146308 ,\n", " 0.17004816, 0.19378832, 0.21752848, 0.24126863, 0.2650088 ,\n", " 0.28874895, 0.3124891 , 0.33622926, 0.35996944, 0.38370958,\n", " 0.40744972, 0.4311899 , 0.45493007, 0.4786702 , 0.50241035,\n", " 0.5261505 ], dtype=float32),\n", " )" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQeUlEQVR4nO3dfYxldX3H8fenoCTVGlEGJAs6QFctNLq2I32wGixtRWlcaFGXGEor6UoLjaaauNhETRsS2oq2aYtkFQJNlIeKKBFsJdRIjI8DrsijLrjKyoYd0YhGQ9312z/mrF6XGe6dOffO7P72/Upu7jm/c86c7+9e+Owv556HVBWSpLb80moXIEkaP8NdkhpkuEtSgwx3SWqQ4S5JDTp4tQsAOOyww2p6enq1y5Ck/cptt932naqaWmjZPhHu09PTzM7OrnYZkrRfSfLNxZZ5WEaSGmS4S1KDDHdJapDhLkkNMtwlqUFDwz3J0Uk+leSeJHcleVPX/owkNyf5evd+6MA2FyTZmuS+JK+YZAckSY83ysh9F/CWqvo14LeB85IcD2wCbqmqtcAt3Tzdsg3ACcApwCVJDppE8ZKkhQ0N96raUVW3d9M/AO4B1gDrgSu71a4ETuum1wNXV9VjVfUNYCtw4pjrliQ9gSUdc08yDbwI+AJwRFXtgPl/AIDDu9XWAA8ObLa9a9v7b21MMptkdm5ubhmlS5IWM/IVqkmeClwHvLmqHk2y6KoLtD3uiSBVtRnYDDAzM7NqTwyZ3nTjsrfddtGpY6xEksZnpHBP8iTmg/2DVfWRrvnhJEdW1Y4kRwI7u/btwNEDmx8FPDSugvcl/sMgaV81ytkyAS4D7qmq9wwsugE4u5s+G/jYQPuGJIckOQZYC3xxfCVLkoYZZeT+EuAs4KtJtnRtbwcuAq5Ncg7wLeA1AFV1V5JrgbuZP9PmvKraPe7CJUmLGxruVfUZFj6ODnDyIttcCFzYoy5JUg9eoSpJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNGuUZqpcn2ZnkzoG2a5Js6V7b9jx+L8l0kh8PLLt0grVLkhYxyjNUrwD+HfjPPQ1V9bo900kuBr4/sP79VbVuTPVJkpZhlGeo3ppkeqFlSQK8Fvj9MdclSeqh7zH3lwIPV9XXB9qOSfLlJJ9O8tLFNkyyMclsktm5ubmeZUiSBvUN9zOBqwbmdwDPrqoXAX8LfCjJ0xbasKo2V9VMVc1MTU31LEOSNGjZ4Z7kYOBPgGv2tFXVY1X1SDd9G3A/8Ny+RUqSlqbPyP0PgHuravuehiRTSQ7qpo8F1gIP9CtRkrRUo5wKeRXwOeB5SbYnOadbtIFfPCQD8DLgjiRfAT4MnFtV3x1nwZKk4UY5W+bMRdr/fIG264Dr+pclSerDK1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0ND7uWsypjfduOxtt1106hgrkdQiR+6S1KBRHrN3eZKdSe4caHtXkm8n2dK9XjWw7IIkW5Pcl+QVkypckrS4UUbuVwCnLND+3qpa171uAkhyPPPPVj2h2+aSPQ/MliStnKHhXlW3AqM+5Ho9cHVVPVZV3wC2Aif2qE+StAx9jrmfn+SO7rDNoV3bGuDBgXW2d22Pk2Rjktkks3Nzcz3KkCTtbbnh/j7gOGAdsAO4uGvPAuvWQn+gqjZX1UxVzUxNTS2zDEnSQpYV7lX1cFXtrqqfAu/n54detgNHD6x6FPBQvxIlSUu1rHBPcuTA7OnAnjNpbgA2JDkkyTHAWuCL/UqUJC3V0IuYklwFnAQclmQ78E7gpCTrmD/ksg14I0BV3ZXkWuBuYBdwXlXtnkjlkqRFDQ33qjpzgebLnmD9C4EL+xQlSerHK1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQUPDPcnlSXYmuXOg7Z+T3JvkjiTXJ3l61z6d5MdJtnSvSydYuyRpEaOM3K8ATtmr7Wbg16vqBcDXgAsGlt1fVeu617njKVOStBRDw72qbgW+u1fbJ6tqVzf7eeCoCdQmSVqmcRxzfwPwiYH5Y5J8Ocmnk7x0sY2SbEwym2R2bm5uDGVIkvY4uM/GSf4O2AV8sGvaATy7qh5J8pvAR5OcUFWP7r1tVW0GNgPMzMxUnzoONNObbuy1/baLTh1TJZL2VcseuSc5G/hj4PVVVQBV9VhVPdJN3wbcDzx3HIVKkka3rHBPcgrwNuDVVfWjgfapJAd108cCa4EHxlGoJGl0Qw/LJLkKOAk4LMl24J3Mnx1zCHBzEoDPd2fGvAz4+yS7gN3AuVX13QX/sCRpYoaGe1WduUDzZYusex1wXd+iJEn9eIWqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWhouCe5PMnOJHcOtD0jyc1Jvt69Hzqw7IIkW5Pcl+QVkypckrS4UUbuVwCn7NW2CbilqtYCt3TzJDke2ACc0G1zyZ4HZkuSVs7QcK+qW4G9H3K9Hriym74SOG2g/eqqeqyqvgFsBU4cT6mSpFEt95j7EVW1A6B7P7xrXwM8OLDe9q7tcZJsTDKbZHZubm6ZZUiSFjLuH1SzQFsttGJVba6qmaqamZqaGnMZknRgW264P5zkSIDufWfXvh04emC9o4CHll+eJGk5lhvuNwBnd9NnAx8baN+Q5JAkxwBrgS/2K1GStFQHD1shyVXAScBhSbYD7wQuAq5Ncg7wLeA1AFV1V5JrgbuBXcB5VbV7QrVLkhYxNNyr6sxFFp28yPoXAhf2KUqS1I9XqEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDhj6JaTFJngdcM9B0LPAO4OnAXwJzXfvbq+qm5e5HkrR0yw73qroPWAeQ5CDg28D1wF8A762qd4+jQEnS0o3rsMzJwP1V9c0x/T1JUg/jCvcNwFUD8+cnuSPJ5UkOXWiDJBuTzCaZnZubW2gVSdIy9Q73JE8GXg38V9f0PuA45g/Z7AAuXmi7qtpcVTNVNTM1NdW3DEnSgGUfcx/wSuD2qnoYYM87QJL3Ax8fwz6e0PSmGye9C0nar4zjsMyZDBySSXLkwLLTgTvHsA9J0hL0Grkn+WXgD4E3DjT/U5J1QAHb9lomSVoBvcK9qn4EPHOvtrN6VSRJ6s0rVCWpQYa7JDXIcJekBo3jVEgdQPqcdrrtolPHWImkJ+LIXZIaZLhLUoMMd0lqkMfcD0DerkFqnyN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqUN/H7G0DfgDsBnZV1UySZwDXANPMP2bvtVX1vX5lSpKWYhwj95dX1bqqmunmNwG3VNVa4JZuXpK0giZxWGY9cGU3fSVw2gT2IUl6An3DvYBPJrktycau7Yiq2gHQvR++0IZJNiaZTTI7NzfXswxJ0qC+d4V8SVU9lORw4OYk9466YVVtBjYDzMzMVM86JEkDeo3cq+qh7n0ncD1wIvBwkiMBuvedfYuUJC3NssM9yVOS/MqeaeCPgDuBG4Czu9XOBj7Wt0hJ0tL0OSxzBHB9kj1/50NV9d9JvgRcm+Qc4FvAa/qXKUlaimWHe1U9ALxwgfZHgJP7FCVJ6scrVCWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJalDfe8tII5vedOOyt9120aljrERqnyN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqUJ9nqB6d5FNJ7klyV5I3de3vSvLtJFu616vGV64kaRR9bj+wC3hLVd3ePSj7tiQ3d8veW1Xv7l+eJGk5+jxDdQewo5v+QZJ7gDXjKkwa1Oe+NOC9aXTgGcsx9yTTwIuAL3RN5ye5I8nlSQ5dZJuNSWaTzM7NzY2jDElSp3e4J3kqcB3w5qp6FHgfcBywjvmR/cULbVdVm6tqpqpmpqam+pYhSRrQK9yTPIn5YP9gVX0EoKoerqrdVfVT4P3Aif3LlCQtRZ+zZQJcBtxTVe8ZaD9yYLXTgTuXX54kaTn6nC3zEuAs4KtJtnRtbwfOTLIOKGAb8MYe+5AkLUOfs2U+A2SBRTctvxxJ0jh4haokNchnqOqA4PNbdaBx5C5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa5KmQ0hCeRqn9kSN3SWqQI3dpghz1a7U4cpekBhnuktQgw12SGuQxd6lBfR8o3oe/FewbHLlLUoMcuUv7qNUcfWv/N7FwT3IK8K/AQcAHquqiSe1L0r7D0z/3DRM5LJPkIOA/gFcCxzP/XNXjJ7EvSdLjTWrkfiKwtaoeAEhyNbAeuHtC+5PUgNUa9bf4A/Skwn0N8ODA/HbgtwZXSLIR2NjN/jDJfROqZdIOA76z2kWsIvtv//eJ/ucfV2W3vfvfs+7nLLZgUuGeBdrqF2aqNgObJ7T/FZNktqpmVruO1WL/7b/93zf7P6lTIbcDRw/MHwU8NKF9SZL2Mqlw/xKwNskxSZ4MbABumNC+JEl7mchhmaraleR84H+YPxXy8qq6axL72gfs94eWerL/Bzb7v49KVQ1fS5K0X/H2A5LUIMNdkhpkuI8gySlJ7kuyNcmmBZY/P8nnkjyW5K2rUeMkjdD/1ye5o3t9NskLV6POSRrhM1jf9X9Lktkkv7cadU7KsP4PrPfiJLuTnLGS9U3aCN//SUm+333/W5K8YzXq/AVV5esJXsz/IHw/cCzwZOArwPF7rXM48GLgQuCtq13zKvT/d4FDu+lXAl9Y7bpX4TN4Kj//DesFwL2rXfdK9n9gvf8FbgLOWO26V/j7Pwn4+GrXOvhy5D7cz26lUFX/B+y5lcLPVNXOqvoS8JPVKHDCRun/Z6vqe93s55m/rqElo3wGP6zu/3LgKex10d5+bmj/O38DXAfsXMniVsCo/d+nGO7DLXQrhTWrVMtqWGr/zwE+MdGKVt5In0GS05PcC9wIvGGFalsJQ/ufZA1wOnDpCta1Ukb9f+B3knwlySeSnLAypS3OcB9u6K0UGjdy/5O8nPlwf9tEK1p5I30GVXV9VT0fOA34h0kXtYJG6f+/AG+rqt2TL2fFjdL/24HnVNULgX8DPjrpooYx3Ic70G+lMFL/k7wA+ACwvqoeWaHaVsqS/huoqluB45IcNunCVsgo/Z8Brk6yDTgDuCTJaStS3eQN7X9VPVpVP+ymbwKetNrfv+E+3IF+K4Wh/U/ybOAjwFlV9bVVqHHSRvkMfjVJuunfYP6Ht1b+kRva/6o6pqqmq2oa+DDw11X10RWvdDJG+f6fNfD9n8h8tq7q9+9j9oaoRW6lkOTcbvmlSZ4FzAJPA36a5M3M/5r+6GrVPS6j9B94B/BM5kdrALtqH71T3nKM+Bn8KfBnSX4C/Bh43cAPrPu1EfvfrBH7fwbwV0l2Mf/9b1jt79/bD0hSgzwsI0kNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg/4fnymLC+eC20AAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyplot.hist(mmds.numpy(), bins=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the empirical distribution, we can compute the threshold for a given probability $p$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor(0.2768)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "torch.quantile(mmds, 0.95) # threshold for 5% significance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or we may assign a $p$-value to our observation:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor(0.0220)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(our_mmd < mmds).float().mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This concludes our short tour of the MMD statistic and the derived two-sample test.\n", "\n", "I hope you enjoyed it, don't hesitate to E-Mail for questions and comments." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.1+" } }, "nbformat": 4, "nbformat_minor": 4 }