init change of variables

This commit is contained in:
ritchie46 committed 2019-09-22 14:12:36 +02:00
1 parent 2d017055c3
commit 057084bc79
3 files changed
+351

No files matched your search

@@ -0,0 +1,348 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import torch\n",
"import torch.distributions as dist\n",
"import numpy as np\n",
"from scipy import stats\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Transforming probability distributions\n",
"Normalizing flows uses transformation to obtain complex posteriors. \n",
"This notebook explores transformations of distributions.\n",
"\n",
"We start with a basis Gaussian $\\mathcal{N}(\\mu=1, \\sigma=0.1)$. These parameters are chose so we can look at transformations in the positive domain and don't have trouble with square rooting negative numbers.\n",
"\n",
"Below we plot the **pdf** and integrate to validate $\\int P(x)dx=1$."
]
},
{
"cell_type": "code",
"execution_count": 221,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 221,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3Bc53nf8e+D+424L0mQAAhSohyRtm6GZF3sWHZtVZYvmsaeWprUStS0jBy7dS5NJ00z9rhtOu1MmjiyUrOK4iaqHTue2rEVl4qtxJItR9EFokRKFCURIgkCIgWAAIg7QFye/rG71BoCiANgd8/Z3d9nZoe7ey77aHXww4v3vOc95u6IiEjuKwq7ABERSQ8FuohInlCgi4jkCQW6iEieUKCLiOSJkrA+uLm52Ts6OsL6eBGRnPTss8+edffYcstCC/SOjg66urrC+ngRkZxkZj0rLVOXi4hInlCgi4jkCQW6iEieUKCLiOSJwIFuZsVm9pyZfX+ZZWZm95pZt5kdNrNr0lumiIisZi0t9M8BR1dY9iFgd+KxD/jKBusSEZE1ChToZtYKfBh4YIVVbgce9LgngXoza0lTjSIiEkDQFvqXgH8PLK6wfDvQm/K6L/HezzCzfWbWZWZdg4ODaypUJNMO953jb198A00pLblq1UA3s48AA+7+7MVWW+a9t/xUuPv97t7p7p2x2LIXOomEYmZugV/9P89yz9ee5a6vPk3v8FTYJYmsWZAW+k3Ax8zsJPBN4P1m9rUl6/QBbSmvW4HTaalQJAu+/tQpzozO8Es37OBgzwi3/NFP+NYzvatvKBIhqwa6u/8Hd2919w7gDuBH7v4vlqz2EHBXYrTL9cCou59Jf7ki6Tc5O8//fLSbmy5t4ou3v50f/uZ7ecf2On7vey8yOTsfdnkiga17HLqZ3WNm9yReHgCOA93AnwK/lobaRLLif//DCYYmz/PvbnkbANvrK/nNWy7j/Pwij72icz2SO9Y0OZe7PwY8lni+P+V9Bz6TzsJEsmF0ao7/9ZPjfODyLVzd3nDh/Ws7GmmqLuNvj7zBh6/QgC3JDbpSVAranz5+nInZeX7rlst+5v3iIuODe7bw6MsDzM4vhFSdyNoo0KWgfe/Q67zvbZu5vKX2Lcv+6du3MjE7zxPdQyFUJrJ2CnQpWG+MztA7PM2NlzQtu/zGS5qoKS/hb198I8uViayPAl0KVlfPMADX7Wxcdnl5STHv/7nNPHK0n/mFla6pE4kOBboUrK6TI1SVFbNnme6WpFvfvpXhyfN09YxksTKR9VGgS8F65uQwV7fXU1K88o/Bey+LUV5SpG4XyQkKdClI4zNzHD0zRueO5btbkqrLS/j5y2L88IjmeJHoU6BLQXru1DkWPT7efDU3XtLE6dEZBsdns1CZyPop0KUgdZ0cprjIuKq9ftV1k33sR06PZboskQ1RoEtBeubkCHtaaqkpX/1i6T3bkoE+mumyRDZEgS4FZ25hked6R+jsaFh9ZWBTRSk7mqp46Yxa6BJtCnQpOEdOjzEztxio/zxp77ZadblI5CnQpeB0nYxfUNS5I1gLHeL96D1DU4zNzGWqLJENU6BLwXnm5DA7mqrYXFsReJu92+oAOKpWukSYAl0KzuG+Ua5pD946h3iXC6B+dIk0BboUlPGZOc6MznDZlk1r2i62qZzmmjL1o0ukBblJdIWZPW1mh8zsiJl9cZl1bjazUTN7PvH4fGbKFdmY7oEJAC7dXLOm7cyMPdvqFOgSaUHuWDQLvN/dJ8ysFPipmT3s7k8uWe9xd/9I+ksUSZ9jiUDfvcZAh3i3ywOPH+f8/CJlJfrjVqInyE2i3d0nEi9LEw9NaiE5qXtggrKSItoaq9a87Z6WWuYWnFf7xzNQmcjGBWpmmFmxmT0PDACPuPtTy6x2Q6Jb5mEz27vCfvaZWZeZdQ0O6ua7kn3dAxNcEquhuMjWvO2FE6PqdpGIChTo7r7g7lcBrcB1Zvb2JascBHa4+5XAl4HvrrCf+9290907Y7HYRuoWWZdjA+Pr6m4B6GiqpqqsWFMASGStqSPQ3c8BjwG3Lnl/LNkt4+4HgFIza05XkSLpMHV+nr6R6XUHelGRcXlLrYYuSmQFGeUSM7P6xPNK4APAy0vW2Wpmlnh+XWK/urOuRMrxwUnc1z7CJdXebbW8dHqMxUWdRpLoCTLKpQX4CzMrJh7U33L375vZPQDuvh/4BPBpM5sHpoE7XHcDkIg5NhA/mbl7y/oD/W1bNzF5foEzYzNsr69MV2kiabFqoLv7YeDqZd7fn/L8PuC+9JYmkl7H+icoKTJ2NFWvex87m+Pbnjw7qUCXyNFgWikYxwYm2NlcTelF7iG6mmSgHz87ma6yRNJGgS4Fo3tgYkPdLQBbNlVQWVrMiUEFukSPAl0KwszcAj1Dk1y6eW1zuCxVVGR0NFdzckiBLtGjQJeCcOLsJIsbHOGStLO5ihPqcpEIUqBLQejewBwuS+1srqZ3eIq5hcUN70sknRToUhCODUxQZG+e1NyInc01zC86fSPTaahMJH0U6FIQugfG2dFUTUVp8Yb3tbM5PrHXibMTq6wpkl0KdCkIx/on0tJ/DvEWOsCJs1Np2Z9IuijQJe8tLDonhybZFdt4dwtAQ1UptRUlaqFL5CjQJe+dGZ1mbsHZ0ZieQDczdsZqNNJFIkeBLnnv1FC8a2RH09pvarGSXc3VnFSXi0SMAl3yXs9wPHjb13GXopV0NFXz+rlpZuYW0rZPkY1SoEve6xmaorTY2JbGybR2Jvrje4bUSpfoUKBL3js1PElrQ9W6bju3kl2J8ew6MSpRokCXvNczNJXW7haADs26KBGkQJe85u6cGppK6wlRgJryEmKbyjmpQJcICXILugoze9rMDpnZETP74jLrmJnda2bdZnbYzK7JTLkia3Nuao7x2fm0t9AhPo2Ahi5KlARpoc8C73f3K4GrgFvN7Pol63wI2J147AO+ktYqRdYpOcJlI3cpWsnOpmpdLSqRsmqge1zyzE9p4rH0fqG3Aw8m1n0SqDezlvSWKrJ2PYl5y9Pd5QLxkS5nJ2YZm5lL+75F1iNQH7qZFZvZ88AA8Ii7P7Vkle1Ab8rrvsR7S/ezz8y6zKxrcHBwvTWLBJa8qCgTXS47Evs8paGLEhGBAt3dF9z9KqAVuM7M3r5kleXGgy1txePu97t7p7t3xmKxtVcrskY9w1NsqS1PyyyLS7UlAr1vRIEu0bCmUS7ufg54DLh1yaI+oC3ldStwekOViaTBqaGptM3hslQy0E8NK9AlGoKMcomZWX3ieSXwAeDlJas9BNyVGO1yPTDq7mfSXq3IGvUMT9Kegf5zgLrKUuoqS+kd1o0uJBpKAqzTAvyFmRUT/wXwLXf/vpndA+Du+4EDwG1ANzAF3J2hekUCm5lboH9s9kJfdya0NVaqhS6RsWqgu/th4Opl3t+f8tyBz6S3NJGNSQZtplroED/Z+vIb4xnbv8ha6EpRyVs9Q5kbg57U1lBF3/A0i4tvGQMgknUKdMlbF8agZ7TLpYrzC4v0j89k7DNEglKgS946NTzFpooS6qtKM/YZyZEuOjEqUaBAl7zVk5iUyyx90+Yu1a6hixIhCnTJW6eGMzcGPWlbfQVm0KtAlwhQoEteWlh0+kamMjrCBaC8pJiW2goFukSCAl3y0htjM8wteEbmcFmqtbGKXl3+LxGgQJe8lGwxtzak7z6iK2lvrFIfukSCAl3yUt9IfNRJW0PmW+htDVX0j80yM7eQ8c8SuRgFuuSl3uEpzGBbfRZa6E3xz0j+EhEJiwJd8lLvyBQttRWUlWT+EE/+FaB+dAmbAl3yUt/wNK1ZOCEKb45F10gXCZsCXfJS78hUVk6IAsQ2lVNeUqRAl9Ap0CXvzM4v8MbYTFZOiAKYGW0a6SIRoECXvHP63Azub86zkg1tDZWaz0VCp0CXvJO8x2dblrpcIN6P3js8RfzWACLhCHILujYze9TMjprZETP73DLr3Gxmo2b2fOLx+cyUK7K6ZEs5qy30xirGZ+cZnZ7L2meKLBXkFnTzwG+5+0Ez2wQ8a2aPuPtLS9Z73N0/kv4SRdamd2SK0mJjS21F1j4z9YbR9VVlWftckVSrttDd/Yy7H0w8HweOAtszXZjIevUOT7GtvpLiosxNm7vUhbHo6keXEK2pD93MOojfX/SpZRbfYGaHzOxhM9u7wvb7zKzLzLoGBwfXXKxIEL0j01kb4ZLU2pi8WlQjXSQ8gQPdzGqAbwO/7u5jSxYfBHa4+5XAl4HvLrcPd7/f3TvdvTMWi623ZpGL6hueoq0xeydEAWorSqmrLNXVohKqQIFuZqXEw/zr7v6dpcvdfczdJxLPDwClZtac1kpFApg6P8/Q5Hlas9xCB2hr1NBFCVeQUS4G/Blw1N3/cIV1tibWw8yuS+x3KJ2FigSRnCArW1eJpmpr0LzoEq4go1xuAj4FvGBmzyfe+12gHcDd9wOfAD5tZvPANHCHa0CuhCB5+X02hywmtTVW8fcvD7C46BRl8YSsSNKqge7uPwUuenS6+33AfekqSmS9LgR6CF0urQ2VnJ9fZHBiNqtDJkWSdKWo5JXekWkqSotorsn+WPDkLxGNdJGwKNAlr/QOT9HaUEXilE5WJUfW6MSohEWBLnmlb2Q6q3O4pGpt0LzoEi4FuuSV3pGpUE6IAlSUFhPbVK6RLhIaBbrkjdGpOcZn5kM5IZrUqml0JUQKdMkbyZZxtq8STdXWUEXfObXQJRwKdMkbp0Icg57U1ljJ6XMzzC8shlaDFC4FuuSNMC8qSmprqGJh0TkzOhNaDVK4FOiSN+JzkZdSW1EaWg3JXyY6MSphUKBL3ghj2tylLlxcpBOjEgIFuuSN3uEp2kPsbgFoqa+gyNRCl3Ao0CUvLC46r49MX7jRRFhKi4toqau8MOujSDYp0CUv9I/PcH5hMfQLine truncated
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(0.01, 2, num=100)\n",
"\n",
"base = stats.norm(1, 0.1)\n",
"\n",
"plt.plot(x, base.pdf(x))\n",
"\n",
"np.trapz(base.pdf(x), x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Transformation\n",
"\n",
"As we apply $f(x) = y = x^2$ we see that the result of the transformation doesn't suffice the requirements of a probability distribution, as $\\int P(y)dy \\neq 1$"
]
},
{
"cell_type": "code",
"execution_count": 222,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.641895835477563"
]
},
"execution_count": 222,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdp0lEQVR4nO3de3Bc53nf8e+zuN8JEguAJEhRpEk5lhzLMqRQdmK7VuzRKI7t6diJNJHjxJ7h2G1TO2nr+jKtp5lOx9Nm1NTppCnHVnyJIsfjWxzHTqzKF7WObpAsWZRoiSJNSQRJYAEQVwKL29M/dg8JwQCx2D2755zF7zNDE9w95HnmGPrx5bvP+77m7oiISPKkoi5ARESKowAXEUkoBbiISEIpwEVEEkoBLiKSULWVvFlXV5fv27evkrcUEUm8xx57bMTd06tfr2iA79u3j4GBgUreUkQk8czshbVe1xSKiEhCKcBFRBJKAS4iklAKcBGRhNowwM3sbjMbNrNjq17/AzN71syeNrP/Wr4SRURkLYWMwD8P3LryBTP7Z8A7gV9292uBPwm/NBERuZINA9zdHwDGVr38IeDT7p7NXzNchtpEROQKip0DPwT8mpk9bGY/MrMb17vQzI6Y2YCZDWQymSJvJ4WaW1jiSw+9wOz8UtSliEiZFRvgtUAncBj4d8BXzMzWutDdj7p7v7v3p9O/sJBIQvYPx87zH755jE984ym017tIdSs2wM8AX/ecR4BloCu8sqRY/3RyBIBv/GSQLz205uItEakSxQb4N4G3AJjZIaAeGAmrKCmOu/Pj50d526t6uOWV3fzx3z3DYy+s/vhCRKpFIW2E9wIPAteY2Rkz+wBwN7A/31r4ZeB9rn+vR+6lsVkGx2f51YNd3PXb17O7s4l/cc/jjM3MR12aiJRBIV0od7j7Tnevc/c+d/+cu8+7+53ufp273+Du369EsXJlwfTJ6w900dFUx12/dT1Dk1nuPz4UcWUiUg5aiVlFfnxylO62Bg6kWwC4fs826mtTnBiejrgyESkHBXiVcHcePDnC6w/sIGgIqkkZB9KtPHt+KuLqRKQcFOBV4sTwNCPT87z+wMubgQ71tHJiSAEuUo0U4FXix8/n579fseNlrx/qaePsxBxTcwtRlCUiZaQArxL/dHKUvdub6etsftnrB7tbATQPLlKFFOBVYGnZeejUKK8/sOMX3rumtw1A0ygiVUgBXgWeOTvJ1NwiN68R4Hs6m2msS/HseY3ARaqNArwKPJ/Jja6v3dXxC++lUsYruls5MawRuEi1UYBXgbPjcwDs3ta05vuHutt4TlMoIlVHAV4FzlyYZUdLPU31NWu+f7CnjaHJLBOz6kQRqSYK8CowOD7L7s61R98A1/TmO1E0ChepKgrwKjB44eK60ycAB7tznSjPDemDTJFqogBPOHfn7Pgcu64Q4Lu3NdFcX6N5cJEqowBPuAsXF5hdWLriCDyVMg52tyrARaqMAjzhBi/MAlxxDhxyH2RqCkWkuhRyoMPdZjacP7xh9Xv/1szczHScWkQGxy8C67cQBg71tDIyneWCDncQqRqFjMA/D9y6+kUz2wO8FXgx5JpkE84EI/ANAzz4IFPTKCLVopATeR4A1jpY8b8DHwV0lFqEzo7P0Vxfw7bmuitedynAtamVSNUoag7czN4BDLr7kyHXI5s0OJ5rIQwOcVhPb3sjdTV2ac5cRJKvdrO/wcyagU8Cbyvw+iPAEYC9e/du9naygY0W8QRSKaOnvZHzEwpwkWpRzAj8AHA18KSZnQb6gMfNrHeti939qLv3u3t/Op0uvlJZ0+CF2Sv2gK+0s6ORcxNzZa5IRCpl0yNwd38K6A5+nQ/xfncfCbEuKcDF+UUuXFzY8APMQG9HEz89M17mqkSkUgppI7wXeBC4xszOmNkHyl+WFOLseG46pK+AKRS4PAJ31+fOItVgwxG4u9+xwfv7QqtGNqXQFsLAzo5G5heXuXBxge0t9eUsTUQqQCsxE2wwPwLfzBw4XB65i0iyKcATbPDCLLX57pJC9Hbkgv68PsgUqQoK8AQ7Oz5Lb0cjNakr94AHghH4uUkFuEg1UIAn2OD4bMHz3wBdrQ3UpEy94CJVQgGeYIMXNhfgNSmjp61BveAiVUIBnlALS8ucn5wraBXmSr0djZoDF6kSCvCEGpqcY9kLbyEM7NzWpAAXqRIK8IQKpkEKbSEM7GzXYh6RaqEAT6jMVBaAdFvDpn5fb0cjswtLTM4ulqMsEakgBXhCjUznAryrdXMBvjPfC35WnSgiiacAT6iRqSwpY9NL4nvzveCaBxdJPgV4QmWm59neUl/wIp7ApcU8CnCRxFOAJ9TIdHbT0yeQmzNPGVrMI1IFFOAJVWyA19WkSGsxj0hVUIAnVC7Ai9sSdmdHE+e1H4pI4inAE2pkar6oETjoaDWRalHIiTx3m9mwmR1b8dp/M7OfmdlPzewbZratvGXKSjPZRWYXlujaZA94QMvpRapDISPwzwO3rnrtPuA6d/9l4Dng4yHXJVdQbA94YGdHI9PZRSbnFsIsS0QqbMMAd/cHgLFVr33P3YOlfA+RO5leKuRygBc3B66DHUSqQxhz4O8Hvrvem2Z2xMwGzGwgk8mEcDvJTM0DpY3AQb3gIklXUoCb2SeBReCe9a5x96Pu3u/u/el0upTbSV4wAt/sPiiB3vZgNaZ6wUWSbMNT6ddjZu8D3g7c4trarqKCAC/2ZPme9kbMNAIXSbqiAtzMbgX+PfAmd78YbkmykZHpLJ3NddTVFPcPqPraFF2tDZwbV4CLJFkhbYT3Ag8C15jZGTP7APA/gTbgPjN7wsz+osx1ygql9IAHetobGJ5SgIsk2YYjcHe/Y42XP1eGWqRAxS6jXynd2sBwfk9xEUkmrcRMoJHpbNGLeALdbY0KcJGEU4An0Mj0fNE94IHu9gZGp7MsLevzZ5GkUoAnzNzCEtPZxZKnULrbGlh2GJ3RKFwkqRTgCXPpLMxS58DzUzDDkwpwkaRSgCfMpWX0baVNoaTbcot5MpoHF0ksBXjCjEyXtow+0J0fgSvARZJLAZ4wpe5EGLg0haJecJHEUoAnzEh+xLyjxC6Uxroa2htr1UookmAK8IQZmc7S3lhLQ21NyX9Wd3ujplBEEkwBnjAj0/MlL+IJdLdpNaZIkinAEyYTwjL6QLpN+6GIJJkCPGFGprMl94AHutsaGJ7Mot2ARZJJAZ4wI1PZkpfRB7rbGskuLjOVXdz4YhGJHQV4gmQXl5icK30ZfaC7XasxRZJMAZ4gYzO5RTw7wpoDb1UvuEiSKcATZPTSKsyQplDatRpTJMkKOZHnbjMbNrNjK17bbmb3mdmJ/M+d5S1TAEYvjcDDCXDthyKSbIWMwD8P3LrqtY8B97v7QeD+/K+lzEYvHWYczhRKbkFQSr3gIgm1YYC7+wPA2KqX3wl8If/1F4B3hVyXrGEs5BG4meV6wSc1By6SRMXOgfe4+zmA/M/d611oZkfMbMDMBjKZTJG3E8itwqyrMdoaNjzKtGBajSmSXGX/ENPdj7p7v7v3p9Ppct+uqo3NZNnR0oCZhfZndrdpPxSRpCo2wIfMbCdA/ufh8EqS9YzNzLO9JZzpk0B3u0bgIklVbIB/C3hf/uv3AX8bTjlyJSPT86HNfwfSrQ1MzC4wt7AU6p8rIuVXSBvhvcCDwDVmdsbMPgB8GnirmZ0A3pr/tZTZ2Mw8O8owAofLB0WISHJs+GmYu9+xzlu3hFyLbGB0OhtaC2GgO98LPjyVpa+zOdQ/W0TKSysxE2JuYYmZ+aXwp1B0Or1IYinAE+JSD3jYUyiXDjdWL7hI0ijAEyLYByXsLpQdrQ2kTMvpRZJIAZ4QozPBYcbhzoHXpIwdrWolFEkiBXhCBCPwsKdQINdKqAAXSR4FeEKEvQ/KSrnFPJoDF0kaBXhCjMxkqa9J0RriPiiB4GxMEUkWBXhCjOVXYYa5D0qgp72RkeksS8s63FgkSRTgCVGOfVAC3e2NLPvl/cZFJBkU4AkxUsYA78n3gg9pGkUkURTgCTE2kw3tNPrVetpzy+mHdLCDSKIowBNidLqMI/AgwNWJIpIoCvAEmFtY4mIZ9kEJdLXWY6YpFJGkUYAnwGiZ9kEJ1Nak6GrV2ZgiSaMAT4CwT6NfS097g+bARRJGAZ4Ao2VchRnoaWvUFIpIwpQU4Gb2h2b2tJkdM7N7zawxrMLksnLugxLobm/UcnqRhCk6wM1sN/CvgX53vw6oAW4PqzC5bKxMOxGu1NPewMj0PAtLy2W7h4iEq9QplFqgycxqgWbgbOklyWqjM/PU16Zoqa8p2z2CVkLtSiiSHEUHuLsPAn8CvAicAybc/XurrzOzI2Y2YGYDmUym+Eq3sNHpebpayrMPSqCnPViNqWkUkaQoZQqlE3gncDWwC2gxsztXX+fuR92939370+l08ZVuYWMz82wv4weYsOJwYwW4SGKUMoXy68DP3T3j7gvA14HXh1OWrFSO0+hX6+0IltNrCkUkKUoJ8BeBw2bWbLl/298CHA+nLFlpdCY3hVJO25vrqU2ZplBEEqSUOfCHga8CjwNP5f+soyHVJSuUcx+UQCpldLc1aAQukiAlHe/i7p8CPhVSLbKGmewiswtLdLWVdwoF1AsukjRaiRlzI/ll9Oky9oAHtJxeJFkU4DGXyfdlV2IE3tOu5fQiSaIAj7lgBN5V5jZCyAX4xOwCcwtLZb+XiJROAR5zmfw+KJWYQunOj/J1Qr1IMijAYy4zlcWMsnehgE7mEUkaBXjMjUxncz3aNeX/v0pnY4okiwI85kamyneY8WqX90PRFIpIEijAYy4znaWrrfzTJwAdTXXU16a0H4pIQijAY25kOluRDzABzEy94CIJogCPuZGp+YpNoYCOVhNJEgV4jFVyGX2gp71RXSgiCaEAj7FLqzArOALvbm9QH7hIQijLine truncated
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def f(x):\n",
" return x**2\n",
"\n",
"transformed = f(base.pdf(x))\n",
"plt.plot(f(x), transformed)\n",
"\n",
"np.trapz(transformed, f(x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Change of variables\n",
"\n",
"We want to transform $P(x)$ to $P(y)$ with $f: \\mathbb{R}^n \\mapsto \\mathbb{R}^n $ and $f(x) = y$\n",
"The transformation blows up or shrinks the area/ volume (in N-D) of the distribution, so we need to determine a transformation so that:\n",
"\n",
"\n",
"$$ \n",
"\\begin{eqnarray}\n",
"\\int P(x)dx &=&\\int P(y)dy = 1 \\\\\n",
"P(x)dx &=&P(y)dy \\\\\n",
"P(y) &=&P(x)\\frac{dx}{dy}\n",
"\\end{eqnarray} $$\n",
"\n",
"So to keep an integral of 1 for $P(y)$ (the definition of a probability distribution), we need to multiply $P(x)$ with its derivative w.r.t. $y$.\n",
"Therefore we need to express $x$ in $y$. This can be done if we only have invertible transformations.\n",
"\n",
"$$\\begin{eqnarray}\n",
"f(x) &=& y \\\\\n",
"f^{-1}(y) &=& x\n",
"\\end{eqnarray} $$\n",
"\n",
"$$\\begin{eqnarray}\n",
" P(y) &=& P(x)\\frac{\\text{d}f^{-1}(y)}{\\text{d}y} \\\\\n",
" &=& P(x)f'^{-1}(y) \\\\\n",
" &=& P(f^{-1}(y))f'^{-1}(y)\n",
"\\end{eqnarray} $$\n"
]
},
{
"cell_type": "code",
"execution_count": 223,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9987379589284238"
]
},
"execution_count": 223,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xc9ZXw/8+ZUS+WbEmWi2zJvYArxsYQiiEQTIkhIU9MEpJNsnFgCSRP2gaSJU/Kbspv8zyEkOA4CbBsKMkGAoYYQi8BbGMbV1xwkWzZMlazqlXn/P64M7Ysq4ykO/28X69BM3PbQb46853vPff7FVXFGGNMYvFEOgBjjDHhZ8nfGGMSkCV/Y4xJQJb8jTEmAVnyN8aYBJQUqQPn5+drSUlJpA5vjDExaePGjVWqWjDU/UQs+ZeUlLBhw4ZIHd4YY2KSiJS5sR/r9jHGmARkyd8YYxKQJX9jjElAlvyNMSYBBZ38RcQrIu+KyDM9LBMRuUdE9orIVhGZ726Yxhhj3DSQlv9XgZ29LFsKTPE/VgD3DTEuY4wxIRRU8heRIuBq4Pe9rLIMeEgda4FcERntUozGGGNcFmzL/27g24Cvl+VjgUNdXpf73zuNiKwQkQ0isqGysnJAgRoTDqu3HGH74bpIh2FMyPWb/EXkGuCYqm7sa7Ue3jtjogBVXaWqC1R1QUHBkG9QM8ZVb+6t4vZH3+Xae//BHU9spbqxNdIhGRMywbT8LwA+KiKlwGPApSLyx27rlAPjurwuAo64EqExYdDpU370zHuMzU3nCxdM4H82lLPkP1/liU3lkQ7NmJDoN/mr6h2qWqSqJcBy4GVV/Uy31VYDn/VX/ZwH1KlqhfvhGhMaf3rnELuONnDnVTP4t2tm8tzXLmTSyCzueGIbja0dkQ7PGNcNus5fRG4WkZv9L9cA+4G9wO+Af3EhNmPCor6lnV88v5tzS4Zz1axRAEwemc13r5pBa4ePF947GuEIjXHfgAZ2U9VXgVf9z1d2eV+BW90MzJhw+fUre6lpbuPBaxYicury1fzxwxmbm85Tm49w/byiCEZojPvsDl+T0A7VNPPAP0r5+PwiZhXlnLbM4xGunTOGN96vsou/Ju5Y8jcJbfWWI7R1+vjGFVN7XL5s7hg6fcqa7db1Y+KLJX+T0N7eV830UdmMzknvcfn0UdlMGZnF05uteM3EF0v+JmG1dnTyTmkN50/K73UdEWHZ3DGsL63hyPETYYzOmNCy5G8S1rsHj9Pa4eP8SXl9rnftnDEAPL3FWv8mfljyNwnrrX3VeAQWThzR53rFeZnMHZfLakv+Jo5Y8jcJa+2+amaNzWFYWnK/6350zhh2HKlnf2VjGCIzJvQs+ZuE1NzWwbuHalncR39/V5dMc8aieqe0JpRhGRM2lvxNQtpQWkt7p/bb3x9QkpdJdloSW8ptxE8THyz5m4T01r5qkr3CgpLhQa3v8QhzinLZWn48xJEZEx6W/E1Cent/NXPH5ZKREvwIJ7OLcthV0UBLe2cIIzMmPCz5m4RT39LOtvLjQff3B8wuyqXDp7xXUR+iyIwJH0v+JuGs31+DTwm6vz9g7rhcALYesq4fE/ss+ZuE89a+alKTPMwbnzug7UblpDEyO5WtdtHXxAFL/ibhbDpYy9xxuaQmeQe87eyiXDbbRV8TByz5m4Ti8ym7jzYwc8ywQW0/d1wO+yubqG9pdzkyY8IrmAnc00RkvYhsEZEdIvKDHta5RETqRGSz/3FXaMI1ZmgO1jRzor2TGaMGl/xnFzldRdut68fEuGDq3FqBS1W1UUSSgX+IyLOqurbbem+o6jXuh2iMe3YddSp1po/OHtT2s/0Tvmwpr+P8yQOrFjImmvSb/P1TNAYGNEn2PzSUQRkTKjsrGvAITBk5uOSfm5FCSV4GW6zix8S4oPr8RcQrIpuBY8ALqrquh9UW+7uGnhWRs3rZzwoR2SAiGyorK4cQtjGDs+toPSX5maSnDPxib8Bsu9PXxIGgkr+qdqrqXKAIWCgiZ3dbZRNQrKpzgF8BT/ayn1WqukBVFxQUFAwlbmMGZffRBqaPGlyrP2B2UQ5H6lo41tDiUlTGhN+Aqn1U9TjwKnBlt/frVbXR/3wNkCwi1iFqokpTawdlNc1MH+TF3oA5J2/2sou+JnYFU+1TICK5/ufpwIeBXd3WGSUi4n++0L/favfDNWbw9nzQgCpDbvmfNWYYXo9Y14+JacFU+4wG/ktEvDhJ/c+q+oyI3AygqiuBG4BbRKQDOAEs918oNiZq7DraAMCM0UNr+WekJDExP5Od/v0ZE4uCqfbZCszr4f2VXZ7fC9zrbmjGuGtXRT1ZqUmMzU0f8r4mj8xityV/E8PsDl+TMHYebWDaqGw8HhnyviaPzKKsppm2Dp8LkRkTfpb8TUJQVXZV1A+5vz9g8sgsOn1KaXWTK/szJtws+ZuEUFHXQn1Lh2vJf1JBFgDvf2ATupvYZMnfJIRA//z0IV7sDZhUkIUI7D1myd/EJkv+JiHs9I/pM82lln96ipei4ensrbTkb2KTJX+TEHZVNDA2N51hacmu7XNyQZa1/E3MsuRvEsKuo/XMGORInr2ZPDKLfZWNdPrslhYTeyz5m7jX2tHJvsqmIQ/r0N3kkVm0dfgor212db/GhIMlfxP3Squa6fQpUwqzXN3vZP+w0Nb1Y2KRJX8T9w5UOcl5Yr7byd/ZnyV/E4ss+Zu4t7/KuRGrJD/D1f3mpCdTkJ3K+5b8TQyy5G/i3oHKJkZmp5LtYqVPgFX8mFhlyd/Evf1VTUzIzwzJviePzGLfsUZsEFsTayz5m7h3oKqJiQWhSf5TCrNoaO3gWENrSPZvTKhY8jdx7XhzGzVNba5f7A2YbGP8mBhlyd/EtcDF3lB2+wDsPWZj+5vYEsw0jmkisl5EtojIDhH5QQ/riIjcIyJ7RWSriMwPTbjGDMyBSn/yD1G3T0F2KtlpSTbGj4k5wUzj2ApcqqqNIpIM/ENEnlXVtV3WWQpM8T8WAff5fxoTUQeqmvB6hHHD3S3zDBARJo+0ih8Te/pt+asjcGYn+x/dSxuWAQ/5110L5IrIaHdDNWbg9lc1Mn5EBilJoevhnGLJ38SgoP4iRMQrIpuBY8ALqrqu2ypjgUNdXpf73+u+nxUiskFENlRWVg42ZmOCtr8ydGWeARMLsqhqbKO+pT2kxzHGTUElf1XtVNW5QBGwUETO7rZKT5OinlH4rKqrVHWBqi4oKCgYeLTGDIDPP81iqJN/SZ7TpVRWZQO8mdgxoO/CqnoceBW4stuicmBcl9dFwJEhRWbMEFXUt9DS7gtZjX9Aif/D5YDN52tiSDDVPgUikut/ng58GNjVbbXVwGf9VT/nAXWqWuF6tMYMwMlKnxC3/ItHOPsvq7Lkb2JHMNU+o4H/EhEvzofFn1X1GRG5GUBVVwJrgKuAvUAz8PkQxWtM0EI1mmd36SleRg1Ls5a/iSn9Jn9V3QrM6+H9lV2eK3Cru6EZMzT7KpvISPFSOCw15Mcqyc+g1Fr+JobYHb4mbh3wD+gm0lM9grsm5GdSVm0XfE3ssORv4taBEI7m2V1xXibVTVbuaWKHJX8Tl1o7OimvbWZimJJ/SV7goq+1/k1ssORv4tLB6mZ8GroxfboLzBJmF31NrLDkb+JSYDTPUFf6BATKPe2ir4kVlvxNXDpwct7e8LT801O8jM5Jo9Ra/iZGWPI3camsuokRmSnkpLs/b29vivOs3NPEDkv+Ji6VVjVTnBeaYZx7MyE/k1Ir9zQxwpK/iUtl1U0nK3DCpSQvk5qmNupOWLmniX6W/E3caWnv5EhdS9hb/sWBck/r9zcxwJK/iTuHapyul3C3/AM3lFnXj4kFlvxN3AkMsxD+lr9zPLvoa2KBJX8TdwLlluFu+acl+8s9LfmbGGDJ38SdsupmhqUlkZsRvjLPgJK8TKv1NzHBkr+JO6XVTZSEaTTP7kqs3NPECEv+Ju6UVTefrLwJt5K8DCv3NDEhmGkcx4nIKyKyU0R2iMhXe1jnEhGpE5HN/sddoQnXmL61dfgor20+Oal6uAWGk7ByTxPtgpnGsQP4hqpuEpFsYKOIvKCq73Vb7w1Vvcb9EI0J3uHjJ/ApjB8RoeTv/8ZxoKqJ2UW5EYnBmGD02/JX1QpV3eR/3gDsBMaGOjBjBuNkpU+YBnTrLlDuedD6/U2UG1Cfv4iU4Mznu66HxYtFZIuIPCsiZ/Wy/QoR2SAiGyorKwccrDH9ORihGv+Ak+WelvxNlAs6+YtIFvA48DVVre+2eBNQrKpzgF8BT/a0D1VdpaoLVHVBQUHBYGM2plel1c6k7QVZoZ+0vTfFeRnW52+iXlDJX0SScRL/w6r6RPflqlqvqo3+52uAZBHJdzVSY4IQqPSJRJlngFPrby1/E92CqfYR4A/ATlX9v72sM8q/HiKy0L/fajcDNSYYpdVNEav0CRifl0FVYyuNrR0RjcOYvgRT7XMBcBOwTUQ2+9+7ExgPoKorgRuAW0SkAzgBLFdVDUG8xvSq06ccqmnmipmjIhpHSZfRPc8akxPRWIzpTb/JX1X/AfT5HVpV7wXudSsoYwbjyPETtHdqxFv+gYvNZdXNlvxN1LI7fE3cODWaZ2TKPAMCx7cxfkw0s+Rv4sapGv/ItvyzUpPIz0qlrMou+proZcnfxI2DNc2kJnkozE6LdCiU5GVYy99ENUv+Jm6UVjVRnJeBxxO5Ms+A4rzMk91QxkQjS/4mbpRWNzF+RGT7+wNK8jI4Wt/CibbOSIdiTI8s+Zu44PMpZdXNTIhwf39AsX9soYM11vo30cmSv4kLR+tbaO3wRbzSJyBQbmr9/iZaWfI3cSGQZCdEaDTP7or93U82uqeJVpb8TVwo9ZdVRmoo5+5yMpIZnpFsLX8TtSz5m7hQWt1ESpKH0cMiX+YZYBU/JppZ8jdxobSqieIR0VHmGWCLine truncated
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def f_i(y):\n",
" return y**0.5\n",
"\n",
"def f_i_prime(y):\n",
" return 0.5*y**-0.5\n",
"\n",
"assert np.allclose(f_i(f(x)), x)\n",
"\n",
"y = f(x)\n",
"px = base.pdf(x)\n",
"\n",
"transformed = px * f_i_prime(f_i(y))\n",
"\n",
"plt.plot(x, px)\n",
"plt.plot(y, transformed)\n",
"plt.xlim(0, 2)\n",
"\n",
"np.trapz(transformed, y)"
]
},
{
"cell_type": "code",
"execution_count": 224,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9987894616589735"
]
},
"execution_count": 224,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhc9X3v8fdX+2ptHtvyvuANYxuM6rAFDCRAWNM8SQlt0zRN6rZZmrZpe9t7+7RJ2/vc7k3am/ZemiYthJBLEmgp2Qhg4xC2COMFr9gGgyzZ2mzJWqxtfvePM2PJsqQ5I82Z5fB5PY8z0sw5Z74Rxx/9/J3f+R1zziEiIuGQl+kCREQkdRTqIiIholAXEQkRhbqISIgo1EVEQqQgiIPOnj3bLV26NIhDi4iE0iuvvNLunIvM9DiBhPrSpUtpbGwM4tAiIqFkZsdTcRy1X0REQkShLiISIgp1EZEQUaiLiISIQl1EJER8hbqZ/baZ7TOz18zsYTMrCbowERFJXsJQN7MFwG8CDc65y4B84MNBFyYiIsnz234pAErNrAAoA5qDK0kkQJ1vwOEfZroKkcAkDHXn3Angb4C3gBagyzn35PjtzGyrmTWaWWNbW1vqKxWZqe4W+Lc74Bs/Bz06RyWc/LRfaoB7gGXAfKDczH5x/HbOufudcw3OuYZIZMZXuoqk1vAAPHwv9Jzyvj+2LbP1iATET/vlPcAbzrk259wQ8ChwTbBliaTY8eehZTfc82UorYUjT2e6IpFA+An1t4CrzKzMzAy4GTgQbFkiKda633tccTOsuBGOPgPRaGZrEgmAn576S8C3gZ3A3tg+9wdcl0hqndoP5RGoiHjB3tsKrfsyXZVIyvma/eKc+xPn3Brn3GXOuY845waCLkwkpVr3w5xLva9X3OQ9qgUjIaQrSiX8olFoOwhz13nfz6qHukug6aeZrUskAAp1Cb/Tb8BQH8xZO/pczTLoejtzNYkERKEu4Rf/kHTOutHnqhfDmbcyU49IgBTqEn6n4qG+ZvS56kXQfxrOdWemJpGAKNQl/Fr3Q81SKCoffa56sfeoFoyEjEJdwq/1wIWtF4DqJd6jWjASMgp1CTfn4MxxqF124fPxkfoZjdQlXBTqEm79p2H4HMyaf+Hz5REoKPECXyREFOoSbmdbvMfK+gufN9MMGAklhbqEW3cs1MeP1EGhLqGkUJdwOxu7n8v4kTpA1SKFuoSOQl3CLT5Sr5x38WvVi6G/EwZ60luTSIAU6hJuZ5uhrA4Kii9+TXPVJYQU6hJu3S1QOUE/HTStUUJJoS7hdrbFW5VxIhVzvMde3a9UwkOhLuF2tmXiD0nBm6sOCnUJFT83nl5tZrvG/Ok2s99KR3EiMzI86AX2RNMZwVsLprBMoS6hUpBoA+fcIeByADPLB04AjwVcl8jM9Zz0HicbqQOUz4be9vTUI5IGybZfbgaOOud0bbVkv6kuPIorj2ikLqGSbKh/GHh4ohfMbKuZNZpZY1ub/pJIFpjqwqM4hbqEjO9QN7Mi4G7gWxO97py73znX4JxriEQiqapPZPp8jdTVfpFwSWak/j5gp3PuVFDFiKTU2WbIL4bSmsm3KZvtjdSdS19dIgFKJtTvY5LWi0hW6mmFyrneioyTKY9AdAjOdaWvLpEA+Qp1MysD3gs8Gmw5IinU2z46F30y5+eqqwUj4eAr1J1zfc65OuechjOSO3rbvPbKVMpnj24rEgK6olTCq69jNLQno6tKJWQU6hJOznktlbK6qbdTqEvIKNQlnAZ7YGQg8Ug9HvrqqUtIKNQlnOIhnainXlAEJdUaqUtoKNQlnPo6vMdEI3XQVaUSKgp1CSe/I3WIhbraLxIOCnUJp75YSJcn+KAUYksFaKQu4aBQl3BKaqSuUJfwUKhLOPW1Q0GJdyOMREpr4dwZiEaDr0skYAp1CafeDm+UPtW6L3FlteCiMNAdfF0iAVOoSzj1tfvrp8PoKo79p4OrRyRNFOoSTr3t/vrpoFCXUFGoSzj1tvubow5jQr0zuHpE0kShLuHUl8xIvdZ77D8TXD0iaaJQl/AZ7IOhPvXU5R1JoS7h05fEHHWA0urYfmq/SO7ze+ejajP7tpkdNLMDZnZ10IWJTFv8wiO/PfX8QiiepZG6hEKBz+2+BPzAOfdBMysCygKsSWRm4ot5+R2pgzdaV6hLCCQMdTObBVwP/DKAc24QGAy2LJEZSGaFxrjSGoW6hIKf9styoA34mpm9amZfMbOLrr02s61m1mhmjW1tWkdDMijeG49/AOpHaY2mNEoo+An1AmAT8M/OuSuAXuAPxm/knLvfOdfgnGuIRBLcwV0kSP2dYHnezS/8Kq3VSF1CwU+oNwFNzrmXYt9/Gy/kRbJTX6cX6HlJTO5S+0VCIuFZ75w7CbxtZqtjT90M7A+0KpGZ6O/0FulKRjzUtVKj5Di/s18+AzwUm/lyDPhYcCWJzFBf5+hVon6NXamxNIm2jUiW8RXqzrldQEPAtYikRn8nzFqQ3D5jrypVqEsO0xWlEj59p5MfqWupAAkJhbqEz3R76vF9RXKYQl3CZeict5hXMnPUQSs1Smgo1CVc+qdx4dHY7dV+kRynUJdwiYdy0u0XrdQo4aBQl3A5v0RAkqGeXwhFlRqpS85TqEu4xNsvyY7UwRutn1NPXXKbQl3CZbojdfCWFtAHpZLjFOoSLhqpyzucQl3Cpa8TCkqhsDT5fUuqNFKXnKdQl3DpPz29UTpopC6hoFCXcJnOYl5x6qlLCCjUJVz6O6EsyQuP4kqrYbgfhgdSW5NIGinUJVxmOlIHjdYlpynUJVyms5hXXHypAPXVJYf5Wk/dzN4EzgIjwLBzTmurS/aJRmProU93pF7lPWqkLjnM752PAG50zrUHVonITA10eXcvmu5IPd5+0UhdcpjaLxIe8XVbpjtSL1VPXXKf31B3wJNm9oqZbZ1oAzPbamaNZtbY1taWugpF/OqLh/o0Z79opC4h4DfUr3XObQLeB3zKzK4fv4Fz7n7nXINzriESiaS0SBFfZrJEAGikLqHgK9Sdc82xx1bgMWBzkEWJTMtMFvMCb/ndwnKN1CWnJQx1Mys3s8r418AtwGtBFyaStJmO1MEbrWukLjnMz+yXucBjZhbf/hvOuR8EWpXIdPSfBmx0auJ0lGj9F8ltCUPdOXcM2JiGWkRmpq/TC/S8/OkfQyN1yXGa0ijhMZOrSeM0Upccp1CX8JjJui9xGqlLjlOoS3hopC6iUJcQ6T89/QuP4kqrYagPhgdTU5NIminUJTz6ZrCYV5yuKpUcp1CXcBgehMGzM2+/6KpSyXEKdQmH/hmu+xKnkbrkOIW6hEMqriaF0V8KGqlLjlKoSzikaqReqpG65DaFuoTDTBfzitN9SiXHKdQlHFLWfomH+umZHUckQxTqEg6pGqlr+V3JcQp1CYf+TsgvgqLymR+rtEbtF8lZCnUJh/jVpN4S0TNTqqUCJHcp1CUcUrGYV1yJFvWS3KVQl3DoPz3zD0njNFKXHOY71M0s38xeNbMngixIZFr6Omc+Rz2utFqzXyRnJTNS/yxwIKhCRGYkFcvuxqn9IjnMV6ib2ULgDuArwZYjMg3OQV8HlM2edJOhkSjtPQP+jldaDcP9MOxze5Es4ufG0wBfBH4fqJxsAzPbCmwFWLx48cwrE/FroBuiw1BWN+HLJ87087GvvczhUz2smVfJltVzuGXdXDYtnqRdM/aq0sq5ARUtEoyEI3UzuxNodc69MtV2zrn7nXMNzrmGSCSSsgJFEurr8B4nCPXXTnTxs1/+CS1d5/jMTZdQXVbIV358jA/80/M80vj2xMeL9+b1YankID8j9WuBu83sdqAEmGVmX3fO/WKwpYn41DtxqG8/1MonH9pJdWkh3/mNa1g11/uHZve5IbY+0MgXHt/HVcvqWFxXduHxtFSA5LCEI3Xn3B865xY655YCHwaeUaBLVplgpN58pp+tD77C0rpyHvvUtecDHWBWSSF/+3OXk5dn/PYjuxgeiV54vBItvyu5S/PUJfedD/XR2S//+MwRnHPc/0tXMndWyUW7LKgu5c/ffxmvHD/NP28/euGLWn5XclhSoe6c2+6cuzOoYkSmJR7q5d7sl+MdvXyr8W3u27yYhTVlk+52z+ULuOfy+Xzx6dfZ/faYANfyu5LDNFKX3NfXEVvMqwKALz39Ovl5xqduvCThrn96z2XMrijiL75/cPTJkirvUSN1yUEKdcl9fR1eP92MI61n+Y9XT/DRa5ZO2HYZr6q0kI9es5QXjnXw+qmz3pP5BVA8Sx+USk5SqEvui4c68PdPvU5pYT6/dv1y37vf27CIovw8Hnzx+OiTuqpUcpRCXXJfXweU1XKgpZvv7mnhV65bRl1Fse/d6yqKuXNDPY/uPEHPwLD3ZGmV2i+SkxTqkvtiI/X/99O3KSrI4xPX+R+lx33k6iX0DAzz2M4m74nSWrVfJCcp1CX39XUwUlrLf+1u5r1r51JVVpj0IS5fVM36BVU88MJxnHPeVaXxW+SJ5BCFuuS2kWHoP8Nb/WV09A7y/isWTOswZsZHrl7C6609vHDMa+ecv5m1SA5RqEtuO3cGcDS2GTVlhdywavrrDt29cT7VZYU8+MLx0fZLNJp4R5EsolCX3Ba78OiFk3DnhvkUFUz/lC4pzOdDVy7kyf2n6CuoAhf1VoAUySEKdcltve0AnBoun3brZay7Ns5nJOp4rTP2V0MtGMkxCnXJbbGRevGsCJsWV8/4cOsXVLGgupQXWlzs+JoBI7lFoS45rbvzJABXXbYKM5vx8cyMW9fN4yfNsV66RuqSYxTqktMOHfOuAn3vlWtTdsz3rZ9H20i5942mNUqOUahLTmtubqLPSllWP/n9SZN15eIa8stja7PrAiTJMQp1yVmt3ecY6W1nuGjmvfSx8vKMq9ctI+qMoZ72lB5bJGgKdclLine truncated
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Same principle, other transformation\n",
"\n",
"def f(x):\n",
" return x**0.5\n",
"\n",
"def f_i(y):\n",
" return y**2\n",
"\n",
"def f_i_prime(y):\n",
" return 2*x\n",
"\n",
"assert np.allclose(f_i(f(x)), x)\n",
"\n",
"y = f(x)\n",
"px = base.pdf(x)\n",
"\n",
"transformed = px * f_i_prime(f_i(y))\n",
"\n",
"plt.plot(x, px)\n",
"plt.plot(y, transformed)\n",
"plt.xlim(0, 2)\n",
"\n",
"np.trapz(transformed, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conditions\n",
"The transformations we've used were strictly increasing. \n",
"\n",
"In the example we've used above $\\frac{dx}{dy}$ was positive, as the transformation was strictly increasing, in the strictly decreasing case, $\\frac{dx}{dy}$ is negative and $P(y)=-P(x)\\frac{dx}{dy}$. If we take the absolute value of the derivative we have one equation for both cases.\n",
"\n",
"$$ \n",
"\\begin{eqnarray}\n",
"P(y)&=&P(x) \\cdot \\left| \\frac{dx}{dy} \\right| \\\\\n",
" &=& P(f^{-1}(y)) \\cdot \\left| f'^{-1}(y)\\right|\n",
"\\end{eqnarray}\n",
"$$\n",
"\n",
"The intuition of taking the modulus is that if the (local) rate of change of $x + dx = y + dy > 0$, than we can ensure ourselves that the rate of amount of probability is increasing, and the total amount of probability is preserved. This is ensured by conditioning $\\frac{dx}{dy} \\geq 0$.\n",
"\n",
"![](img/absolute.png)\n",
"\n",
"## Multiple dimensions\n",
"In N-D $\\frac{dx}{dy}$ is expressed in the Jacobian matrix, i.e. the matrix with first order partial derivatives $\\frac{\\partial(x)}{\\partial{y}}=\\mathbf{J}$. This absolute value, is now a determinant, such that;\n",
"\n",
"$$ \n",
"\\begin{eqnarray}\n",
"P(y)&=&P(f^{-1}(y)) \\cdot \\left| \\det \\frac{\\partial{ f(y)^{-1} }}{\\partial{y}} \\right| \\\\\n",
"P(y)&=&P(f^{-1}(y)) \\cdot \\left| \\det \\mathbf{J}_{f^{-1}}(y) \\right| \\\\\n",
"\\end{eqnarray}\n",
"$$\n",
"\n",
"The geometric meaning of the determinant is a scaling factor. $\\det A$ tells us how much an n-dimensional volume is scaled by applying the transformation of $A$."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD5CAYAAAADZljUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9WYxcWXrn9zt3v7EvuXIna2FVdVX1KnWjpda0NGrIhoWxDAgYG7bhJ+t1DBgwYL/Ms18MG/CLBY9hDGSMH2R5xvDYFnrGVm+j7lZtzaquKrKKTDL3JfaIuy/HDycyKpPMJJPMrCqyef9AocjM4L03Is53zrf8v/8npJQUKFDg+YD2ZT9AgQIFvjgUBl+gwHOEwuALFHiOUBh8gQLPEQqDL1DgOUJh8AUKPEcwzupCQggdeAvYkFL+8cNeOzc3J69cuXJWty5QoMABvP322x0p5fxRvzszgwf+EfARUHvUC69cucJbb711hrcuUKDAPoQQ94773Zm49EKIC8C/A/yPZ3G9AgUKfD44qxj+vwX+CyA/o+sVKFDgc8CpDV4I8cfArpTy7Ue87s+EEG8JId7a29s77W0LFCjwBDiLE/53gH8ghLgL/K/AHwgh/uL+F0kp/1xK+S0p5bfm54/MJxQoUOBzxqkNXkr5X0opL0gprwD/PvD/Sin/o1M/WYECBc4cRR2+QIHnCGdZlkNK+TfA35zlNQs82xj4MSsdj1GQUHNNrs6VaZSsL/uxnlsUJ3yBzw0DP+bd1T5xmtMsWcRpzrurfQZ+/GU/2nOLwuALfG5Y6XiULIOSZSCEmP15peN92Y/23KIw+AKfG0ZBgmvqh37mmjqjIPmSnqhAYfAFPjfUXJMgyQ79LEgyaq75JT1RgcLgC3xuuDpXxo9T/DhFSjn789W58pf9aM8tCoMv8LmhUbL4+qUmlqHR92MsQ+Prl5pFlv5LxJmW5QoUuB/K6AsDf1pQnPAFCjxHKAy+QIHnCIXBFyjwHKEw+AIFniMUBl+gwHOEwuALFHiOUJTlCpw5ig65pxfFCV/gTFF0yD3dKAy+wJmi6JB7ulEYfIEzRdEh93SjMPgCZ4qiQ+7pRmHwBc4URYfc043C4AucKYoOuacbpy7LCSEc4MeAPb3eX0op//Fpr1vg2UXRIff04izq8BHwB1LKiRDCBH4qhPi/pZQ/P4NrFyhQ4AxxaoOXUkpgMv2rOf1Pnva6BQoUOHuc1fRYXQjxHrAL/FBK+YuzuG6BAgXOFmdi8FLKTEr5NeAC8NtCiNfvf00xTLJAgS8fZ5qll1IOUJNn/q0jflcMkyxQ4EvGWWTp54FESjkQQrjAHwL/9amf7DlH0YBS4PPAWZzwy8D/J4S4AfwdKob/P8/gus8tigaUAp8XziJLfwP4+hk8S4EpDjagALP/r3S8or5d4FQo+uGfQoyChOYB930cJmz0fXZGIUDh3hd4YhTU2qcQBxtQxmHCze0RXpSyWHMK977AqVCc8E+AzzuhdnWuzLurfQA2+j4CkAjON0uFe1/gVChO+MfEF5FQO9iAsjMKKdsG15eqVB3VYlr0lxd4UhQn/GPii0qoHWxAidN8dh8o+ssLPDkKg39M3J9QA3Xi9j+nmPqge++aOkGS4ccp15ean8v9nmYU3ITTozD4x8R+Qu2LOnH33fuVjkffj6m5JteXPp/+8ic1qC/CEPdDqZJl0CxZBEnGu6v9otf+MVHE8I+JL0PRZd/o/971hc9tgT9pbuKLIgkV4phng8LgHxNfhKLLvhH96ObuF1KCG/gxf/3BNrd2Jqz1fCZRemKD+qIMsRDHPBsULv0T4PNUdHmU63rW7vP+/QZBzFzZJs4kN7fHXF+qUrGNR+YmjspppFnOrzcHZ+rif9Gh1G8qihP+KcPDTszPw33ev1+7bBNnOY6p45g6m4PgRAZ1v0rtOEy4sTHE0vUnfsajPJxCHPNsUBj8U4aHua431gdsDAI+2hpza2dMlstTu8/79zvXcAmTjDBJsXRB14tPZFD3G+KdzgSB5Np85Ylc/OM2NaAQxzwDFC79U4bjXFch4L3VAXMVm5qjE6U5N7fHvLxYIbxPB/5J7ld1TK4v1dgcBHS8iIZrncig7q8iRGnOG+cbM5IQPF7Z8uE8h2bBLjwlCoN/ynBc3V3XBK2yhRAghMCZegErXY+vnKufyf0qtsHFVol25WTGvo+DOY2aaxKn+aHfP06s/UXzHJ43FC79U4bjqgBSKuPcd7ullEiZ05tEp4pjT1N1eJxYu1W2TlR5KCbXfL4QSnT2i8W3vvUt+dZbb33h932WsR/XZrlkcxAwiVIMTXCpXeL3Xl74wp/nYDXhoCfy9UuKAXiwktAqW9zZmxz52vs3loddt4jXTwYhxNtSym8d9bvCpX9GsO96lyyDlxerM0N480LjS3mex4m195/7JP0HXySz8HlEYfBfEE5bP39aDGH/ffzskz0Waw7nm6VDXXxHxdqPG5cXk2s+PxQG/wXgrHjgX7YhHHwfizUHL0q5uT3i+lKNqmMeG2sXpJmnB4XBfwF4FjTqTuKBHHwfddfi4+0xUZLS9WK+frGJpnFkF1/R8ff04Cxkqi8C/xRYAnLgz6WU/91pr/ubhCcpNT1JCHCabreTeCD772McJmwMAs41XIZ+zM44ZKXr8YPXFo+83xcdjhRttMfjLE74FPjPpZTvCCGqwNtCiB9KKT88g2v/RuBxXdonCQH2/02eQ8+L+Hh7xNv3+vzgtUUut8uz1xxlCCsdjzxn1jhTsdV97/dA9t/H5iCYUXBLls6ldpmLrRI9L57d6348Khw5KyMt2mgfjlPX4aWUW1LKd6Z/HgMfAedPe93fJJyUB76/WP/qnTU2+gFZLk9MT9032tWeR5pL5isOhoAffrjNwI8fysPfHATc605IspyaYzIMYn726S7/6sOdQzXz/ffR9WIsXRAmKWGSUXNMVrsqkfck3P6z7BEo2mgfjjMl3gghrqA06othkgdwEnLLwUWvoaEJuLk9Yhyq9s9HtYKOgoSeF01PXrXYa65FlisjeJghTMIUTQgc08CPM7ZHEXkOAg4Z3/77aLgmHS/C1DXON0psDIJTqeqexEhP2jJctNE+HGeWtBNCVID/DfjPpJSjI37/Z8CfAVy6dOmsbvvM4FEu7cFFX3EMkkzimILNQcC5BtzZmxBn2bHubs01+Xh7xHzFmf0sSnNaJXO22I/LI1QcAy9Wp/XeOERIgWFolB3jgQRjo2TxR68vzdzm1a6HQCLhiVV1H5XjeBw3/XHCp+cx1j+rcdEmytj/FynlXx31mmKY5MNx8GTa71yTEnbHITfW+4yjlOuLtWNP0KtzZXRNYxTESCmnFNyMVtmm5prHUla9OOXXmyNurA94516Xu10P1xIs113mKzbw4Al5tKpu7USquked1I+i0z6Om/644dPzNs7rLLL0AvgnwEdSyv/m9I/0bONJT42DJ5PqXKtypzOhM465tlDm2lzlUAfa/Sdoo2Txg9cW+T/e22Cl62FqGgt1Gz9J+eacKn/dXxpb6Uy41/GouibtskWWSQbDgFyCJtTGA0efkE+iqnvcSX1tvsKdvcmhZztYtjtJlePg565rgp4XszkMEEheXKg+8CzPQqn088BZuPS/A/zHwPtCiPemP/uvpJT/1xlc+6nHwYUmBAyDhIWq80jX8/6NYZ9vDmox65rgfMOlWbK42Cyh9tWHj52qTzeYZskiyXJMXUNM73dUacyPVdxdL9nUXIvOOGSSZNzpeHztYhMpJb9a79ObRHztUnMWxx/E49TYjzOynhc/tGz3KDf9/o1kbxxxa3vE6+cbzFftI7+HUZBgaIKb26NZZWK57pyq1fhZwFkMk/wpzNbVc4X7F9qN9QHjKKVdtmeuJzx4ahx10t3Zm3BtvkLPiw8t+pWON1vs+2OnBBxKkO0v5JWOx3zV4XK7MruXH6eH4u+Dz/GzT/ZYrqtTvGwZlNsVLjZL3Ol4mIbg7Xt9WmWLb15uYejakZvX49TYH3ZSPyzHcdymslSv8O5qn3dW+9iGxrU5JbrR92Pqrknfj3Etpd7T9WJ2RxF/9PoSjZJqM35/Y0Ddtag5JlGa8f7GgNdO0Wr8LKBojz0F7o8t01xSd0w2B8HsNUfFs8fFpD1PtZfWXJVoW+l4tMrWLA49auzUwVj2cTPU81VnVgXYxzhMuNRS/PhvX2vz1YtNaq710Lj5pKq6T9r6elSVYz8MUFUN0BDc3B4zDhMmkRL02B1H3Nwek2SSubLFIIgPxekSwWdnlZj+/TcbBbX2FLj/xKrYBnGaMYk+GwR5VHb9uJNure8z8ONjT/6dUfjQhpXHJfh854U2//yddQCqjsk4TBgGCb//6iJ3O96MVbffjlu2dEq2MWuBPQmOC3kel2J7vwdwsAOv6pgkUz2+zUFAxdYZhwlBnNIsWTimTpgoz2t/05IS3jxfZ2sYMgpTKrbOm+frpPkX3y7+RaIw+FPgfgM713C5sd6n4piMgpj3NwZIBG+erx9yv48zzEmYMlexj41x4eEJslbZ4ocf7pDlOa2SRatso2nMXN/7E4mX22X+5BsX+PntLlvDgPmqw++/qph5PS9mbxyx2vNwTJ3a9D0Nw+TIWB6Oz0sc3MAECVGqKginodge3DTPNVxubo+wDY1JKLnYKrHR99E1cYggdLldnm2Q+8o815dqs2v6cUrJ/s12eguDPwXujy11TXChWVI18Z0xFcc8Mrt+XExacYwjXfL9E/xhCbKBH3Nnb8LVdpmeF9HzE4ZhwrevtmdGZ2iCX28M+cmtPb52qcGbFxpcbpePpMNenSvz9r0+hgDbUBp6EsHVdvnITPZReYkffrjD1Xb50AY2X3VmxKPjcJJKx/2bpqYJPt4eY+hwbeGzjazjRbTLNpfbZaqOSlTuX/N5bOgpDP4UOCph9bsvzdMoWUipiC772XW4P0HLine truncated
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"m = dist.Normal(torch.zeros(2), scale=torch.ones(2))\n",
"s = m.sample((1000,))\n",
"\n",
"plt.scatter(s[:, 0], s[:, 1], alpha=0.2)\n",
"plt.axis('scaled');"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" torch.exp(x) + 2\n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Binary file not shown.

After

Width:  |  Height:  |  Size: 7.0 KiB

@@ -0,0 +1,3 @@
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg xmlns="http://www.w3.org/2000/svg" style="background-color: rgb(255, 255, 255);" xmlns:xlink="http://www.w3.org/1999/xlink" version="1.1" width="383px" height="207px" viewBox="-0.5 -0.5 383 207"><defs/><g><path d="M 1 172 L 161 172" fill="none" stroke="#000000" stroke-width="2" stroke-miterlimit="10" pointer-events="none"/><path d="M 41 172 L 41 98.37" fill="none" stroke="#000000" stroke-miterlimit="10" stroke-dasharray="3 3" pointer-events="none"/><path d="M 41 93.12 L 44.5 100.12 L 41 98.37 L 37.5 100.12 Z" fill="#000000" stroke="#000000" stroke-miterlimit="10" pointer-events="none"/><path d="M 121 172 L 139.46 98.18" fill="none" stroke="#000000" stroke-miterlimit="10" stroke-dasharray="3 3" pointer-events="none"/><path d="M 140.73 93.08 L 142.43 100.72 L 139.46 98.18 L 135.64 99.03 Z" fill="#000000" stroke="#000000" stroke-miterlimit="10" pointer-events="none"/><path d="M 1 92 L 161 92" fill="none" stroke="#6c8ebf" stroke-width="2" stroke-miterlimit="10" pointer-events="none"/><g transform="translate(21.5,185.5)"><switch><foreignObject style="overflow:visible;" pointer-events="all" width="39" height="12" requiredFeatures="http://www.w3.org/TR/SVG11/feature#Extensibility"><div xmlns="http://www.w3.org/1999/xhtml" style="display: inline-block; font-size: 12px; font-family: Helvetica; color: rgb(0, 0, 0); line-height: 1.2; vertical-align: top; width: 40px; white-space: nowrap; overflow-wrap: normal; text-align: center;"><div xmlns="http://www.w3.org/1999/xhtml" style="display:inline-block;text-align:inherit;text-decoration:inherit;white-space:normal;"><span class="MathJax_Preview" style="color: inherit;"></span><div class="MathJax_Display" style="text-align: center;"><span class="MathJax" id="MathJax-Element-28-Frame" tabindex="0" style="text-align: center;"><nobr><span class="math" id="MathJax-Span-147" style="width: 0.814em; display: inline-block;"><span style="display: inline-block; position: relative; width: 0.595em; height: 0px; font-size: 132%;"><span style="position: absolute; clip: rect(1.56em, 1000.55em, 2.338em, -1000em); top: -2.165em; left: 0em;"><span class="mrow" id="MathJax-Span-148"><span class="mi" id="MathJax-Span-149" style="font-family: MathJax_Math; font-style: italic;">x</span></span><span style="display: inline-block; width: 0px; height: 2.165em;"></span></span></span><span style="display: inline-block; overflow: hidden; vertical-align: -0.086em; border-left: 0px solid; width: 0px; height: 0.741em;"></span></span></nobr></span></div><script type="math/tex; mode=display" id="MathJax-Element-28">x</script></div></div></foreignObject><text x="20" y="12" fill="#000000" text-anchor="middle" font-size="12px" font-family="Helvetica"><span xmlns="http://www.w3.org/1999/xhtml" class="MathJax_Preview" style="color: inherit;"></span><div xmlns="http://www.w3.org/1999/xhtml" class="MathJax_Display" style="text-align: center;"><span class="MathJax" id="MathJax-Element-29-Frame" tabindex="0" style="text-align: center;"><nobr><span class="math" id="MathJax-Span-150" style="width: 0.814em; display: inline-block;"><span style="display: inline-block; position: relative; width: 0.595em; height: 0px; font-size: 132%;"><span style="position: absolute; clip: rect(1.56em, 1000.55em, 2.338em, -1000em); top: -2.165em; left: 0em;"><span class="mrow" id="MathJax-Span-151"><span class="mi" id="MathJax-Span-152" style="font-family: MathJax_Math; font-style: italic;">x</span></span><span style="display: inline-block; width: 0px; height: 2.165em;"></span></span></span><span style="display: inline-block; overflow: hidden; vertical-align: -0.086em; border-left: 0px solid; width: 0px; height: 0.741em;"></span></span></nobr></span></div><script xmlns="http://www.w3.org/1999/xhtml" type="math/tex; mode=display" id="MathJax-Element-29">x</script></text></switch></g><g transform="translate(101.5,178.5)"><switch><foreignObject style="overflow:visible;" pointer-events="all" width="38" height="27" requiredFeatures="http://www.w3.org/TR/SVG11/feature#Extensibility"><div xmlns="http://www.w3.org/1999/xhtml" style="display: inline-block; font-size: 12px; font-family: Helvetica; color: rgb(0, 0, 0); line-height: 1.2; vertical-align: top; width: 38px; white-space: nowrap; overflow-wrap: normal; text-align: center;"><div xmlns="http://www.w3.org/1999/xhtml" style="display:inline-block;text-align:inherit;text-decoration:inherit;white-space:normal;"><span class="MathJax_Preview" style="color: inherit;"></span><div class="MathJax_Display" style="text-align: center;"><span class="MathJax" id="MathJax-Element-30-Frame" tabindex="0" style="text-align: center;"><nobr><span class="math" id="MathJax-Span-153" style="width: 3.899em; display: inline-block;"><span style="display: inline-block; position: relative; width: 2.922em; height: 0px; font-size: 132%;"><span style="position: absolute; clip: rect(1.362em, 1002.87em, 2.463em, -1000em); top: -2.219em; left: 0em;"><span class="mrow" id="MathJax-Span-154"><span class="mi" id="MathJax-Line truncated

After

Width:  |  Height:  |  Size: 34 KiB