Gaussian KDE with periodic boundary conditions

gaussian_kde_angular_metric.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian KDE with periodic boundary conditions\n",
    "\n",
    "We can run a Gaussian KDE with periodic boundary conditions using ``scikit-learn``'s ``KernelDensity`` method using an angular metric. Unfortunately, the built-in angular metric only works in 2D (lon,lat), but we can hack it to work in 1D...\n",
    "\n",
    "First we install ``scikit-learn``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: scikit-learn in /Users/bovy/miniforge3/envs/py39/lib/python3.9/site-packages (1.1.3)\n",
      "Requirement already satisfied: scipy>=1.3.2 in /Users/bovy/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-learn) (1.8.1)\n",
      "Requirement already satisfied: numpy>=1.17.3 in /Users/bovy/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-learn) (1.22.4)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/bovy/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-learn) (3.1.0)\n",
      "Requirement already satisfied: joblib>=1.0.0 in /Users/bovy/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-learn) (1.2.0)\n"
     ]
    }
   ],
   "source": [
    "!pip install scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we make some mock data that is periodic in angle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOsUlEQVR4nO3dX4xc5X3G8ecJdlQXiHyxq2KtHTaVrF5QKQWtHJClyKJNxT+VXnABUkHhxgoiFaiRKoeLoPbKuUFN4hTLApqgUlAVCLWKaYoEFX8kKGvXBIyTyqJELLjyAuKPiyXk9unFHNrNeGZnzvqMz553vx9ptHPOeffM79gzz77zznvOOIkAAN33ubYLAAA0g0AHgEIQ6ABQCAIdAApBoANAIda19cBTU1OZnZ1t6+EBoJMOHjz4bpLpQdtaC/TZ2VnNz8+39fAA0Em2fzVsG0MuAFAIAh0ACkGgA0AhCHQAKASBDgCFGBnotrfYfsb2UdtHbN8xoM0O2x/aPlzdvjOZcgEAw4wzbfG0pG8lOWT7QkkHbT+V5PW+ds8lua75EgEA4xjZQ09yPMmh6v7Hko5Kmpl0YQCAemqNoduelXSppJcGbL7C9iu2n7R9yZDf32l73vb84uJi/WoBAEONfaao7QskPSrpziQf9W0+JOniJCdtXyPpcUlb+/eRZJ+kfZI0NzfHN2sAKML23U/r7Q9OnbF+ZuMGvbDrynNWx1iBbnu9emH+UJLH+rcvDfgkB2z/te2pJO82VyoArE5vf3BKb+6+9oz1s7ueOKd1jDPLxZLul3Q0yT1D2lxUtZPtbdV+32uyUADA8sbpoW+XdLOkV20frtbdJemLkpRkr6QbJN1m+7SkU5JuDF9WCgDn1MhAT/K8JI9os0fSnqaKAgDUx5miAFAIAh0ACkGgA0AhCHQAKASBDgCFINABoBAEOgAUgkAHgEIQ6ABQCAIdAApBoANAIQh0ACgEgQ4AhSDQAaAQBDoAFIJAB4BCEOgAUAgCHQAKQaADQCEIdAAoBIEOAIUg0AGgEAQ6ABSCQAeAQhDoAFAIAh0ACkGgA0AhCHQAKASBDgCFINABoBAEOgAUgkAHgEKMDHTbW2w/Y/uo7SO27xjQxra/b/uY7Z/bvmwy5QIAhlk3RpvTkr6V5JDtCyUdtP1UkteXtLla0tbq9hVJ91Y/AQDnyMgeepLjSQ5V9z+WdFTSTF+z6yU9mJ4XJW20vanxagEAQ9UaQ7c9K+lSSS/1bZqR9NaS5QWdGfoAgAkaO9BtXyDpUUl3Jvmof/OAX8mAfey0PW97fnFxsV6lAIBljRXotterF+YPJXlsQJMFSVuWLG+W9E5/oyT7kswlmZuenl5JvQCAIcaZ5WJJ90s6muSeIc32S7qlmu1yuaQPkxxvsE4AwAjjzHLZLulmSa/aPlytu0vSFyUpyV5JByRdI+mYpE8k3dp4pQCAZY0M9CTPa/AY+dI2kXR7U0UBAOrjTFEAKASBDgCFINABoBAEOgAUgkAHgEIQ6ABQCAIdAApBoANAIQh0ACgEgQ4AhSDQAaAQBDoAFIJAB4BCEOgAUAgCHQAKQaADQCEIdAAoBIEOAIUg0AGgEAQ6ABSCQAeAQhDoAFAIAh0ACkGgA0AhCHQAKASBDgCFINABoBAEOgAUgkAHgEIQ6ABQCAIdAApBoANAIQh0ACjEyEC3/YDtE7ZfG7J9h+0PbR+ubt9pvkwAwCjrxmjzI0l7JD24TJvnklzXSEUAgBUZ2UNP8qyk989BLQCAs9DUGPoVtl+x/aTtS4Y1sr3T9rzt+cXFxYYeGgAgNRPohyRdnOTLkn4g6fFhDZPsSzKXZG56erqBhwYAfOasAz3JR0lOVvcPSFpve+qsKwMA1HLWgW77Ituu7m+r9vne2e4XAFDPyFkuth+WtEPSlO0FSXdLWi9JSfZKukHSbbZPSzol6cYkmVjFAICBRgZ6kptGbN+j3rRGAECLOFMUAApBoANAIQh0ACgEgQ4AhSDQAaAQBDoAFIJAB4BCEOgAUAgCHQAKQaADQCEIdAAoBIEOAIUg0AGgEAQ6ABSCQAeAQhDoAFAIAh0ACkGgA0AhCHQAKASBDgCFINABoBAEOgAUgkAHgEIQ6ABQCAIdAApBoANAIQh0ACgEgQ4AhSDQAaAQBDoAFIJAB4BCEOgAUAgCHQAKMTLQbT9g+4Tt14Zst+3v2z5m++e2L2u+TADAKOP00H8k6apltl8taWt12ynp3rMvCwBQ18hAT/KspPeXaXK9pAfT86KkjbY3NVUgAGA8TYyhz0h6a8nyQrXuDLZ32p63Pb+4uNjAQwMAPtNEoHvAugxqmGRfkrkkc9PT0w08NADgM00E+oKkLUuWN0t6p4H9AgBqaCLQ90u6pZrtcrmkD5Mcb2C/AIAa1o1qYPthSTskTdlekHS3pPWSlGSvpAOSrpF0TNInkm6dVLEAgOFGBnqSm0Zsj6TbG6sIALAinCkKAIUg0AGgEAQ6ABSCQAeAQhDoAFAIAh0ACkGgA0AhCHQAKASBDgCFINABoBAEOgAUgkAHgEIQ6ABQCAIdAApBoANAIQh0ACgEgQ4AhSDQAaAQBDoAFIJAB4BCEOgAUAgCHQAKQaADQCEIdAAoBIEOAIUg0AGgEAQ6ABSCQAeAQqxruwAAKNXMxg2a3fXEwPUv7Lqy8ccj0AFgQoaF9qCQbwJDLgBQCHroADCm7buf1tsfnDpj/czGDS1Uc6axAt32VZK+J+k8Sfcl2d23fYekf5D0H9Wqx5L8ZXNlAkD73v7glN7cfW3bZQw1MtBtnyfph5K+JmlB0su29yd5va/pc0mum0CNAIAxjDOGvk3SsSRvJPlU0iOSrp9sWQCAusYJ9BlJby1ZXqjW9bvC9iu2n7R9yaAd2d5pe972/OLi4grKBQAMM06ge8C69C0fknRxki9L+oGkxwftKMm+JHNJ5qanp2sVCgBY3jiBviBpy5LlzZLeWdogyUdJTlb3D0hab3uqsSoBACONE+gvS9pq+0u2Py/pRkn7lzawfZFtV/e3Vft9r+liAQDDjZzlkuS07W9K+pl60xYfSHLE9jeq7Xsl3SDpNtunJZ2SdGOS/mEZAMAEjTUPvRpGOdC3bu+S+3sk7Wm2NABAHZz6DwCFINABoBAEOgAUgkAHgEIQ6ABQCAIdAApBoANAIQh0ACgEgQ4AheAr6ACgz2r/qrlhCHQA6LPav2pumE4G+nJ/PV/YdWULFQFA+zoZ6MP+em7f/bRmdz1xxnqCHsBa0MlAH2ZYaA8KeQAoDbNcAKAQBDoAFKKoIZdhZjZuYGwdQPHWRKAztg5gLVgTgQ5gbVsrU50JdADFW8lU5y4i0AGsWSX1ziVmuQBAMdZ0D53ZL0BZunpRraas6UCf9OyXpj6IWSsf6AD9hj33h5nZuKGTF9VqypoO9Ekb9kFM3T8YTe0H6JquXvWwLQR6A9b62zygH+8q20Gg17Dck5ReBNaiuq+JYe8q6RQ1g0Cvgbd/KF3dnnXd18RyExF4bZ09Ar1AvN3FSk368xqef5NFoHdY3d4OH6KuPfxxX1sI9AGWC8rVpO4Lknn3a89aOeUdPQT6AG2F26Q/GBp2XHx13+Sttp5yk50BrB4EegtW2wdDXTnBqisGHW/Xh8FK/H8q0ViBbvsqSd+TdJ6k+5Ls7tvuavs1kj6R9PUkhxqutRhdeXHUHaKpO4Vttb0zaOoPT52ZH03+G3fleYXJGRnots+T9ENJX5O0IOll2/uTvL6k2dWStla3r0i6t/qJDqvbc687ha3uENAwTYVZ3fHm5eoZ10qGwRgTxzDj9NC3STqW5A1Jsv2IpOslLQ306yU9mCSSXrS90famJMcbrxitm/R4at1wnnRPvwtj3PTOIY0X6DOS3lqyvKAze9+D2sxI+rVAt71T0s5q8aTtX9aq9v9N+bt6d4W/u1pMSZ0+hjPq/5Ukf7udYgYZo54zjsHfnWBBzev6c0jq/jGsuP6zeK5dPGzDOIHuAeuygjZKsk/SvjEec/mC7Pkkc2e7nzZ1/Ri6Xr/U/WPoev1S949htdU/zhdcLEjasmR5s6R3VtAGADBB4wT6y5K22v6S7c9LulHS/r42+yXd4p7LJX3I+DkAnFsjh1ySnLb9TUk/U2/a4gNJjtj+RrV9r6QD6k1ZPKbetMVbJ1eypAaGbVaBrh9D1+uXun8MXa9f6v4xrKr63ZuYAgDoOr4kGgAKQaADQCE6F+i2r7L9S9vHbO9qu566bD9g+4Tt19quZSVsb7H9jO2jto/YvqPtmuqw/Ru2/9X2K1X9f9F2TStl+zzb/2b7H9uupS7bb9p+1fZh2/Nt17MS1QmUP7H9i+r1cEXrNXVpDL26DMG/a8llCCTd1HcZglXN9lclnVTvzNrfbbueumxvkrQpySHbF0o6KOmPu/J/UF136PwkJ22vl/S8pDuSvNhyabXZ/jNJc5K+kOS6tuupw/abkuaSdPakIts/lvRckvuqGYC/meSDNmvqWg/9/y5DkORTSZ9dhqAzkjwr6f2261ipJMc/u/Bako8lHVXvrOBOSM/JanF9detOr6Zie7OkayXd13Yta5HtL0j6qqT7JSnJp22HudS9QB92iQG0wPaspEslvdRyKbVUQxWHJZ2Q9FSSTtVf+StJfy7pf1quY6Ui6Z9tH6wuCdI1vy1pUdLfVMNe99k+v+2iuhboY11iAJNn+wJJj0q6M8lHbddTR5L/TvJ76p3RvM12p4a+bF8n6USSg23Xcha2J7lMvSu13l4NRXbJOkmXSbo3yaWS/ktS65/pdS3QucTAKlCNPT8q6aEkj7Vdz0pVb5H/RdJV7VZS23ZJf1SNQz8i6Urbf9tuSfUkeaf6eULST9UbTu2SBUkLS97d/US9gG9V1wJ9nMsQYIKqDxXvl3Q0yT1t11OX7WnbG6v7GyT9gaRftFpUTUm+nWRzkln1XgNPJ/mTlssam+3zqw/UVQ1T/KGkTs36SvKfkt6y/TvVqt/Xr19SvBWd+gq6YZchaLmsWmw/LGmHpCnbC5LuTnJ/u1XVsl3SzZJercahJemuJAfaK6mWTZJ+XM2Y+pykv0/SuWl/Hfdbkn7a6xtonaS/S/JP7Za0In8q6aGqc/mGJn/Jk5E6NW0RADBc14ZcAABDEOgAUAgCHQAKQaADQCEIdAAoBIEOAIUg0AGgEP8L1XiUFgfw9IkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy\n",
    "from sklearn.neighbors import KernelDensity\n",
    "%pylab inline\n",
    "xs= 2*numpy.random.uniform(size=1001)-1\n",
    "xs= xs**3.\n",
    "thetas= (2*numpy.arccos(xs)-numpy.pi-0.2) % (2*numpy.pi)\n",
    "hist(thetas,bins=51,density=True,histtype='step');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute the regular Gaussian KDE as well as one with an angular metric. Note that for the angular metric, we have to augment the data with a second, latitude, dimension that we just set to zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "kde_g = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(numpy.atleast_2d(thetas).T)\n",
    "\n",
    "thetas2D= numpy.tile(thetas,(2,1)).T\n",
    "thetas2D[:,1]= 0.\n",
    "kde_a = KernelDensity(kernel='gaussian', bandwidth=0.2,metric='haversine').fit(thetas2D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compare the result from both to the histogram of the data; again, we have to augment the input data for the angular-distance KDE. It also does not come out to be normalized, so we normalize it by hand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x16d1128e0>]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAh2klEQVR4nO3deZycVZ3v8c+v1t6XdHcWOp10ErJAImuEQFhiQISAwgAqDIOj9/Uargw6OnpV1FHH0bmTK+KGCoMgKqAQZBEkgIEQCDiBJKwJnY0spLP1vi/VVXXuH1WJWbrTW3VX6uH7fr3qla7nOfXUrzrd337q1HnOMeccIiKS+XzpLkBERFJDgS4i4hEKdBERj1Cgi4h4hAJdRMQjAul64tLSUldZWZmupxcRyUhr166tc86V9bYvbYFeWVnJmjVr0vX0IiIZycx29LVPXS4iIh6hQBcR8QgFuoiIRyjQRUQ8QoEuIuIR/Qa6mVWY2fNmVmVm683sC720WWBmzWb2RvL27ZEpV0RE+jKQYYtR4MvOudfMLB9Ya2bLnHPvHNZupXPustSXKCIiA9HvGbpzbo9z7rXk161AFVA+0oWJiMjgDKoP3cwqgVOBV3rZfZaZvWlmT5nZ7D4ef4OZrTGzNbW1tYOvVkRE+jTgK0XNLA94GPiic67lsN2vAZOdc21mtgh4DJh++DGcc3cCdwLMnTtXK2uIiCfMX7ycXU2dR2wvL8rm5ZsXjlodAwp0MwuSCPP7nXOPHL7/4IB3zi01s1+aWalzri51pYqIHJt2NXWyffGlR2yvvPnJUa1jIKNcDLgbqHLO/aiPNuOT7TCzM5LHrU9loSIicnQDOUOfD1wPvG1mbyS3fQOYBOCcuwO4GrjRzKJAJ3CN02KlIiKjqt9Ad869BFg/bX4O/DxVRYmIyODpSlEREY9QoIuIeIQCXUTEIxToIiIeoUAXEfEIBbqIiEco0EVEPEKBLiLiEQp0ERGPUKCLiHiEAl1ExCMU6CIiHqFAFxHxCAW6iIhHKNBFRDxCgS4i4hEKdBERj1Cgi4h4hAJdRMQjFOgiIh6hQBcR8QgFuoiIRyjQRUQ8QoEuIuIRCnQREY9QoIuIeIQCXUTEIxToIiIeoUAXEfEIBbqIiEco0EVEPEKBLiLiEf0GuplVmNnzZlZlZuvN7Au9tDEz+5mZbTGzt8zstJEpV0RE+hIYQJso8GXn3Gtmlg+sNbNlzrl3DmpzCTA9eTsTuD35r4iIjJJ+z9Cdc3ucc68lv24FqoDyw5pdDvzOJawCisxsQsqrFRGRPg2qD93MKoFTgVcO21UO7DzofjVHhr6IiIygAQe6meUBDwNfdM61HL67l4e4Xo5xg5mtMbM1tbW1g6tURESOakCBbmZBEmF+v3PukV6aVAMVB92fCOw+vJFz7k7n3Fzn3NyysrKh1CsiIn0YyCgXA+4GqpxzP+qj2ePAp5KjXeYBzc65PSmsU0RE+jGQUS7zgeuBt83sjeS2bwCTAJxzdwBLgUXAFqAD+EzKKxURkaPqN9Cdcy/Rex/5wW0ccFOqihIRkcHTlaIiIh6hQBcR8QgFuoiIRyjQRUQ8QoEuIuIRCnQREY9QoIuIeIQCXUTEIxToIiIeoUAXEfEIBbqIiEco0EVEPEKBLiLiEQp0ERGPUKCLiHiEAl1ExCMU6CIiHqFAFxHxCAW6iIhHKNBFRDxCgS4i4hEKdBERj1Cgi4h4hAJdRMQjFOgiIh6hQBcR8QgFuoiIRyjQRUQ8QoEuIuIRCnQREY9QoIuIeIQCXUTEIxToIiIe0W+gm9mvzazGzNb1sX+BmTWb2RvJ27dTX6aIiPQnMIA2vwF+DvzuKG1WOucuS0lFIiIyJP2eoTvnXgQaRqEWEREZhlT1oZ9lZm+a2VNmNruvRmZ2g5mtMbM1tbW1KXpqERGB1AT6a8Bk59zJwG3AY301dM7d6Zyb65ybW1ZWloKnFhGR/YYd6M65FudcW/LrpUDQzEqHXZmIiAzKsAPdzMabmSW/PiN5zPrhHldERAan31EuZvYHYAFQambVwHeAIIBz7g7gauBGM4sCncA1zjk3YhWLiEiv+g1059y1/ez/OYlhjSIikka6UlRExCMU6CIiHqFAFxHxCAW6iIhHKNBFRDxCgS4i4hEKdBERj1Cgi4h4hAJdRMQjFOgiIh6hQBcR8QgFuoiIRyjQRUQ8QoEuIuIRCnQREY9QoIuIeIQCXUTEIxToIiIeoUAXEfEIBbqIiEco0EVEPEKBLiLiEQp0ERGPUKCLiHiEAl1ExCMU6CIiHqFAFxHxCAW6iIhHKNBFRDxCgS4i4hEKdBERj1Cgi4h4hAJdRMQj+g10M/u1mdWY2bo+9puZ/czMtpjZW2Z2WurLFBGR/gzkDP03wMVH2X8JMD15uwG4ffhliYjIYPUb6M65F4GGozS5HPidS1gFFJnZhFQVKCIiA5OKPvRyYOdB96uT245gZjeY2RozW1NbW5uCpxYRkf1SEejWyzbXW0Pn3J3OubnOubllZWUpeGoREdkvFYFeDVQcdH8isDsFxxURkUFIRaA/DnwqOdplHtDsnNuTguOKiMggBPprYGZ/ABYApWZWDXwHCAI45+4AlgKLgC1AB/CZkSpWRET61m+gO+eu7We/A25KWUUiIjIkulJURCQFWiItPLr5Ueo769NWQ79n6CIicnT+cDXXPHYlOzv3MSFnPL+88HaOLz5+1OvQGbqIyFB1NvHokispqryNnZ37CMfj7OnYy/VPXc+qPatGvRwFuojIEP3xmc/z7c7NRHzGVWNO4bmis/lwewdtPW3cuOxG/HnvjGo9CnQRkSGItuzhV/VrAejaexn//tF7Kbz4Fn7Y4ecfollEXZSssUuJu/io1aRAFxEZgmeXf53dAT+TcibQ03h2YmN2Eb6P/F++vHMTxwXy8YXreLH6xVGrSYEuIjJIrnk3v6tJ9JFf/4H/xSFROucqAtMW8vf1+wC49517R60uBbqIyCC9seI7vB0OUhjM52PTPnboTjO49FaubGwgEPfx6t5X2dCwYVTqUqCLiAxG8y5+u3sFAJ+YdQ05wZwj24yZSv7EM/hwS6L/fLTO0hXoIiKDsHPNf7M8O0zQF+DaWUe5kP6Ey/h86y4M46ltT1HXWTfitSnQRUQG4cFtS3FmLJpyKWU5R5kGfOYiKqIxFuZV0hPv4YEND4x4bQp0EZEBcvVbWUY7AFfPuProjUumsSFewXUtifZ/3vpnElNfjRwFuojIAG1887fsDgYoDRdzUtlJ/bZ/Jj6X0957jZJwMbvado34h6MKdBGRAXpu+zMAfGjyhfis//j8S2wufhfngtzJADz73rMjWp8CXURkIFp281ysEYCFkxYO6CHrXSUUVnBBSxMAz+5QoIuIpN3ON+9ncyhEXiCHM8efOcBHGcy6lA9ue5WCUAFbm7eytWnriNWoQBcRGYDl7z4BwLkV5xP0Bwf+wFmXEox2sSB/GgDPvffcSJQHKNBFRPrXXsdz3XsBuGDSBYN77KSzIZTPhyOJES7LdixLdXUHKNBFRPpRt+4h3giHCPmCnFN+zuAe7A/ApDM5a89GcgI5VDVUsatt14jUmXGB3tDVwGNbHhuVq65ERACe3/InnBnzJpxFbjB38AeYPJ9w7UbOS/a9j9SHoxkX6P/+13/nWy9/i+XvLU93KSLyfhCNsKJtGzDw0S1HmDwfgAvC44CR60fPuEBfULEAgBeqX0hvISLyvtC1/UVeDSWWXz6/4vyhHeS4UyGQzXnNjcwaM4szxp8xIleNZtwi0edNPA+AVbtX0dHT0ftMZyIiKbL6nSV0+XycWDyL0uzSoR0kEIKKD5KzcxUPffal1BZ4kIw7Qy/NLuWk0pOIxCNpWYRVRN5fVu5bDSSGKw7L5HNg7zrobBp+UX3IuECHv73tUbeLiIwk1/geK60TgHMnnju8g00+G3Dw3sidiGZmoE9MBvrOF0Z1AVYReX/Z9s4SqoNBikMFzCmZM7yDTZwL/hDseDk1xfUiIwN9RvEMjss9jvquetbVrUt3OSLiUSu3JSbjmj/xPPw+//AOFsyG8tMV6IczswPdLit2rkhrLSLiUdFuVra/B8C55cPsbtlv8tmw+w3obkvN8Q6TkYEOsGDiAgBWVK9Iax0i4k3t7y5nbTiID+Ps484e0jHKi7KpvPnJA7frnwuCi/GvP7wzxdUmZNywxf3mjp9LbjCXzY2b2dW2i/K88nSXJCIesqrqQaJmnFI6h6KsoiEd4+WbD7sQqfs8WHwLUzvfHH6BvcjYM/SQP3Tgr6a6XUQkpZxjZc1aAM6r+FDqjhvOh8t+xNOxM1J3zINk7Bk6JC7DXbZjGct2LOO6E65Ldzki4hGxvW+zwh8D/IcMV5y/eDm7mjqPaF9elD3wg5/+adY/9GQKqjzSgALdzC4Gfgr4gbucc4sP278A+BOwLbnpEefcf6SuzN4tmLiAkC/Ea/teo7aj9ugrcIuIDNCbb/2O+oCf8pzxzCyeeWD7rqZOti++NI2VHV2/XS5m5gd+AVwCnAhca2Yn9tJ0pXPulORtxMMcIC+Ux9nlZ+NwI75Wn4i8fzyXHGxxQeVFmFl6ixmEgfShnwFscc5tdc5FgAeAy0e2rIG7aPJFADyTXLxVRGQ4XNNOnnPtwBAWs0izgQR6ObDzoPvVyW2HO8vM3jSzp8xsdm8HMrMbzGyNma2pra0dQrlHWlCxgKAveKDbRURkODa9dS+7ggHGhAo5uezkdJczKAMJ9N7ebxw+7+NrwGTn3MnAbcBjvR3IOXenc26uc25uWVlq+rvzQ/nML5+vbhcRSYlntz0NwIcmXzj8q0NH2UACvRqoOOj+RGD3wQ2ccy3Oubbk10uBoJkNcZ7Jwdvf7fKX7X8ZracUES/qbOK5SA0AF06+MM3FDN5AAn01MN3MpphZCLgGePzgBmY23pKfHJjZGcnj1qe62L7s73ZZu2+tlqYTkSHbuX4Jm0NB8vzZnJlcLi6T9Bvozrko8DngGaAKWOKcW29mnzWzzyabXQ2sM7M3gZ8B17iRWI6jD/mhfOYfl+h2GckVtUXE257b8BCQmPs86A+muZrBG9CVos65pc65Gc65ac65/0xuu8M5d0fy658752Y75052zs1zzv11JIvuzcVTLgbg0c2PjsjSTiLicW01/KUjMRnXBRnY3QIZfOn/4S6cfCFF4SKqGqo0pa6IDNrm1bfzdjhEbiA7dbMrjjLPBHrYH+aK468AYMmmJektRkQyi3M8vOUxAC6d+tGMXavYM4EO8PEZHwfg6W1P09zdnOZqRCRTdL+3iif83QBcNeOqNFczdJ4K9EkFkzhrwll0xbp4/N3H+3+AiAiwbPVPaPH7ObF4JieW9DazSWbwVKADfHLmJwFYsnGJPhwVkf5FOni4ITE/+VUzP5HmYobHc4F+fsX5jM0ey/aW7azeuzrd5YjIMW7b6/ewJhwk2xdi0ZRF6S5nWDwX6AFf4EAf2G/W/ya9xYjIsS0e55G37wHgkqmXkhfKS3NBw+O5QIdEt0tOIIeVu1bqLF1E+tT81u95xNcBwFUzrk5zNcPnyUAvyS7h03M+DcBP1v5EfekicqRYD3e/egstfj9njPsgHyj9QLorGraMXoLuaP7xxH/kwQ0P8lbdWyzbsYyLKi9Kd0kicgzZ/covuT8YA4wvffDLhyxkkZKl5tLAs4GeE8zhxpNv5PuvfJ+fvf4zPjTpQwR9mTc3g4iMgJ5Obnv7TiJZPhZNWcTskkOXcDjWl5rrS8YGem9/QcuLsnn55oUH7l8540rurbqXHS07eHjTw1wz65rRLlNEjkFVL/4nf87yEbQA/3Lav6S7nJTJ2EDv7S/o/MXLqbz50NW0A/nnkD1xBz9e+2PmTZhHZWHlKFYpIseaWPVqbtnyEGSF+PsTrqM8r7cF2DJTxgZ6bw4+O9/PuUWccNvbdBS+xVde/Ar3LbqPsD+chupEJO3a6/j5E59hdU6IolAB/3TSP6W7opTy5CiXg5kZXXuvpCK/gg0NG/jh6h+muyQRSYdYlKV//AR35Rh+83HLglspDBemu6qU8nygAxDP4ofn/5CgL8gDGx/gme3PpLsiERlN8Tjr//zPfDu+D4CvfvBrzJswL81Fpd77ItDLi7JZdMs2WncnFsH48vNfY9r3bmX+4uVprkxERlx3K2/84e+4qe4lun0+rpp+FdfOujbdVY0IT/Wh92V/37pzi/ivV/P4w4Y/UFB5H3t3xIEj+91FxBtc/VYeePjj/CDUTTTgZ96EeXzzzG8eMubcS7wX6LEe2L4SNi+DthroboWeDig4Dhs3h6+PO5PATMe9Gx8ge+J9PL39A1xceXG6qxaRVOpqYe/KH/DjLUtYmhMGjEj9OSyruoTpy/9yoNnhQ50znXcCvXE7rLwVqp6AzkYIZEPBBAjlQTAHtr0Ibz2IAV8JZBOYMod7ovv4ygtf4fV9r/OluV/S6BeRTNfRQN2aX3H3untYkh0gkhMm2x+m4b2/Y8s3vnlE896GOsOxf0VoXzI/0CMd8NKP4eWfgs8PJ3wUTrwcpi2E4GH/Ke31sOcNbMOT/Os7f6I00M2PxxTx+w2/Z83uv/KDD/2EaUXT0vM6RGRoerroeXc5L7/+3zzRuJ4V2VlEchNXhX+k8iN87pTPseC/1vf6UC+dnUOmB/r2l+CR/w0t1fCBj8OF34XCo1wkkFsCx18Ax1+AXfIDnv3WD7hv/Gq+2lHFppbtXP2nK7hq3Nl89tzvUZo7dvReh4gMTuN2Gjc9xcubH+elli28nBWkye+H3GwMY8HE87np1JuYNWZW8gG9B7rXZGagO8c/+f8Mv30QxkyBzzwFk88e3DH8Ad4tmMelb53MeN8eZo2/n00FtTy476/8aclCriuazdWnf56Jk84ZmdcgIgMTj+Fqqti79TnerF7J2saNrPVF2BwKJfbnZgEwrXAq1TtnU79vNk9UFfHEsneBd4HM7UIZLEvX1LJz5851a9asGfwDu1rgT/+c6Cs/4WNw+S8gqyAlNW2p38BPX/oWK5o2HNg2LxbgY+PO5KyZV1JauQACoQEf72gztg3mrV6qjiNyzIt242o3sXvnX9mwexUrdrxFXbCVqnCA+oD/0LZxP9GOKUTbZxJtm4WLlFJelJMRvxOVNz855Mm/zGytc25ub/sy7wx941LYsJTv91zHv33iF5DC4UfHl8zitssfYup3f8knz93EsuoVrPJHWVX3MtS9zIzno8wNFjOjoJLpZScxrXweuaUzIG88+I4c0t/XjG29fQhzNKk6jsgxo6OBzpp32LlnLTvq1rGjZTvbOmvZGu9kazBAx/7fpyKAxGCFwmA+c8pO4rRxp3H6uNOZUzpHAxkOk3mBftIn4bhTuevWLfzbCI0ljXdO5r8W/jM3dzezdNtSXtj+LGtrX2dTCDbRCq1vJ25b76cgFmNCNE5R1Ed2NIQ/mgPRPLqjY/jomONo3BKmqOwELH984kNbkfeB+YuXU9vUwiTfLiaENlMYfo9AqJZYuJXWcDc7fVATOCx+ggCJd8BjQgXMLDmRF9aF+PEVlzKndA4V+RWeHT+eKpkX6GZQNhPYMuJPVRgu5NpZ13LtrGvpjnVz3k/voTH2Lr7wHsLhXbhQAy1+aPH7kycRMaA1edsDrOe8l5cRjscZH4szngDjA7lcOi7Aw48/TFnecZQUVlBaOI2CguPIyhmLZReCX/O2S2aYv3g5LU01TA5uoSxrKzmh3cRDDRxf1E5eWZw9fj/7jgjhROwEzEd5zjgqi47nfzb6aWktJt49lnhkLK2xXHaQvMp76rHfhXKsyLxAT5OwP8y+mgq2L/7sgW3OORq6Gtjbvpc97Xuo7ayltiNxq2nbRU37HvZ11tFKFzt8PnYA0A5j4MXG1dAI7PzbcwScIy8eJ8s5Qs4IAX6MU6c6rr37a/gxAhgB85FlfhZOjPMf991CYTCX0qxSSnLHMa5oChPHn0bpuFOwoN6OyuAc+XmNI59WTiiq5ZsXZbO3cSt7W3eyu30vuyKNjCnppHucHfTzvZ8P8BHAx8ScsUwumsakoqk89D+dNDQXEo+U4HoKacTPOhLBvTYD+r6PdQr0YTAzSrJLKMkuYXbp7D7btfe0Hwj9ve17+cYTK/nEGbnUte+hvr2Guu5GWqOddBNNDL06KkfinUAM8mF1rB5i9dD1HjQBu4D1kBWPM74HiiI5+LtKaeueTGfgTJ782mdS9vrFG5qbdlC96xV21q5neuBVLjk5xr5IM/tindQQo8FnVJnxD28f9kAfEEr8vJYE85icP4nK4hlUFk1lSuEUKgsqKc8vP2SlsK9+cPRe1/uRAn0U5AZzmVY07cBFS/+nNpvvn3vkh5yRWITWSCvdsW4isQiReIRYPMalt73IYzedRczFiMaj9MR66Ip18bVHV9PU1Yr5OvAFWggHmggG63GhFrr8EbaHgXAX5FcD1cDLLPz1rZwYKGB20TTmlJ/NnOmXUVxQMarfDxk9B59xB6yDyVnrKMveSHbOHiy/g3fjnTT4DuoSKYO1keTXfgAfBoSiYaYXT2Bc3nGMK6igPL+C8rxybvj1Nt765rXkhfJG+ZVJbxTovSgvyk7L5cAhf4iS7JIjtse7tnJS2UlHbF/9hb7forZEWtjevJ3NjZvZUl/FhprXebNhK7X+KC+4Vl5ofAMa34B1v2RsD5ycU8rs4lnMKj+TmZUfprTAO6u4vB81d9SxcdtzTPD/kfNOaaaqs4ZtFqXGjJqDG/oMf9xHKJJPrKcEn03kiwvPZlzOOMbmjGVczjhKsks4/hvPsLaXkVbx7icV5scQBXov0jWONZUrjReECjip7KRD/hDEXZydrTtZX7eedXvXsH7vGqra3qMmGGNZTx3Lal6Cmpfg9VspjDpKo1mcXjqJyuJpHDdmJhNKZzO+ZAZFWcX47H0x83LKJP5vO4B4cosPsGFdS9DR08Hepq3srX2H6rr1bGvYxLb2Xbzb08xeSz7PWNjYnXg6Pz6OzypjVskJzBx3GsePmcHUwqmMzx3f7/9nuk5yZHAU6GlwtF+OkVxp3Gc+JhdMZnLBZBZNXQRANB5lW/M21tet553dq9hYt45N7btpDvTQHOjm3bbN0LYZdj79t+M4KMAosgDZFiDLFyTsCxA0P358+M2H78DIBuO9+g6iccBZ4oYfnJ+QP8y5x5eTHcwjN5RPfk4pBTljKcibQHHBRIqyiikIFRwTQ9XiLk5rpJXGznqaWqtpbNlFU9seWjpqae5qoLm7ieZIK1saGum2CD2+GD0WJ+JzRMwRHQ/5Ew59Heag0znO+o0RwgjjI2w+wuYnaL4D30uAmIsTJU5XPEar66GVOJ19fVsMsvAzPWccm3cX85WP/B0nls5mevH0IY/bzoSLdWSAgW5mFwM/JdGrdpdzbvFh+y25fxHQAXzaOfdaimv1jGPplyPgCzC9eDrTi6dzxfQrgMTond3tu7nyV4/RGtlOYXgngWA9LthKV6CTbn+cJhxN9IDrgVhn4nPavuT2vatq3/aj1+dgjPkp8WUxJphHSbiQkuSInuK8CRTlT6S4oIL8rCLyQnnkh/IJ+UJ9/hFwzhGNR2nvaae1p5XW7haa2/bQ1FLN7c+/SiRaRyDQBv52Yv4uIv4InYEY7T5HfCB/V3o9Yf3bA/0kPtbe//etx4yeRGUc+LDb9STu9iV5OH/cyImG8ffkEqSU6+YtYErpiUwpnMKk/En4fX4qb36ST8wauZMEObb0G+hm5gd+AXyYxCdrq83scefcOwc1uwSYnrydCdye/FcykJlRnlfOK/96U6/7K29+nLXfmU9zdzOd0U66Yl10Rbv41D3/w68+dRoxF8M5h0umknOOuIsTJ040GiEa66anp5NItJ2O7lY6I620R1po7W6mJdLKprpaunzddPt6iPgdNcSoibdDdzt074OWTUev30EYI2SGD2N/JT04ut1Rgrno6N+XcMxHOBYgEAtjsSxcLJdoLI9IrICuWCGdsQLGZBdx53XnkBvMJSeYQ3Ygmyx/FkF/kIAFDvyh2f896Yn3EIlHiMQidMe6E7doN5/5zSpq2joxi4EDhx+cj7F5eTx244UUhArIDmQfON78xcv53pJOoBN4J3lTl8j7zUDO0M8AtjjntgKY2QPA5ez/iUm4HPidS0wMs8rMisxsgnNuT8orlrQrL8rj9O+u6mX7aSyclNp3H92xbho6G6jvqqe+s46GlmrqW3ZS37abxq56mrqb2dxYTzcRIv4Y3b44cYMuHF29neZa4oc+Bz/5vhD5/jD5gRyKw4W8tq2Hq06bTXHeBMYUVFCcO5bicDFjssZQlFV0yPC74TIz/ObH7/OTRdYR+1d9+YRBHe9Yetcn6TOQQC/nkMtfqObIs+/e2pSTuFzyADO7AbghebfNzDYOqtq/KQXq7P8N8dHHhlKgLt1FDMMR9e8A7OvpKWaIjngNK3gsPZUMTab/DEHmv4Yh1z+M/Jrc146BBHpvb1APP/UZSBucc3cCdw7gOY9ekNmavmYbyxSZ/hoyvX7I/NeQ6fVD5r+GY63+gYw9qwYOvvJkIrB7CG1ERGQEDSTQVwPTzWyKmYWAa4DHD2vzOPApS5gHNKv/XERkdPXb5eKci5rZ54BnSIy6+rVzbr2ZfTa5/w5gKYkhi1tIDFsc6QlDht1tcwzI9NeQ6fVD5r+GTK8fMv81HFP1p23FIhERSS1dvy0i4hEKdBERj8i4QDezi81so5ltMbOb013PYJnZr82sxszWpbuWoTCzCjN73syqzGy9mX0h3TUNhpllmdmrZvZmsv7vprumoTAzv5m9bmZ/TnctQ2Fm283sbTN7w8yGsFp8+iUvoPyjmW1I/j6clfaaMqkPPTkNwSYOmoYAuPawaQiOaWZ2HtBG4sraOemuZ7DMbAIwwTn3mpnlA2uBKzLl/yA571Cuc67NzILAS8AXnHNHXvp6DDOzLwFzgQLn3GXprmewzGw7MNc5l7EXFZnZb4GVzrm7kiMAc5xzTemsKdPO0A9MQ+CciwD7pyHIGM65F4GGdNcxVM65PfsnXnPOtQJVJK4KzgguoS15N5i8Zc5ZDWBmE4FLgbvSXcv7lZkVAOcBdwM45yLpDnPIvEDva4oBSQMzqwROBV5JcymDkuyueAOoAZY55zKqfuAnwFf52+TqmcgBfzGztckpQTLNVKAWuCfZ9XWXmR1lXtHRkWmBPqApBmTkmVke8DDwRedcS7rrGQznXMw5dwqJK5rPMLOM6foys8uAGufc2nTXMkzznXOnkZip9aZkV2QmCQCnAbc7504F2oG0f6aXaYGuKQaOAcm+54eB+51zj6S7nqFKvkVeAVyc3koGZT7wsWQf9APAQjO7L70lDZ5zbnfy3xrgURLdqZmkGqg+6N3dH0kEfFplWqAPZBoCGUHJDxXvBqqccz9Kdz2DZWZlZlaU/DobuBDYkNaiBsE593Xn3ETnXCWJn//lzrl/SHNZg2JmuckP1El2U1wEZNSoL+fcXmCnmc1MbrqAQ6cUT4uMWoKur2kI0lzWoJjZH4AFQKmZVQPfcc7dnd6qBmU+cD3wdrIfGuAbzrml6StpUCYAv02OmPIBS5xzGTn0L4ONAx5NLs4RAH7vnHv66A85Jn0euD95crmVkZ/ypF8ZNWxRRET6lmldLiIi0gcFuoiIRyjQRUQ8QoEuIuIRCnQREY9QoIuIeIQCXUTEI/4/EsfxphrE8xAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist(thetas,bins=51,density=True,histtype='step');\n",
    "vs= numpy.linspace(0,2*numpy.pi,100)\n",
    "plot(vs,numpy.exp(kde_g.score_samples(numpy.atleast_2d(vs).T)))\n",
    "vs_a= numpy.tile(vs,(2,1)).T\n",
    "vs_a[:,1]= 0.\n",
    "yvs_a= numpy.exp(kde_a.score_samples(vs_a))\n",
    "yvs_a/= numpy.sum(yvs_a)*(vs[1]-vs[0])\n",
    "plot(vs,yvs_a,lw=2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the Gaussian KDE with angular metric is indeed periodic."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('py39')",
   "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.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "a82f9eb2f595400cd0f024a877274681853baaeedd921355b08c40b0424842e5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}