working 2d change of variables
This commit is contained in:
ritchie46 committed
1 parent
ef1043675c
commit
631e96e41a
1 file changed
+166
-28
@@ -2,7 +2,7 @@
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
@@ -28,7 +28,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 221,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -37,7 +37,7 @@
|
||||
"1.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 221,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
@@ -75,7 +75,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 222,
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -84,7 +84,7 @@
|
||||
"5.641895835477563"
|
||||
]
|
||||
},
|
||||
"execution_count": 222,
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
@@ -137,15 +137,14 @@
|
||||
"\\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(y) &=& P(f^{-1}(y))\\frac{\\text{d}f^{-1}(y)}{\\text{d}y} \\\\\n",
|
||||
" &=& P(f^{-1}(y))f'^{-1}(y)\n",
|
||||
"\\end{eqnarray} $$\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 223,
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -154,7 +153,7 @@
|
||||
"0.9987379589284238"
|
||||
]
|
||||
},
|
||||
"execution_count": 223,
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
@@ -194,16 +193,16 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 224,
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.9987894616589735"
|
||||
"0.9987894616589734"
|
||||
]
|
||||
},
|
||||
"execution_count": 224,
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
@@ -230,7 +229,7 @@
|
||||
" return y**2\n",
|
||||
"\n",
|
||||
"def f_i_prime(y):\n",
|
||||
" return 2*x\n",
|
||||
" return 2*y\n",
|
||||
"\n",
|
||||
"assert np.allclose(f_i(f(x)), x)\n",
|
||||
"\n",
|
||||
@@ -262,7 +261,7 @@
|
||||
"\\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",
|
||||
"The intuition of taking the modulus is that if the (local) rate of change of $x + dx = y + dy > 0$, then 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",
|
||||
"\n",
|
||||
"\n",
|
||||
@@ -272,7 +271,7 @@
|
||||
"$$ \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",
|
||||
"P(y)&=&P(f^{-1}(y)) \\cdot \\left| \\det \\mathbf{J}_{f^{-1}(y)} \\right| \\\\\n",
|
||||
"\\end{eqnarray}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
@@ -281,12 +280,29 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"execution_count": 259,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"(50, 50, 2)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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": [
|
||||
"<matplotlib.contour.QuadContourSet at 0x7f5f90e20ac8>"
|
||||
]
|
||||
},
|
||||
"execution_count": 259,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbL0lEQVR4nO3d/48cd53n8efLyTiOsOOc5HhtxfkicdFFBxKQs0y+SCsfghUxFvklkbISoM0vViJOAt0ipPADK+4PQMDuJb4RrJYIWIQOkosiJ1pLtwj4IUGx1yFkbZ3CriFzts+QFXasRLZC3vfDdJN2u2e6uupTVZ+qej2kkWema6Y/U55+1qc/XdOtiMDMzPpvQ9sDMDOzZjj4ZmYD4eCbmQ2Eg29mNhAOvpnZQDj4ZmYDUSj4kk5KelnSMUkvzrhckr4h6VVJv5B0R/qhmplZFVcvsO1/jojfrXHZvcBto7cPA4+P/jUzs0ykWtK5D3giVj0PXC9pZ6LvbWZmCRSd4QfwD5IC+B8RsTx1+Y3AaxMfr4w+d3pyI0kHgAMAV3H1f3rP1f+u1KDNzIbq/Nu//V1E3FDma4sG/56IOCVpO3BY0omI+MnE5ZrxNVc8Z8PoQLEMsHVpe9y97YGFB2xmNmTPnXns12W/ttCSTkScGv17FngS2DO1yQpw08THu4BTZQdlZmbpzQ2+pPdI2jJ+H/gz4JdTmz0NfGZ0ts6dwLmIOI2ZmWWjyJLOnwBPShpv/72IeE7SwwARcRA4BOwDXgXeBB6qZ7hmZlbW3OBHxL8AH5jx+YMT7wfw2bRDMzOzlPyXtmZmA+Hgm5kNhINvZjYQDr6Z2UA4+GZmA+Hgm5kNhINvZjYQDr6Z2UA4+GZmA+Hgm5kNhINvZjYQDr6Z2UA4+GZmA+Hgm5kNhINvZjYQDr6Z2UA4+GZmA+Hgm5kNROHgS7pK0j9JembGZXslnZN0bPT25bTDNDOzqoq8iPnY54DjwHVrXP7TiNhffUhmZlaHQjN8SbuATwDfrHc4ZmZWl6JLOl8Dvgi8s842d0l6SdKzkt5XfWhmZpbS3OBL2g+cjYgj62x2FLglIj4A/DXw1Brf64CkFyW9eOmdt0oN2MzMyikyw78H+KSkk8D3gY9I+s7kBhFxPiIujN4/BCxJ2jb9jSJiOSJ2R8TujRuurT56MzMrbO6DthHxKPAorJ6NA3whIj41uY2kHcD/i4iQtIfVA8nr6YdrVs2l23cl+14bT6wk+15mTVjkLJ3LSHoYICIOAvcDj0h6G3gLeDAiIs0QzdaWMuB1X7cPENY2tdXlrUvb4+5tD7Ry3dYtbUa9ST4gWBHPnXnsSETsLvO1pWf4ZqkNJexrmfXz+yBgKTn41oqhx70oHwQsJQffGuHApzO9L30AsKIcfKuFA98cHwCsKAffknHk8+ADgK3FwbfSHPhumPx/cvyHzcG3hTn03eX4D5uDb4U48v0z/j91+IfDwbc1OfLD4Fn/cDj4dhlHftgc/35z8A1w6O1KXvLpHwd/4Bx6m8fh7w8Hf6D6Hvpz772m1evf+quLrV5/HRz+7nPwB6YvoW876POsN76uHwwc/u5y8Aeiq6HPPexlrPUzde1A4PB3j4Pfc10LfR8DX9T0z96VA4DD3x0Ofk91JfRDDvw8XTsAOPz5c/B7JvfQO/DldeUA4PDny8HvkVxj78jXY3K/5hj/S7fvcvQzUzj4kq4CXgT+b0Tsn7pMwNeBfcCbwF9ExNGUA7W15Rh6R75Zucbfs/28LDLD/xxwHLhuxmX3AreN3j4MPD7612qWU+wd+TzkGH+HPw8bimwkaRfwCeCba2xyH/BErHoeuF7SzkRjtBku3b4rm9ife+81jn2mcvu/yeV3dqiKzvC/BnwR2LLG5TcCr018vDL63OnJjSQdAA4AbNqweaGB2rtyudHkFJL1vHGLar+OLb+O2q+jipxm/Z7tt2du8CXtB85GxBFJe9fabMbnrrgFRMQysAywdWl73reQDOUQ+hwj30TQq44hpwPC+P8wh/A7+s0qMsO/B/ikpH3AJuA6Sd+JiE9NbLMC3DTx8S7gVLphWtuxzyX0OcS9jFnjbvsgkEP4Hf1mzQ1+RDwKPAowmuF/YSr2AE8D/0XS91l9sPZcRJzGkmgz9m2GvqtxLyqXg0Db4fcST3NKn4cv6WGAiDgIHGL1lMxXWT0t86Ekoxu4IYa+75GfZ/Lnbzr+OYTf0a+XItq5W7l1aXvcve2BVq67C9qKfdOhH3rgF9H0AaCt8Dv663vuzGNHImJ3ma8tdFqmNWsIsX/jFjn2Cxrvs6b2W1v38tp+vKrP/NQKmWnjl72pG7YDn854X9Y9629rmcfr+vVw8DPSx9g78vVqas2/zfA7+uk4+JloOvZDCP3Fmy81fp3X/GZj49c51sSs/9x7r3H0O8zBb1nfZvVthL6NsK9lrbE0eSCoO/xtzPYd/TQc/Bb1aVbfVOhzivsiZo277oNAE+F39LvFwR+IumLfROi7Gvl5Jn+uOuNfZ/gd/W5x8FvS5Oy+jtjXHfq+Rn4t0z9vHQeAusLf9BKPo1+eg9+CpmLftVn90CK/nvG+6Fr4235CNluf//CqYV2OfR1/9HPx5kt/fLMr1bl/6jhwN/U3Hf7jrHIc/AZ1NfZ1ht6Kq2Of1fF/6+jny0s6Dely7FNqI/K37vpt7ddxcuWG2q9jrI7lnjduUdIlnqaWd7yevxgHvwGOfTOhbyLsi153nQeC1OF39PvPwe+JlLHvSujbDHxR02Os4wCQMvypH9B19PPi4Nesidl9jrGvK/RdiPx66jwApA5/16Jv8/lB2xo59mncuuu3f3zrmzp+rlQP7qa8p9fEA7l+EHc+z/A7LLfYpwx9H+O+nsmfN9Ws/+LNlyrP9rs20/fSzvo8w69J3bONvsa+rzP5RaTcBylm+12b6dva5gZf0iZJP5f0kqRXJH1lxjZ7JZ2TdGz09uV6htsNXbprWfXGnGr5wKG/UurwV5HyfP26o9+l21/TiizpXAQ+EhEXJC0BP5P0bEQ8P7XdTyNif/oh2rRUN5gUsa+qych/bMeJ5N/z8Jnbk3/PaeN9VHWpJ8UST1d4aWe2ucGP1Vc5vzD6cGn01s4rn3dAV5Zy2o593aGvI+5Fr6eug0CK8FeNfqo1fZ+5045CD9pKugo4Avx74L9HxAszNrtL0kvAKeALEfHKjO9zADgAsGnD5tKDzpVjX0xdsW8q8vNMjyP1AeDWXb919AvwLP9KhYIfEX8APijpeuBJSe+PiF9ObHIUuGW07LMPeAq4bcb3WQaWAbYubfe9hBa0Gfs6Qp9L5NczOcZU8a862696zn5Xom+XW+gsnYj4PfBj4ONTnz8fERdG7x8CliRtSzXILujC7L4vsf/YjhN/fOua1OOuul+r/J/m8LrF8/gB3MsVOUvnhtHMHknXAh8FTkxts0OSRu/vGX3f19MPd5hyOJUth9h3NfKzpPxZqp7N03b0c/j9HooiSzo7gW+P1vE3AD+IiGckPQwQEQeB+4FHJL0NvAU8OHqwdxC6MIuocsNsO/Z9ifws458txVJP1bX9vvJa/rvUVpe3Lm2Pu7c90Mp1p1Zn8Nteyikbe4e+nBThLxv9Kg/kpljPr3Mtv0/Bf+7MY0ciYneZr/Vf2vZcF2Pfp6WbRaX42cvu/7aXdurUhXvhTfBz6VSU++y+rDZj37T7rztaaLv/ef6Omkfyro/tOFFptl92eafNP87yGTv1c/B7LPdZ17S6Y1807It+fV0Hgq5FP/ULqFh6Dn4Fnt1frsrsvq7YV438oteROv5VH9Tt2gO5dc7y/eCtg99bZWf3Tce+jtA3Efki150y/lVn+4vyLL+f/KBthrp2XnIusb//uqOtxn5a6vGU3V9dexbSrv3+d4mD30NNz+7bllvop3U1+mV/H3J+7GjoZ+t4Saekof/ijLU5u8858tPGY23yTB+zaZ7hG1BuNufYLy7FuIcwy/eyTj0c/MxU/UXP+e40DDv2YymWoIb6h2lWjYNvpbT1QGDXY59Smeg3OcvP1ZCXYx38Evr2C9PUDbrqrLRvse/bzzMt93ubQ+Tg90ifb2B9jWMbSztdOU3T6/jpOfgZ6coveJlgVJnd9zX2Y134+fq2rDNUDv7A+YZsNhwOvlkGqszyfcaOFeXgW+28nJOfrqzjW1oOvi3EobBF5HoiQd/OtCuqyIuYb5L0c0kvSXpF0ldmbCNJ35D0qqRfSPLfjzcs1xtWFUOb3Q/t57XmFXkunYvARyLigqQl4GeSno2I5ye2uRe4bfT2YeDx0b9mZpaJuTP8WHVh9OHS6G36Ca/vA54Ybfs8cL2knWmHamZmVRRaw5d0laRjwFngcES8MLXJjcBrEx+vjD43/X0OSHpR0ouX3nmr7JjNzKyEQsGPiD9ExAeBXcAeSe+f2mTWAvIVL3sTEcsRsTsidm/ccO3iozUzs9IWOksnIn4P/Bj4+NRFK8BNEx/vAk5VGpmZmSVV5CydGyRdP3r/WuCjwPSJ1U8DnxmdrXMncC4iTicfra2pj68jOrQXCxnaz2vNK3KWzk7g25KuYvUA8YOIeEbSwwARcRA4BOwDXgXeBB6qabzWspMrN/hcfOu8jSdW2h5CK+YGPyJ+AXxoxucPTrwfwGfTDs3Muq6P9zy7zH9pa7U7fOb20l87lGWOKj9nmf17cuWG0tdn3eXgD9w1v9nY9hDMrCEOfka2/upi20MopMzs0LP8tXXh5/PEoB8c/B7p83ppF6JYRtWfq8/LOV2ZAHWJg1/CUB/hr6rKLB/6F/2+/TzT+jwB6SoH30rdXe/KLLHPqh5Ai+rbcs6QJ2wOfs/kPqtKMcvv+sw4xc9Qdj/6QD1sDn5m2lq3bHKWn2Jm2tXotznuMv9fZWf3VSceXr+vR5G/tLUZNp5YGeyr5uRiHM8uvHBIytA3tZRj/eMZfg+VnV1Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
@@ -298,30 +314,152 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"m = dist.Normal(torch.zeros(2), scale=torch.ones(2))\n",
|
||||
"s = m.sample((1000,))\n",
|
||||
"x1 = np.linspace(1, 5, num=50)\n",
|
||||
"x2 = np.linspace(1, 5, num=50)\n",
|
||||
"x1_s, x2_s = np.meshgrid(x1 ,x2)\n",
|
||||
"\n",
|
||||
"plt.scatter(s[:, 0], s[:, 1], alpha=0.2)\n",
|
||||
"plt.axis('scaled');"
|
||||
"base = stats.multivariate_normal(mean=[3, 3], cov=0.5)\n",
|
||||
"assert base.pdf([3, 3]) > base.pdf([1, 1])\n",
|
||||
"\n",
|
||||
"# field of all possible events X\n",
|
||||
"x_field = np.concatenate([x1_s[..., None], x2_s[..., None]], axis=-1)\n",
|
||||
"print(x_field.shape)\n",
|
||||
"\n",
|
||||
"plt.contourf(x1_s, x2_s, base.pdf(x_field))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Jacobian matrix.\n",
|
||||
"\n",
|
||||
"Assume the function $f: \\mathbb{R}^n \\mapsto \\mathbb{R}^m$. The jacobian is a 2D matrix, that stores on the partial derivatives of all the outputs and all the inputs. The width of the jacobian is equal to the functions input $n$. The height of the jacobian is equal to the functions output $m$.\n",
|
||||
"\n",
|
||||
"$$ \\mathbf{J} = \\begin{bmatrix}\n",
|
||||
"\\frac{\\partial f_1}{\\partial x_1} & \\dots & \\frac{\\partial f_1}{\\partial x_n} \\\\[6pt]\n",
|
||||
"\\vdots & \\ddots & \\vdots \\\\[6pt]\n",
|
||||
"\\frac{\\partial f_m}{\\partial x_1} & \\dots & \\frac{\\partial f_m}{\\partial x_n} \\\\[6pt]\n",
|
||||
"\\end{bmatrix}$$\n",
|
||||
"\n",
|
||||
"We will map 2D events to another domain. So we will have two inputs $n$ and two outputs $m$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 23,
|
||||
"execution_count": 275,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
|
||||
" import sys\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def f(x1, x2):\n",
|
||||
" return torch.exp(x1 / 3), x2**2\n",
|
||||
"\n",
|
||||
"def f_i(y1, y2):\n",
|
||||
" return 3 * torch.log(y1), y2**0.5\n",
|
||||
"\n",
|
||||
"x_field = torch.tensor(x_field)\n",
|
||||
"# Transform x events to y events\n",
|
||||
"y_field = np.concatenate(f(x_field[..., 0, None], x_field[..., 1, None]), axis=-1)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 276,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def f(x):\n",
|
||||
" torch.exp(x) + 2\n",
|
||||
" \n"
|
||||
"def create_det_jac(y_field):\n",
|
||||
" # create for every y1, y2 combination the determinant of the jacobian f_i(y1, y2)\n",
|
||||
" det_jac = np.zeros((y_field.shape[0], y_field.shape[1]))\n",
|
||||
" for i in range(y_field.shape[0]):\n",
|
||||
" for j in range(y_field.shape[1]):\n",
|
||||
" y_field = torch.tensor(y_field)\n",
|
||||
" y_field.requires_grad_(True);\n",
|
||||
" fiy = torch.cat(f_i(y_field[..., 0, None], y_field[..., 1, None]), dim=-1)\n",
|
||||
" fiy[i, j].sum().backward()\n",
|
||||
" \n",
|
||||
" # Ouputs of the partial derivatives are independent. I.e. f1 is dependent of y1 and not y2\n",
|
||||
" # and vice versa f2 is dependent of y2 and not y1\n",
|
||||
" # therefore the multiplication w/ 0 \n",
|
||||
" row1 = y_field.grad[i, j].data.numpy() * np.array([1., 0.])\n",
|
||||
" row2 = y_field.grad[i, j].data.numpy() * np.array([0., 1.])\n",
|
||||
" \n",
|
||||
" det = np.linalg.det(np.array([row1, row2]))\n",
|
||||
" det_jac[i, j] = det\n",
|
||||
" return det_jac\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 277,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
|
||||
" \"\"\"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.contour.QuadContourSet at 0x7f5f90d534a8>"
|
||||
]
|
||||
},
|
||||
"execution_count": 277,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD8CAYAAABuHP8oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWAElEQVR4nO3df6ik1X3H8c/XuGuCGmNYfyyuuiCh2zQQYxexEYJtakjTUg0kJUJT2wY2LSkYEiipfzS2pZA/0rS0kBZTQyxNLKFqEoKRLNZWAm3Kam007NIkZbUbV1ebVlcqLuq3f9xn4uzs/Hh+nPOcH8/7BYv3zsy9c+7xzvueOfPMjLm7AADlOS31AAAA/RBwACgUAQeAQhFwACgUAQeAQhFwACjUxoCb2cVmdr+ZHTSz75rZTc3pt5jZD83s4ebfe+IPFwAwY5uOAzeznZJ2uvtDZna2pAclXS/pVyQ97+6fjj9MAMCi0zddwN2PSjrafHzczA5Kuij2wAAA621cgZ90YbPdkh6Q9BZJH5P065Kek3RA0sfd/X+WfM0+Sfsk6TU6/afPPP3coWMGgEl57qWnn3H38xZPbx1wMztL0j9J+mN3v8vMLpD0jCSX9Efa2mb5zXXf45xt5/vbd7y/8+ABYMruffKzD7r73sXTWx2FYmbbJN0p6YvufpckuftT7v6yu78i6XOSrgw5YADAem2OQjFJt0k66O6fmTt959zF3ivp0fDDAwCssvFBTElXS/qgpEfM7OHmtJsl3WBml2trC+WwpA9HGSEAYKk2R6F8S5ItOeue8MMBALTFMzEBoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKRcABoFAEHAAKtTHgZnaxmd1vZgfN7LtmdlNz+hvNbL+Zfa/577nxhwsAmGmzAn9J0sfd/SclXSXpI2b2ZkmfkHSfu79J0n3N5wCAkWwMuLsfdfeHmo+PSzoo6SJJ10m6vbnY7ZKujzVIAMCpOu2Bm9luSW+T9G1JF7j7UWkr8pLODz04AMBqrQNuZmdJulPSR939uQ5ft8/MDpjZgROvvNBnjACAJVoF3My2aSveX3T3u5qTnzKznc35OyUdW/a17n6ru+91973bT3tdiDEDANTuKBSTdJukg+7+mbmzvibpxubjGyV9NfzwAACrnN7iMldL+qCkR8zs4ea0myV9StKXzexDkh6X9P44QwQALLMx4O7+LUm24ux3hh0OAKAtnokJAIUi4ABQKAIOAIUi4ABQKAIOAIUi4ABQKAIOAIUi4ABQKAIOAIVq81R6FOLEnl2jXM/2Q0dGuR4A6xHwQo0V6y7XTdiBcRHwAqSMdReEHRgXAc9UKdFuY9nPQtSB4Qh4RmqK9ibzPysxB/oh4IlNKdqrLM4BQQfaIeCJEO7VWJ0D7RDwkRHubmbzRciBUxHwkRDuYViVA6ci4JER7vCIObCFgEdCuMdBzDFlvBZKBMQ7jRN7djH3mBQCHhgBSY+QYyrYQgmEYOSHI1hQOwIeAPHOG/vkqBUBH4Bwl4dVOWrCHnhPxLts7JOjBgS8B2749SDkKBkB74gbe50IOUpEwDvgBl4/Qo6S8CBmSyXfqJ+97IxeX3fOD14MPJJy8GAnSkDAW8g93n0DHeL71h75E3t2EXFki4BvkGO8YwW7j2VjqS3qrMaRKwK+Ri7xzinYbSyOt5agE3LkhoCvkDrepUV7ndqCzrYKckHAl0gZ75rCvcr8z1hqzFmNIwcEfEGKeE8h2quUHnNCjpQI+Jyx4z3lcC9TcszZVkEKG5/IY2afN7NjZvbo3Gm3mNkPzezh5t974g4zPuKdl2cvO6O4OeJJQBhbm2difkHSu5ec/qfufnnz756ww6pXiWFKaTZfJc0ZEcdYNm6huPsDZrY7/lDSGeMGV1KAcjWbwxK2V9gbxxiGvBbK75jZd5otlnNXXcjM9pnZATM7cOKVFwZcXRzEuzwlrchZjSOmvgH/S0mXSbpc0lFJf7Lqgu5+q7vvdfe92097Xc+riyP2jauk0JSolO0V9sYRS6+Au/tT7v6yu78i6XOSrgw7rPLlHpXalBJyIKRehxGa2U53P9p8+l5Jj667fI5i3phyCMnxSy3o9zv7MQ/6/WLJfZ+cww0R0saAm9kdkq6RtMPMjkj6pKRrzOxySS7psKQPRxxjcLHinSrcoWPd5TpyDfuzl52RdcQlHuDEcG2OQrlhycm3RRhL0caM9xjBbmtxLDkFndU4aje5Z2LGWH2PEe+cor1OjkHPOeREHENMKuClxbuUaK8z/zOkjnmuIWdLBX1NKuChxYp3DeFeJpfVea7746zG0dVk3tQ49Oo7RryPX2rVxnuZlD9vrocdcqghumAF3kPoG/6Uor1Mym2WHLdVWImjrUmswEOuakLGe2or7jZSzUluq3GevYk2JhHwUELHG6ulCHmO2ypEHOtUH/DcbgCsurtJFfKc5PY7jHywB95SiBt1qnC/eMmJXl93xuPbA4+kv9ncjbVHntuRKuyLY5mqAx5q5VJKvPuGus/3SxX3MUOe2wOcRByLqg54CEPjHTPcoYM95LrHDvrYISfiyFG1AQ+x+s4x3imjvU6qoI8V8pxW40QcM9U/iJlKyHi/eMmJH/8rxdhjHuvxhVwe4OSBTUiVBjz16jtUTEqL9ipjxXysI1aIOHJRZcCHSh3vWsK9zFghj42IIwfVBTzlL/TQcNQc7kWxf9YxVuNEHKlVF/Ch+t4oQ8R7imJvr4wR8RxCTsSniYAHMCQSU1p1bxJrLqayGifi01PVYYRDf4HHvhGOGe7du57udPnDR86LNJLNXrzkRJTDEI9falEPN8zheHEOMZyWqgI+xNhbJ7Hi3TXUXb/PWGGfzU/okBNx1KSagKe4+9gn3qHDHSrYQ64vZtRjhDz2k3+IOMbCHrj6rb5Tx3v3rqdHj/cqs7HEHE+M/fGY++LsiWMM1azAx5Qq3rkEe535McZYmYfeH4+5pZLDShx1q2IFPmSlMcZKaWi8c1ptdxFr3KzE22MVXrcqAj6mrjf2IbEpNdyLYvwcobdUiDhKNOmAd71hjR3v2sQKeSg1v1MSEa9T8QGv7RezllX3OqF/xpCr8VgRT70Kl+q7raCCgI9ljNV37eFeFCPkIdQccdRlskehxLwx5Rrvay881Opy+5/cE3kkJ5v97CGOWgl1lEqso1NSH5nC8eF1mWzAu+iyIssh3m1D3fXrY4d9966niTjQQdEB77unl9Nd2VDxHhrtPtcRI+ihVuNEfDVW4fUoOuBjiLn6HhrvMaLd9vpDxzzEajz3iKdExOvAg5iBjBnvay88lDzei2KMKcS9k5wf2MzpniDKNLmAd7nRxDoaYWi8cxY65EQ8Hg4rLF+xWyg5/fJ1CUTfIOUe7kWz8YbYWgmxLx7rNcZLx1ZK2Sa3Ak9pKvGeF3JFPnQ1HmIlXtsqHGXbGHAz+7yZHTOzR+dOe6OZ7Tez7zX/PTfuMMOIsX0S8111ctzr7ouIr8dWCvposwL/gqR3L5z2CUn3ufubJN3XfI41uoanlnDPC/UHaWrPWAVW2Rhwd39A0o8WTr5O0u3Nx7dLuj7wuNaKvWIIvfom3idLHXFW4adiFV6mvnvgF7j7UUlq/nv+qgua2T4zO2BmB0688kLPq5uO2uM9k/rnzDXiQBfRj0Jx91sl3SpJ52w7P9mzIUKvbmKtvkN73+sfGvT1f//cFYFGcqprLzw06CiVoU/2yfHIFJ6hiS76BvwpM9vp7kfNbKekYyEHlVLKVVWoVenQaK/6XjFinjriQ9X4LE2Uo+8Wytck3dh8fKOkr4YZTl26rL5DxPt9r38oaLzH+v5Df/bU++GhsReOttocRniHpH+W9BNmdsTMPiTpU5KuNbPvSbq2+XwyQt/oQ8V7LDFCnnpPfAj2wpHKxi0Ud79hxVnvDDyWqoy19z1muFddd6itlSHbKUO2UtgLR6mKeyZmn7t4be+SplhJDVl5poz3vJDjGDIfKR8wrmkVzjZKOYoLeGq57JnmEu+Z3MbTVS7/X4EuCHgEbVeCfVebucYy1LhKXYWHxmukYBMCXphc4z2T+/jWGboKr2kbBWUg4AUpJY4hxlnyUSk1YB+8DAQ8sNjbJ6VI+cempm0UYB0C3gEPdI2rxD9yobdR2AfHOgS8kfv+ZSnbJ/NKHDN/pFESAo4qsY2CKag+4DneBS1xa6CvoavwKc0V0FX1AQeAWhHwApS4lwwgPgIOAIUi4ABQKAIOAIUi4AWI+b6UAMpVfcBzfFH8Ie8BWZqhf3ymNFdAV9UHHNOU8o2OgbEQcAAoFAFvnP2Yb7xMyvdNLHEfvMTtk6H/j9v8HnWR4xYg8kHAA2t7171PnEqKeMqxsn2CqSDgCC5EvHnwMq3th46kHgJaIOCFyX0Vnvv41slt+wTYhIB31OZGHnMbRdqKZI6hDDWmIavvmrZP2P/GJpMIeNsbQooV1JBY5RLxkH9QUsWb1fer2D4pR3EBL+WXa6yVYOqIh7z+VPveKY8uWoXVN9o4PfUASnTG49uDvvXW/if3DHrjgllEx3rZ2Rh/NIbGO+XWSU2rb5SFgC84+zEP9v6Yh4+c1/qtvYZGXDo5rKFjHnOlnzLerL5PVso9XGyZTMDP+cGLQd9eLfQqXAoT8ZllwW0b9TG3ZUqPN6tvpDSZgKfSZRUuhY34otT75fNC7HfXGG9W3+iiuAcxpfi/aG1vmG0j0DU0tT+JJXW8c0W80dWkVuCht1G66LMSl+p6V/Ycwh1qz5vXPEEOilyBjyH0KlzqF5/9T+4pfkULine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"px = base.pdf(x_field)\n",
|
||||
"transformed = px * np.abs(create_det_jac(y_field))\n",
|
||||
"plt.contourf(y_field[..., 0], y_field[..., 1], transformed)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 278,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.9907110850291531"
|
||||
]
|
||||
},
|
||||
"execution_count": 278,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"np.trapz(np.trapz(transformed, y_field[:, 0, 1], axis=0), y_field[0, :, 0])"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
||||
Reference in new issue
Block a user