{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "- Currently using a pretty naive minimization algorithm. It will approximate the gradient and move in that direction. This requires lots of function calls. \n", "- We have no way to get the Hessian - which is needed for the variance estimates of our parameters. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Automatic differentiation - how to differentiate an algorithm\n", "\n", "Symbolic differentiation is what we learned in school, and software like SymPy, Wolfram and Mathematica do this well. But I want to differentiate the following:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5b348c83O4EACQkBSSBhE1EUMIDsCGix14p70S7gUtqrttpWW3vtdv3VXrdetdWqVBHqAi6tilYEZJUiS5B9D2ENSQh7IOsk398fc/COMQHinMnkJN/36zWvnDnnOTPfnDyZ7zznec55RFUxxhjTfEWEOwBjjDHhZYnAGGOaOUsExhjTzFkiMMaYZs4SgTHGNHNR4Q7g60hOTtaMjIxwh2GMMZ6yevXqQ6qaUnO9JxNBRkYG2dnZ4Q7DGGM8RUT21LbeTg0ZY0wzZ4nAGGOaOUsExhjTzFkiMMaYZs4SgTHGNHOuJAIRmSoiB0VkYx3bRUT+LCI5IrJeRPoHbJsoIjucx0Q34jHGGHPu3GoRTAPGnWH7VUAP5zEZeB5ARJKA3wGDgIHA70Qk0aWYjDHGnANXriNQ1SUiknGGIuOBv6v/ntfLRaStiHQERgHzVPUIgIjMw59QZrgRlzHNXbmviryjpeQdK+XIqQpOlPk4UVqJr0qpVkVViYmKoEVMFPExkbRtEU1yQiztWsbQoU0c8TGevNTI1FND/ZU7AfsCnu931tW1/itEZDL+1gSdO3cOTZTGeFhJhY+Vu46wfv9x1u8/zuYDx8k/UUYwU460axlDWlI8me3i6ZGaQPf2rejVIYHOSfGIiHvBm7DyTLpX1SnAFICsrCybTccYoKi4nH+tP8CCbUUszz1Mha8aEeia3JIBmUlkJrckPTGetMQWtGsVS+sWUbSOiyYmMgIREBEqfNWUVlRxqsLHsZJKDp0s59DJcvKPl7H/aAn7jpSyavdR3lt74Iv3TYiL4sLzWnNJWlv6d0nk0i6JJLeKDeORMMFoqESQB6QHPE9z1uXhPz0UuH5RA8VkjCdVVyuf5hxi5sq9zNtciK9a6ZrSku9d1oVR56fQN70tCXHR5/x6MVERxERF0CY+mvPatqiz3MlyHzkHT7I1/wQb8o6zMe84r/x7Ny8uyQUgM7kll3VNYnC3ZAZ3bUdKgiUGr2ioRDALuEdEZuLvGD6uqvkiMgf4Y0AH8ZXArxooJmM8RVVZsPUgT87dzpb8EyTGRzNpSAbfHpBOj9SEkL9/q9go+qa3pW96WyY468oqq9h04Dir9xxl5a4jfLgunxkr/Wd7e3dszcjzUxjZM4WsLolERdpo9cZK3JizWERm4P9mnwwU4h8JFA2gqi+I/2Tis/g7gkuA21Q129n3duC/nJd6RFVfOdv7ZWVlqd10zjQna/cd478/2MSavcfonBTPvWN6cPUlHYmNigx3aF/iq6pm04ETLM05xJLtRazecxRftdKmRTSjzk9h7AWpXN6rPa1iPXNWukkRkdWqmvWV9V6cvN4SgWkuyn1VPPPJDl5YvJOUhFjuHdOTm7LSiPbIt+viskqW7jjEJ1sOsmBrIUdLKomJimBEjxS+2acDY3un0roep7FMcOpKBJaWjWmkthac4N4Za9lWWMzNWWn8+urenvvQTIiL5qo+HbmqT0eqqpXVe44ye2M+H28s4JMthcRERXD5+Slcc0knxlzQnrjoxtXCaS6sRWBMI7R4exF3v/45LWIieeyGPozulRrukFxVXa2s2XeMD9cf4MP1+RQVl5MQG8VVfTpwff80BmYkERFhw1PdZqeGjPGIN1ft5b/e3UjP1ARemTSADm3iwh1SSFVVK8tzD/Pumjxmb8jnVEUV6UktuLF/Ojdc2om0xPhwh9hkWCIwxgOeW5jDE3O2MaJnCs/d2q9ew0CbgpIKH3M3FfL26n38O+cwIjCsezK3DuzM2N6pnukbaawsERjTyL2VvY9fvLOe8X3P48mbLmn2H3r7jpTwzur9vJW9j/zjZSS3iuXbA9K4dVAXOp3hegdTN0sExjRii7cXcfu0VQzu2o6pkwYQE9W8k0Cgqmpl8faDvLFiHwu2FgIwuld7vjc4g+Hdk60voR5s1JAxjdTGvOPc9dpqeqYm8Px3+1sSqCEyQhjdK5XRvVLJO1bKjBV7mblqL59sWUnX5JZMHJLBDZem2bUJQbAWgTFhVFLhY9zTn+Krqubdu4eS2rppdwy7pdxXxewNBUxbtpu1+46REBvFtwekM3FIBulJ1rlcF2sRGNMIPf7xNvYeKeHNyZdZEqiH2KhIru3XiWv7dWLtvmNMXbqLact2M/Xfuxh3UQd+MLwr/Trb1CbnyloExoTJ8tzDTJiynElDMvj9NReGOxzPyz9eyvRle3hjxR5OlPkYmJHED0Z0ZUyv9taP4LDOYmMakdOnhERg9r3DbQIYF50s9/HWqn28vHQXecdK6d6+FT8c0ZXxfTs1+/6XuhJB8z4qxoTJk3O2s/dICY/fcLElAZe1io3i9mGZLH5gFM9M6EtUhPDAO+sZ+cRCpi7dRUmFL9whNjqWCIxpYPnHS3l1+W5uGZjOoK7twh1OkxUVGcH4vp2Yfe9wpt02gM5J8Tz84WaGPbaQv8zfwfHSynCH2GjYVxFjGtiUJbmowl2juoc7lGZBRBh1fntGnd+e7N1HeG5hDn+at50pn+YyaUgGtw/NJLFlTLjDDCtrERjTgIqKy5mxci/X9utkwxzDICsjiVduG8iHPx7G0G7J/GVBDsMeW8Cjs7dy+GR5uMMLG0sExjSgl5fuotxXzV2juoU7lGbtok5teOF7lzL3pyO4vFd7Xlyyk+GPL+R/Zm/hyKmKcIfX4FxJBCIyTkS2iUiOiDxYy/anRGSt89guIscCtlUFbJvlRjzGNEbHSip49bPdXH3xeXRNaRXucAzQMzWBZ2/tz9z7RjD2glSmLMll+GMLeGLOVo6VNJ+EEPTwURGJBLYDVwD7gVXALaq6uY7yPwb6qertzvOTqlqv/wobPmq86OlPtvP0Jzv4+L7h9OrQOtzhmFrsKCzmmfk7+NeGfFrGRHHHsEzuGJ7puQmB6hLK4aMDgRxVzVXVCmAmMP4M5W8BZrjwvsZ4RlW18tryPYzp1d6SQCPWw2khfHzvCIZ1T+aZ+TsY/thCnluYw6nypjvs1I1E0AnYF/B8v7PuK0SkC5AJLAhYHSci2SKyXESuretNRGSyUy67qKjIhbCNaTjZu49w6GQF1/ar9V/DNDLnd0jghe9dyoc/HsalXRJ5Ys42Rj6xkGn/3kW5ryrc4bmuoTuLJwDvqGrgkeziNFVuBZ4WkVp70VR1iqpmqWpWSkpKQ8RqjGtmbyzwz8/bq324QzH1cFGnNkydNIB//OcQurdvxe8/2MzoJxfzdvY+qqq9d1eGuriRCPKA9IDnac662kygxmkhVc1zfuYCi4B+LsRkTKOhqszZVMCIHil2q2SPurRLIjN+cBmv3TGI5FYxPPDOesY9vYQ5mwrw4m16anIjEawCeohIpojE4P+w/8roHxHpBSQCnwWsSxSRWGc5GRgK1NrJbIxXrdt/nPzjZYy7qEO4QzFBEBGG9UjmvbuH8vx3+lOlyg9fXc31zy9j5a4j4Q4vKEEnAlX1AfcAc4AtwFuquklEHhaRawKKTgBm6pfT5wVAtoisAxYCj9Y12sgYr5q9MZ+oCOGKC1LDHYpxgYhwVZ+OzL1vBI9e34cDx0q5+cXPuGPaKrYXFoc7vK/F7j5qTAipKpc/uYj0pHhevWNQuMMxIVBaUcW0Zbv56yL/yKKbLk3np1f0pEObxje/hN191Jgw2FpQzO7DJVx1Ucdwh2JCpEVMJP85qhtLHric24dm8u6aPEY9uZA/zd3GSY8MObVEYEwIzd5YgAhc0dtOCzV1iS1j+PXVvZn/85Fc2bsDf1mQw6gnFvLa8j34qqrDHd4ZWSIwJoTmbipgQEYSKQmx4Q7FNJD0pHj+fEs/3r97KF1TWvHr9zbyjaeXsGBrYaMdYWSJwJgQKa2oYlthMUO62ZwDzdEl6W15c/JlTPnepajC7dOy+e7LK9h84ES4Q/sKSwTGhEjOwZOowvmpCeEOxYSJiHDlhR2Y89MR/P5bvdl04AT/8ZdP+eU76zlYXBbu8L5gicCYEDk9lLBnB0sEzV10ZASThmay+H5/h/I/1+zn8icW8dzCHMoqw3/LCksExoTI9sJiYiIj6GIT0BhHm/hofnN1b+b+dCRDuifzxJxtjPnTYv61Pj+s/QeWCIwJkW2FxXRr34qoSPs3M1+WmdySv30/izfuHERCXBR3v/E5335xORvzjoclHquhxoTIjsKT9Ey1CWhM3YZ0T+ZfPxnOH6/rQ07RSb717FIe/Md6DjXwtJmWCIwJgeKySvKOldLTOorNWURGCLcO6szC+0dxx9BM3lnt7z/425JcKnwNc/2BJQJjQmDHwZMAlgjMOWvTIppfX92bOT8dwaUZiTzy0RbGPbOERdsOhvy9LREYEwLbC/wjhmzoqKmvbimtmHbbQKZOykIVJr2yijunr2LP4VMhe09LBMaEwLbCYlpER5KW2CLcoRiPGt0rlY/vG86DV/Xis52HueJ/l/DEnK2UVLh//yJLBMaEwI7Ck/RIbUVEhIQ7FONhsVGR/GhkNxbcP4r/uLgjLy7OZd+RUtffxxKBMSGwrbDY+geMa1Jbx/HUt/uy6IFRnB+CCxRdSQQiMk5EtolIjog8WMv2SSJSJCJrncedAdsmisgO5zHRjXiMCaejpyooKi63oaPGdWmJobk4MegJVEUkEngOuALYD6wSkVm1zDT2pqreU2PfJOB3QBagwGpn36PBxmVMuHxxawlrERiPcKNFMBDIUdVcVa0AZgLjz3HfbwDzVPWI8+E/DxjnQkzGhM12GzpqPMaNRNAJ2BfwfL+zrqYbRGS9iLwjIun13NcYz9heUExCbBQdG+FUhcbUpqE6iz8AMlT1Yvzf+qfX9wVEZLKIZItIdlFRkesBGuOWbYXF9OyQgIiNGDLe4EYiyAPSA56nOeu+oKqHVfX0zTNeAi49130DXmOKqmapalZKSooLYRsTGjsKi62j2HiKG4lgFdBDRDJFJAaYAMwKLCAigTN3XwNscZbnAFeKSKKIJAJXOuuM8aSqauVoSSWpre20kPGOoEcNqapPRO7B/wEeCUxV1U0i8jCQraqzgJ+IyDWADzgCTHL2PSIi/w9/MgF4WFWPBBuTMeFS7vNPMhIXHRnmSIw5d0EnAgBV/Qj4qMa63wYs/wr4VR37TgWmuhGHMeFWXum/W2RslF2rabzDaqsxLir3nU4E1iIw3mGJwBgXnZ5/Ni7a/rWMd1htNcZF1iIwXmSJwBgXne4stj4C4yVWW41xUZnTWWyjhoyXWCIwxkVftAisj8B4iNVWY1xkw0eNF1ltNcZFZXZBmfEgSwTGuMhaBMaLrLYa46KyL0YNWYvAeIclAmNcVP7FqCH71zLeYbXVGBfZBWXGiywRGOOi07eYsD4C4yVWW41xUbmvmpjICCIibHYy4x2WCIxxUbmvyloDxnOsxhrjorLKamLtGgLjMa4kAhEZJyLbRCRHRB6sZfvPRGSziKwXkfki0iVgW5WIrHUes2rua4yXWIvAeFHQM5SJSCTwHHAFsB9YJSKzVHVzQLE1QJaqlojIfwKPA992tpWqat9g4zCmMSj3Vdt9hoznuFFjBwI5qpqrqhXATGB8YAFVXaiqJc7T5UCaC+9rTKNTXllFnA0dNR7jRiLoBOwLeL7fWVeXO4DZAc/jRCRbRJaLyLV17SQik51y2UVFRcFFbEyIWIvAeJErk9efKxH5LpAFjAxY3UVV80SkK7BARDao6s6a+6rqFGAKQFZWljZIwMbUU3lltfURGM9xo8bmAekBz9OcdV8iImOBh4BrVLX89HpVzXN+5gKLgH4uxGRMWJT5quzOo8Zz3EgEq4AeIpIpIjHABOBLo39EpB/wIv4kcDBgfaKIxDrLycBQILCT2RhPsRaB8aKgTw2pqk9E7gHmAJHAVFXdJCIPA9mqOgt4AmgFvC0iAHtV9RrgAuBFEanGn5QerTHayBhP8Q8ftRaB8RZX+ghU9SPgoxrrfhuwPLaO/ZYBfdyIwZjGoKyy2u48ajzHaqwxLrIWgfEiSwTGuKjcZy0C4z1WY41xiapSVmktAuM9lgiMcYmvWqlWm4vAeI/VWGNccnp2MruOwHiNJQJjXPLF7GTWR2A8xmqsMS75v/mK7d/KeIvVWGNcUu60COzUkPEaSwTGuKSs0loExpusxhrjknKf00dgw0eNx1giMMYlX/QRWGex8Rirsca45ItRQ9YiMB5jicAYl9ioIeNVVmONcUmZjRoyHmWJwBiXWIvAeJXVWGNcYreYMF7lSiIQkXEisk1EckTkwVq2x4rIm872FSKSEbDtV876bSLyDTfiMSYcyu0WE8ajgq6xIhIJPAdcBfQGbhGR3jWK3QEcVdXuwFPAY86+vfHPcXwhMA74q/N6xniOnRoyXuVGjR0I5KhqrqpWADOB8TXKjAemO8vvAGPEP3nxeGCmqpar6i4gx3k9YzynvLIKEYiJtERgvMWNGtsJ2BfwfL+zrtYyquoDjgPtznFfAERksohki0h2UVGRC2Eb464yXzWxURH4v+MY4x2e+eqiqlNUNUtVs1JSUsIdjjFfUW6zkxmPciMR5AHpAc/TnHW1lhGRKKANcPgc9zXGE2y+YuNVbtTaVUAPEckUkRj8nb+zapSZBUx0lm8EFqiqOusnOKOKMoEewEoXYjKmwdl8xcarooJ9AVX1icg9wBwgEpiqqptE5GEgW1VnAS8Dr4pIDnAEf7LAKfcWsBnwAXeralWwMRkTDuVOH4ExXhN0IgBQ1Y+Aj2qs+23AchlwUx37PgI84kYcxoST/9SQtQiM99jXF2Nc4j81ZP9Sxnus1hrjknJftV1VbDzJaq0xLin3VRFnncXGgywRGOOSskprERhvslprjEvKfTZ81HiTJQJjXFJeaReUGW+yWmuMS+yCMuNVlgiMcYmNGjJeZbXWGBeoqnNlsbUIjPdYIjDGBTYpjfEyq7XGuMDmKzZeZonAGBeU+5z5iq1FYDzIaq0xLiivtFNDxrus1hrjgtMtAjs1ZLzIEoExLiizFoHxMKu1xrjgiz4CaxEYDwoqEYhIkojME5Edzs/EWsr0FZHPRGSTiKwXkW8HbJsmIrtEZK3z6BtMPMaEy+k+gjhrERgPCrbWPgjMV9UewHzneU0lwPdV9UJgHPC0iLQN2P6AqvZ1HmuDjMeYsCizFoHxsGATwXhgurM8Hbi2ZgFV3a6qO5zlA8BBICXI9zWmUbFRQ8bLgq21qaqa7ywXAKlnKiwiA4EYYGfA6kecU0ZPiUjsGfadLCLZIpJdVFQUZNjGuMsuKDNedtZEICKfiMjGWh7jA8upqgJ6htfpCLwK3Kaq1c7qXwG9gAFAEvDLuvZX1SmqmqWqWSkp1qAwjUtZpV1QZrwr6mwFVHVsXdtEpFBEOqpqvvNBf7COcq2BfwEPqerygNc+3ZooF5FXgPvrFb0xjYS1CIyXBfv1ZRYw0VmeCLxfs4CIxADvAn9X1XdqbOvo/BT8/Qsbg4zHmLCwW0wYLwu21j4KXCEiO4CxznNEJEtEXnLK3AyMACbVMkz0dRHZAGwAkoE/BBmPMWFhF5QZLzvrqaEzUdXDwJha1mcDdzrLrwGv1bH/6GDe35jGotxXRVSEEBVpicB4j9VaY1xQXlltrQHjWVZzjXFBma/KLiYznmWJwBgXlFdW2+0ljGdZzTXGBf6J661FYLzJEoExLiirrLI+AuNZVnONcYG1CIyXWSIwxgXlPmsRGO+ymmuMC8ps+KjxMKu5xrig3Fdt9xkynmWJwBgX2Kkh42VWc41xgf/KYmsRGG+yRGCMC8p9VcRF27+T8Saruca4wFoExsssERjjgjJrERgPs5prTJCqqpXKKrUWgfGsoBKBiCSJyDwR2eH8TKyjXFXApDSzAtZnisgKEckRkTed2cyM8ZQKZ5rKWGsRGI8KtuY+CMxX1R7AfOd5bUpVta/zuCZg/WPAU6raHTgK3BFkPMY0uEMnywFo2yI6zJEY8/UEmwjGA9Od5en45x0+J848xaOB0/MY12t/YxqLHQeLAeiR2irMkRjz9QSbCFJVNd9ZLgBS6ygXJyLZIrJcRE5/2LcDjqmqz3m+H+gUZDzGNLhtBScB6JGaEOZIjPl6zjpnsYh8AnSoZdNDgU9UVUVE63iZLqqaJyJdgQXOhPXH6xOoiEwGJgN07ty5PrsaE1LbC4s5r00crePs1JDxprMmAlUdW9c2ESkUkY6qmi8iHYGDdbxGnvMzV0QWAf2AfwBtRSTKaRWkAXlniGMKMAUgKyurroRjTIPbXlhsrQHjacGeGpoFTHSWJwLv1ywgIokiEussJwNDgc2qqsBC4MYz7W9MY1ZVrew4eJLzO1giMN4VbCJ4FLhCRHYAY53niEiWiLzklLkAyBaRdfg/+B9V1c3Otl8CPxORHPx9Bi8HGc8ZlVVW4c8/xrhjz+FTVPiq6dHeOoqNd5311NCZqOphYEwt67OBO53lZUCfOvbPBQYGE0N9/Oa9jeQdK+V337rQvsEZV2wv9HcUW30yXtasroC5OL0tmw6c4Jt//pTfz9rE8ZLKcIdkPG57YTEi0N1aBMbDmlUi+N5lXVh0/ygmDEhn+me7ufxPi3hjxV6qqu10kfl6thcWk54YT3xMUI1rY8KqWSUCgMSWMTxyXR8+uGcY3VJa8l/vbuCaZ5eyaveRcIdmPGh7YTE9bcSQ8bhmlwhOu6hTG9764WD+fEs/jpyq4KYXPuPHM9aQd6w03KEZj6jwVZNbdIqedkWx8bhmmwgARIRrLjmP+T8fyU/G9GDupgLG/GkRT83bTmlFVbjDM43c7sOn8FWrdRQbz2vWieC0+JgofnZFT+b/fCRjLkjlmfk7GP2nRby/Ns+Gm5o6bStw7jHU3hKB8TZLBAHSEuN57tb+vPXDwbRrFcO9M9dy/fPL+Hzv0XCHZhqh7YXFREYIXVNahjsUY4JiiaAWAzOTmHX3MB6/8WL2Hy3l+r8u4yfWf2Bq2F5YTEa7eOKibUIa422WCOoQESHcnJXOovtH8ePR3ZmzqYDRTy7i8Y+3Ulxm1x8Y/8VkNmLINAWWCM6iZWwUP7/yfBbcP4qrLurAXxft5PInF/H6ij34qqrDHZ4Jk7LKKvYcPmWJwDQJlgjOUae2LXh6Qj/ev3somckteejdjVz1zKcs3HrQOpSboc/3HKVa4cLzWoc7FGOCZomgni5Jb8tbPxzMC9+9FF+1ctu0VXznpRVszKvX9ArG4z7eVEBcdATDeiSHOxRjgmaJ4GsQEcZd1IG5Px3Bf19zIVvyT3D1X5Zy38w17D9aEu7wTIhVVysfbyxgVM/2dmsJ0yRYIghCdGQEE4dksPgXl3PXqG7M3ljA6CcX84cPN3P0VEW4wzMhsmbfUQ4WlzPuotom7jPGeywRuKB1XDS/GNeLRQ+M4tp+5zH137sY8fhCnluYY1coN0GzNxQQHSmMvqB9uEMxxhWWCFzUsU0LHr/xEj6+bwSDurbjiTnbGPnEQl5fsYdKG2HUJKgqH28qYFj3ZJuj2DQZQSUCEUkSkXkissP5mVhLmctFZG3Ao0xErnW2TRORXQHb+gYTT2PRMzWBlyZm8faPBtM5KZ6H3t3IFf+7mFnrDlBtt7z2tE0HTrD/aKmdFjJNSrAtggeB+araA5jvPP8SVV2oqn1VtS8wGigB5gYUeeD0dlVdG2Q8jcqAjCTe/tFgpk7KIi46kp/MWMM3//wp87cU2pBTj5q9MZ/ICOGK3pYITNMRbCIYD0x3lqcD156l/I3AbFVtNkNrRITRvVL56CfDeWZCX0orq7hjejbXP7+Mf+ccCnd4pp4+3ljAoMwkklrGhDsUY1wTbCJIVdV8Z7kASD1L+QnAjBrrHhGR9SLylIjE1rWjiEwWkWwRyS4qKgoi5PCIiBDG9+3EJz8byR+v60PB8TK+89IKJkz5zCbF8YjNB06ws+iUnRYyTY6c7RSFiHwC1FbzHwKmq2rbgLJHVfUr/QTOto7AeuA8Va0MWFcAxABTgJ2q+vDZgs7KytLs7OyzFWvUyiqrmLlyL88u3Mmhk+UM75HMfWN7cmmXWg+faQR+MmMN87cUsvSXo0m0FoHxIBFZrapZNdef9WoYVR17hhctFJGOqprvfKgfPMNL3Qy8ezoJOK99ujVRLiKvAPefLZ6mIi46kklDM7l5QDqvLd/DC4tzueH5ZYzsmcK9Y3vQv7MlhMZk16FTfLj+AD8Y3tWSgGlygj01NAuY6CxPBN4/Q9lbqHFayEkeiIjg71/YGGQ8nhMfE8XkEd349BeX84tx57N+/zGu/+syvj91Jav32DwIjcXzi3KIjozgjuGZ4Q7FGNcFmwgeBa4QkR3AWOc5IpIlIi+dLiQiGUA6sLjG/q+LyAZgA5AM/CHIeDyrZWwUd43qztJfjubBq3qxMe84Nzy/jO+8tJzluYfDHV6ztv9oCf/8PI9bBnamfUJcuMMxxnVn7SNojJpCH8HZlFT4eGPFXl5YnMuhk+UMyEjkrsu7M6pnCv4GlGkov3lvIzNX7WXxA5dzXtsW4Q7HmK+trj4Cu7K4kYqPieLO4V1Z+svL+f23epN3tJTbXlnF1X9ZyofrD1BlF6Y1iMITZbyZvY8b+qdZEjBNliWCRu50p/KiBy7n8RsvprSiinveWMOYP/knxymrtHsZhYqq8pv3/N1Wd43qHuZojAkdSwQeERMVwc1Z6cz72Uie/05/2rSI5qF3NzLssYU8u2AHx0rsbqdum7XuAHM3F/LzK3rSuV18uMMxJmSsj8CjVJXPcg8zZUkui7YV0SI6kpuz0rh9WCZd2rUMd3ied7C4jCufWkJmckve+dEQIiOsX8Z439e+jsA0TiLCkG7JDOmWzLaCYv72aS5vrNzL35fvYewFqdw+NJPLuiZZx/LXoKr8+t2NlFRU8cSNl1gSME2eJYIm4PwOCTx50yX84hvn8+ryPby2fA/zNhfSq0MCk4ZkMNEQuJsAAApYSURBVL5vJ1rERIY7TM94fcVe5m4u5FdX9aJ7+1bhDseYkLNTQ01QWWUV763JY9qy3WwtKKZtfDQ3Z6XznUGd7bTRWczfUsgP/p7N8B4pTJ00wFoDpkmp69SQJYImTFVZuesI0z/bzdxNhfiqlRE9U/jOoM6M7tWe6EgbKxBo3b5jTJiynG7tW/Lm5MG0jLUGs2larI+gGRIRBnVtx6Cu7Sg8UcbMlfuYsXIvP3x1Ne0TYrkpK42bs9KtlQDkFp3k9mmraNcqhqmTBlgSMM2KtQiaGV9VNYu2FTFj5V4WbjtItcKgzCRuzkrnqj4diI9pfh+Aq/cc4c7p2YgIb/9oMN1SrF/ANE12ash8Rf7xUv75eR5vZe9jz+ES4mMiGXdRB67vl8bgbu2axfnxjzbkc9+bazmvTRzTbhtIRrK1jkzTZYnA1Ol0X8K7a/L41/p8ist9tE+I5T8u7sg1l5xH3/S2TW4YaoWvmmcX7OAvC3Po3zmRv30/y2YdM02eJQJzTsoqq/hkSyEfrDvAwq1FVFRV06ltC666qANX9elAv/REIjzeUth84AT3v72OzfknuKF/Go9cdxFx0Ta81jR9lghMvZ0oq2TOxgJmbyxg6Y5DVFRVk5IQy5he7Rl7QSpDuyd76vqEYyUVvLA4l5eX5tKmRTR/vK4PV15o006a5sMSgQnKibJKFmw5yLzNhSzeXsTJch+xUREM6tqOET2SGdkzhe7tWzXKU0jFZZVMXbqblz7N5WSFj+v6duLXV/e2U0Gm2bFEYFxT4atmee5hFm0rYvH2g+wsOgVAcqtYBndrx2Vdk8jqkkSP9q3CdhpJVVm77xgzV+7jg/UHKKmo4sreqfzsyp706tA6LDEZE24huY5ARG4Cfg9cAAxU1Vo/nUVkHPAMEAm8pKqnZzLLBGYC7YDVwPdU1W6j2cjFREUwomcKI3qmAL3Zd6SEZTsP8dnOwyzbeZgP1h0AICEuir7pbenTqQ19OrXhok5t6NS2RciSQ7mvilW7jrJg60EWbjvIrkOniI+J5FsXn8d3L+tCn7Q2IXlfY7wuqBaBiFwAVAMvAvfXlghEJBLYDlwB7AdWAbeo6mYReQv4p6rOFJEXgHWq+vzZ3tdaBI2XqrLncAmr9xwle89R1uw9Ss7Bk/iciXRaREfSvX0rurdvRXpSPOmJLUhLjCclIZaUVrG0bhF1xtNLqkpZZTUHjpey70gJ+46WsiX/BBv2H2drwQkqq5SYqAgGd23HuIs6cPXFHUmIi26oX9+YRi0kLQJV3eK8+JmKDQRyVDXXKTsTGC8iW4DRwK1Ouen4WxdnTQSm8RIRMpJbkpHckhsuTQP8I5G2FRSzOf8EOwpPsuNgMStyD/Pe2jxqfg+JihBaxkYRHxNJi5hIBFCFKlVOlfs4Ueqjoqr6S/skxEVxcVob7hjWlQEZiQzp5q1ObGPCrSEuI+0E7At4vh8YhP900DFV9QWs71TXi4jIZGAyQOfOnUMTqQmJuOhILklvyyXpbb+0vsJXzYFjpeQdK+XQyXKKiss5fKqCknIfJRVVlDizr0WIECHQMjaK1nHRtG4RRYfWcU6LIp72CbGeH9JqTDidNRGIyCdAbWPsHlLV990PqXaqOgWYAv5TQw31viZ0YqIivmg9GGPC56yJQFXHBvkeeUB6wPM0Z91hoK2IRDmtgtPrjTHGNKCGuA/xKqCHiGSKSAwwAZil/l7qhcCNTrmJQIO1MIwxxvgFlQhE5DoR2Q8MBv4lInOc9eeJyEcAzrf9e4A5wBbgLVXd5LzEL4GfiUgO/j6Dl4OJxxhjTP3ZBWXGGNNM1DV81KaoMsaYZs4SgTHGNHOWCIwxppmzRGCMMc2cJzuLRaQI2PM1d08GDrkYjlssrvqxuOrH4qqfphpXF1VNqbnSk4kgGCKSXVuvebhZXPVjcdWPxVU/zS0uOzVkjDHNnCUCY4xp5ppjIpgS7gDqYHHVj8VVPxZX/TSruJpdH4Exxpgva44tAmOMMQEsERhjTDPX5BOBiDwhIltFZL2IvCsibesoN05EtolIjog82ABx3SQim0SkWkTqHA4mIrtFZIOIrBWRkN9prx5xNfTxShKReSKyw/mZWEe5KudYrRWRWSGM54y/v4jEisibzvYVIpIRqljqGdckESkKOEZ3NlBcU0XkoIhsrGO7iMifnbjXi0j/RhDTKBE5HnCsfhvqmJz3TReRhSKy2flfvLeWMu4eL1Vt0g/gSiDKWX4MeKyWMpHATqArEAOsA3qHOK4LgPOBRUDWGcrtBpIb8HidNa4wHa/HgQed5Qdr+zs62042wDE66+8P3AW84CxPAN5sJHFNAp5tqPoU8L4jgP7Axjq2fxOYDQhwGbCiEcQ0CvgwDMeqI9DfWU4Attfyd3T1eDX5FoGqztX/mxd5Of6Z0GoaCOSoaq6qVgAzgfEhjmuLqm4L5Xt8HecYV4MfL+f1pzvL04FrQ/x+Z3Iuv39gvO8AY0Qk1BMrh+Pvck5UdQlw5AxFxgN/V7/l+Gcv7BjmmMJCVfNV9XNnuRj/PC4153N39Xg1+URQw+34s2hNnYB9Ac/389UDHy4KzBWR1SIyOdzBOMJxvFJVNd9ZLgBS6ygXJyLZIrJcREKVLM7l9/+ijPNF5Dj+yZdC6Vz/Ljc4pxPeEZH0WraHQ2P9HxwsIutEZLaIXNjQb+6cUuwHrKixydXjddY5i71ARD4BOtSy6SFVfd8p8xDgA15vTHGdg2Gqmici7YF5IrLV+SYT7rhcd6a4Ap+oqopIXeOeuzjHqyuwQEQ2qOpOt2P1sA+AGapaLiI/xN9qGR3mmBqrz/HXp5Mi8k3gPaBHQ725iLQC/gHcp6onQvleTSIRqOrYM20XkUnA1cAYdU6w1ZAHBH4zSnPWhTSuc3yNPOfnQRF5F3/zP6hE4EJcDX68RKRQRDqqar7TBD5Yx2ucPl65IrII/7cptxPBufz+p8vsF5EooA1w2OU46h2XqgbG8BL+vpfGICR1KhiBH76q+pGI/FVEklU15DejE5Fo/EngdVX9Zy1FXD1eTf7UkIiMA34BXKOqJXUUWwX0EJFMEYnB37kXshEn50pEWopIwull/B3ftY5waGDhOF6zgInO8kTgKy0XEUkUkVhnORkYCmwOQSzn8vsHxnsjsKCOLyENGleN88jX4D//3BjMAr7vjIa5DDgecCowLESkw+l+HREZiP/zMtTJHOc9Xwa2qOr/1lHM3ePV0D3iDf0AcvCfS1vrPE6P5DgP+Cig3Dfx987vxH+KJNRxXYf/vF45UAjMqRkX/tEf65zHpsYSV5iOVztgPrAD+ARIctZnAS85y0OADc7x2gDcEcJ4vvL7Aw/j/8IBEAe87dS/lUDXUB+jc4zrf5y6tA5YCPRqoLhmAPlApVO/7gB+BPzI2S7Ac07cGzjDSLoGjOmegGO1HBjSQMdqGP6+wfUBn1vfDOXxsltMGGNMM9fkTw0ZY4w5M0sExhjTzFkiMMaYZs4SgTHGNHOWCIwxppmzRGCMMc2cJQJjjGnm/j+wSaPC1yIHwQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "\n", "def f(x):\n", " y = 0.\n", " for i in range(100):\n", " y = np.sin(y + x)\n", " return y\n", "\n", "x = np.linspace(-2, 2, 100)\n", "plt.plot(x, [f(x_) for x_ in x])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mathematically, it's something like:\n", "\n", "$$f(x) = \\sin(x + \\sin(x + \\sin(x + ...\\sin(x)))...)$$\n", "\n", "Good luck differentiating that and getting a nice closed form. If this is not complicated enough for you, feel free to add some `if` statements. \n", "\n", "We can use `autograd`, an automatical diff package, to compute _exact, pointwise_ derivatives." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.26242653107144726" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from autograd import grad\n", "from autograd import numpy as np\n", "\n", "def f(x):\n", " y = 0.\n", " for i in range(100):\n", " y = np.sin(y + x)\n", " return y\n", "\n", "\n", "grad(f)(1.)\n", "# magic! " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, you can string together these pointwise derivatives into a function:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3Bc5Zkm8Oftu64tWWpLsixjG5uLMdgwwtxhwiUFhODNJjskmbAhhDFkk2wyla1MZrOzNUnVzO5saphkJpNMGGCADSGZcEkMFRIgmIAXbJDBsjHGxpaNLevWlq27+v7uH31aluRuqSX16T7n+PlVqXSkPup+1ZIeff197zlHVBVERGRdrlIXQEREM2NQExFZHIOaiMjiGNRERBbHoCYisjiPGXdaX1+vy5cvN+OuiYgcaceOHcdVNZTtNlOCevny5WhrazPjromIHElEPsx1G6c+iIgsjkFNRGRxDGoiIotjUBMRWRyDmojI4hjUREQWx6AmIrI4U/qoiaxkNJrAI68fRjSeBACU+z2468rlCHjdJa6MKD8ManK8rQeO43u/2zflcxc2B3HVqvoSVUQ0N5z6IMeLGCPp33/jOvz6y1dN+RyRHTCoyfGiiRQAwOd2wedJ/8rHjM8R2QGnPsjxMkHt97gQS7qmfI7IDmYdUYvIuSKyc9LbkIh8vRjFERVCbCKo3fBzRE02NOuIWlX3AVgPACLiBnAMwDMm10VUMJlQ9nlc8GVG1EkGNdnHXOeobwBwUFVzno6PyGomB7Xf7Z7yOSI7mGtQfxrAE9luEJFNItImIm3hcHjhlREVSDSRhNslcLtkYjExmmDXB9lH3kEtIj4AtwP4ZbbbVfUBVW1V1dZQKOtFCohKIpZITcxNs+uD7GguI+pbALytqr1mFUNkhlgyNRHQbpfA4xIGNdnKXIL6M8gx7UFkZbFECj73qV91n8fFoCZbySuoRaQCwE0Anja3HKLCiyZOjaiBdFCzj5rsJK8DXlR1FECdybUQmWLyHDVgHPjCoCYb4SHk5HjpEfWpM+X5PC7E2EdNNsKgJsebvJgIpM/5wRE12QmDmhwvGk/CP2Ux0c0+arIVBjU5XiyZgt87dY6ai4lkJwxqcjy255HdMajJ8WLT2vP8XEwkm2FQk+Od1kftdiEaZ1CTfTCoyfFO66P2ckRN9sKgJsdjex7ZHYOaHC+9mDjtgBcGNdkIg5ocL5pIZjnXB/uoyT4Y1ORoqZQintRp5/pwc0RNtsKgJkfLLBpOH1FzMZHshEFNjpYJZP+0xcR4UpFKaanKIpoTBjU5WqZfevqIGgBH1WQbDGpytGwjav/EBW4Z1GQPDGpytMyi4fRDyCffRmR1+V6Kq0ZEnhSR90Vkr4hcYXZhRIWQacOb3kc9+TYiq8vrUlwAfgDgt6r6KRHxASg3sSaigsk2ovZxRE02M2tQi0gQwLUA7gIAVY0BiJlbFlFhZMJ4eh81wMVEso98pj5WAAgD+DcReUdEHjSuSj6FiGwSkTYRaQuHwwUvlGg+so6o3RxRk73kE9QeAJcA+LGqXgxgFMC3pu+kqg+oaquqtoZCoQKXSTQ/0RmmPtj1QXaRT1B3AuhU1e3Gx08iHdxEljcR1G7OUZN9zRrUqtoD4KiInGt86gYA75laFVGBZOahA16255F95dv18VUAjxsdHx0AvmBeSUSFMzFHnbU9j0FN9pBXUKvqTgCtJtdCVHATfdRZj0xkHzXZA49MJEfL3vXhnnIbkdUxqMnRsvZRe3lSJrIXBjU5GvuoyQkY1ORo0UQKIoDHJROf42Ii2Q2DmhwtlkzB53ZB5PSg5oia7IJBTY4WS6SmzE8D6dG1SxjUZB8ManK0aCIFn8c95XMiwusmkq0wqMnRoonkaSNqIL2gGI2zj5rsgUFNjhZLpKZ0fGT4PG6OqMk2GNTkaNnmqIF0XzW7PsguGNTkaLFk9hG13+PiYiLZBoOaHC0aT005xWmGjyNqshEGNTlarhG1jyNqshEGNTnaTHPUDGqyCwY1OVrurg/2UZN9MKjJ0aKJ5GkHvABGHzXPR002kdeFA0TkMIBhAEkACVXlRQTIFmKJ3IuJnPogu8j3UlwA8BFVPW5aJUQmiCVTE+efnszvcTOoyTY49UGOFuWImhwg36BWAC+IyA4R2ZRtBxHZJCJtItIWDocLVyHRAkRzdH2wj5rsJN+gvlpVLwFwC4Avi8i103dQ1QdUtVVVW0OhUEGLJJoPVc3d9eHmiJrsI6+gVtVjxvs+AM8A2GBmUUSFEE8qAGTvo/a6EGV7HtnErEEtIhUiUpXZBvBRAO+aXRjRQmX6pLOe68MYUatqscsimrN8uj4aADxjXMrIA+BnqvpbU6siKoDM+aZzLSYCRldIlj5rIiuZNahVtQPAuiLUQlRQp0bUWQ54mXTdRAY1WR3b88ixMouF2c/14Z6yD5GVMajJsTIhnOtcHwB4vg+yBQY1OVZ0pqA25q2jcQY1WR+DmhxrxqDmiJpshEFNjjXzHLVryj5EVsagJsfKjJZzHUIOgIeRky0wqMmxTvVR527P4zmpyQ4Y1ORYEyPqHKc5BTj1QfbAoCbHmmjPy3JkIueoyU4Y1ORY7KMmp2BQk2Oxj5qcgkFNjjVje56XI2qyDwY1OdZMpznNjKg5R012wKAmx4rOsJjo42Ii2QiDmhwrmkjC53bBOJf6FOyjJjthUJNjxXJc2Bbg1AfZC4OaHCvXhW0BQETSVyLnYiLZQN5BLSJuEXlHRJ4zsyCiQpkpqIFT100ksrq5jKi/BmCvWYUQFVp0lqD2eVw8KRPZQl5BLSJLAXwMwIPmlkNUODPNUQPp/mqOqMkO8h1Rfx/ANwHk/K0WkU0i0iYibeFwuCDFES1ELDn7iJpBTXYwa1CLyG0A+lR1x0z7qeoDqtqqqq2hUKhgBRLNVyyRytpDncGgJrvIZ0R9FYDbReQwgJ8DuF5EfmpqVUQFEE0k85ijZh81Wd+sQa2qf6mqS1V1OYBPA3hZVT9nemVEC5Seoz79ogEZfo+b5/ogW2AfNTnWrF0fbM8jm/DMZWdVfQXAK6ZUQlRgs/VR+zwujI4lilgR0fxwRE2OFU2k4OdiIjkAg5ocK5ZMZb1eYgb7qMkuGNTkWPm05/HIRLIDBjU51mzteX4GNdkEg5oca9bFRLcLMfZRkw0wqMmREskUUoqZ+6i97KMme2BQkyPNdL3EjEwftaoWqyyieWFQkyNF47mvl5jh87iQUiCRYlCTtTGoyZHyGlHzArdkEwxqcqRM+M52PurJ+xJZFYOaHCnTdpfXiJoLimRxDGpypMzpS2caUWfmrzPz2URWxaAmR4rNaUTNXmqyNgY1OdKpOeqZz0cNgEcnkuUxqMmR8un64GIi2QWDmhwp3z5qgCNqsr58Lm4bEJE3RaRdRPaIyHeKURjRQrCPmpwknyu8RAFcr6ojIuIFsFVEnlfVbSbXRjRv7KMmJ5k1qDV9IoQR40Ov8cZjbsnSMu15M89RczGR7CGvOWoRcYvITgB9AF5U1e3mlkW0MMdHYgCAugp/zn1qK7wAgP7RaFFqIpqvvIJaVZOquh7AUgAbRGTt9H1EZJOItIlIWzgcLnSdRHPSMxhBdcCDMl/u9rz6Cj88LkH3YKSIlRHN3Zy6PlR1AMAWADdnue0BVW1V1dZQKFSo+ojmpWcogsZgYMZ9XC7B4io/ehnUZHH5dH2ERKTG2C4DcBOA980ujGgheociaKieOagBoCEYQM8Qg5qsLZ8RdROALSKyC8BbSM9RP2duWUQL0zMYQdMsI2oAaGJQkw3k0/WxC8DFRaiFqCASyRSOj0TRmM+IujqAP+zjmgpZG49MJMcJj0SR0vS0xmwaqwMYjSUxHIkXoTKi+WFQk+NkujjyGVFnFhx7uKBIFsagJsfJdHHktZho7MN5arIyBjU5TiZ0Z2vPA06NujmiJitjUJPj9AxF4HO7sKjcN+u+mTDv5YiaLIxBTY7TOxjB4mo/XC6Zdd+A142aci+nPsjSGNTkOD1DkbwWEjMaqwPoGeT5Psi6GNTkOD2Dkbxa8zIaqgPoGRo3sSKihWFQk6OoKkfU5DgManKUofEEIvFUXoePZzQGA+gfjSKe5HmpyZoY1OQomUXBfHqoMxqDAagCfcMcVZM1MajJUebSQ53BXmqyOgY1OUrPYHpRcC5z1A0MarI4BjU5SmZRcHF17ktwTTdxvg/2UpNFMajJUXqGIqir8E1cuDYfteVe+DwuHp1IlsWgJkfJ98ouk4mI0aLHoCZrYlCTo/QMzn6txGwaq3mlF7KufK6Z2CIiW0TkPRHZIyJfK0ZhRPMxnxE1kL7IAKc+yKryGVEnAHxDVdcAuBzAl0VkjbllEc1dNJFE/2hsTh0fGY3VfnQPRqCqJlRGtDCzBrWqdqvq28b2MIC9AJrNLoxorvqG0h0fjcH8Oz4yGqoDiCVSGBjjJbnIeuY0Ry0iy5G+0O32LLdtEpE2EWkLh3mxUCq+Uwe7lM35a5uMr+E8NVlR3kEtIpUAngLwdVUdmn67qj6gqq2q2hoKhQpZI1FeeuZwrcTpMqNwBjVZUV5BLSJepEP6cVV92tySiOanex5HJWZkRuFdAzzdKVlPPl0fAuAhAHtV9X7zSyKan93HhtBYHUCw3Dvnr22qDqAq4MG7x057sUhUcvmMqK8CcCeA60Vkp/F2q8l1Ec1Z+9EBrG+pmdfXulyCdUtr0H50oMBVES2cZ7YdVHUrgNkvPkdUQv0jURw5MYbPXrZs3vexvqUGP/7DQYzHkijz5X8IOpHZeGQiOUJ7Z3okPN8RdeZrkynF7mODhSqLqCAY1OQIO48MwCXAhc3Bed/H+mXpkN959GShyiIqCAY1OcI7RwdwTkMVKvyzzublVF/pR8uiMuzkPDVZDIOabC+VUrQfHcDFy+Y/7ZGxvqUWO48wqMlaGNRke4f6RzEUSSxofjpjfUsNugYj6OOBL2QhDGqyvcwIeH1L7YLvKxP273D6gyyEQU22t/PoACp8bqxaXLng+7pgSTW8buE8NVkKg5psb+fRAVy0tAZu18Lb/QNeN85vquY8NVkKg5psLRJPYm/30ERrXSGsb6nBrs4BJFM8NzVZA4OabG1P1yASKS3IQmLG+pYajMaSONA3UrD7JFoIBjXZ2jvGFMXFBQ5qgAe+kHUwqMnWfr+3DytDFVg8j1Ob5rKivgIN1X68tLevYPdJtBAMarKtnsEIth3qx+3rlhT0fkUEH79oCV7Z14dBXpqLLIBBTbb13K4uqKLgQQ0AG9c3I55UPP9ud8Hvm2iuGNRkW5vbu3BhcxArQwvvn55ubXM1VtRXYHN7V8Hvm2iuGNRkSx3hEezqHMTG9YUfTQPp6Y/b1y3BGx396OXh5FRiDGqypc3tXRABbrvInKAGgNvXL4Eq8CxH1VRi+Vwz8WER6RORd4tRENFsVBWb27tw2YpFaAwWrttjurNDlVjbXM2gppLLZ0T9CICbTa6DKG97uobQER7FxvXNpj/WxnXNaO8cxKHjo6Y/FlEuswa1qr4K4EQRaiHKy5M7OuF1C25Z22j6Y922rgkiwNNvd5r+WES5FGyOWkQ2iUibiLSFw+FC3S3RFH3DEfz8rSP4+EVLUFPuM/3xmoJluOn8Bjzy+mH2VFPJFCyoVfUBVW1V1dZQKFSouyWa4kdbDiKeVPzXG1YX7TH//KZzMBxJ4MGtHUV7TKLJ2PVBttE9OI6fbT+CT12yFMvrK4r2uOc3VeNjFzbh4a2HcGI0VrTHJcpgUJNt/PDlA1AovnrDqqI/9tdvXI2xeBI/efVg0R+bKJ/2vCcAvAHgXBHpFJEvml8W0VRHT4zh39uO4o5LW7C0trzoj7+6oQob1y3BY69/iPBwtOiPT2e2fLo+PqOqTarqVdWlqvpQMQojmuzvX9gHEcFXPlK8uenpvnbjOYglU/jB7/eXrAY6M3Hqgyzvt+9241c7u3DvtStNPcBlNivqK/D5K5bjp9uO4LUP2NlExcOgJkvrG47gL5/ejQubg0Xt9Mjlmzefi1WLK/HfftmOgTEuLFJxMKjJslQVf/HkLozFkviHO9bB6y79r2vA68b371iP/pEY/urXe0pdDp0hSv+bT5TDT7cfwZZ9Yfz3W8/HqsVVpS5nwtrmIP78pnPwbHsXfvXOsVKXQ2cABjVZ0qv7w/jO5j247pwQ7rz8rFKXc5p7r12JS5fX4i+e2oW2wzzDApmLQU2Ws6tzAPf9dAdWN1Thnz57MVwuKXVJp/G4XfjJna1orinD3Y+8hf29w6UuiRyMQU2Wcuj4KL7wb29hUYUPj37hUlQHvKUuKadFFT48evcGBLxufP7hN9E1MF7qksihGNRkGR/0DuNP/3UbFMBjd28o6JXFzdKyqByP3r0BI5EEPvuv2/BhP0+HSoXHoCZL2N7Rj0/++HXEU4rH7t5gynUQzXJ+UzUeuXsDBsbj+I8/eh3tRwdKXRI5DIOaSu7Z9i7c+dCbCFX58fSXrsTa5mCpS5qzPzqrFk996UqU+9349APb8NJ7vaUuiRyEQU0lE4kn8Ve/ehdffeIdXLQ0iKe+dCVaFhX/PB6FcnaoEk996UqsWlyJex5rw9/+Zi9iiVSpyyIHYFBTSezvHcbGH/4//N9tH+LPrlmBx//ssqJcCMBsi6sC+Pd7r8CfXrYMD7zagU/++HVexosWjEFNRRWJJ3H/C/tw2z9uRf9oFI984VJ8+2Nr4Pe4S11awZT53PibT1yIf/ncH+HIiTHc8oNX8c9bDnB0TfPmKXUBdGZQVbz8fh/++tk9OHpiHBvXL8H/+NgahKr8pS7NNDevbcS6liC+s/k9fO93+/DU25347u1rcfXq+lKXRjYjqlrwO21tbdW2traC3y/Z0/aOfvz9C/vx5uETWLW4Et/deAGuPPvMCqst+/rw15v34MP+MVy1qg7f+Oi5uGRZbanLIgsRkR2q2pr1NgY1mUFVsfXAcfzkDx3YeuA4Flf58dXrV+GOS5fB5zkzZ9wi8SQe334EP9pyAP2jMXzk3BDuve5sXLZiEUSsd/QlFdeCg1pEbgbwAwBuAA+q6v+eaX8G9ZlrPJbEs+1deGjrIezrHUZ9pR/3XbcSn7v8LAS8zpmHXojRaAKPvnEYD76Wvgbj2uZqfPHqFbhlbROfozPYgoJaRNwA9gO4CUAngLcAfEZV38v1NQzqM4uqYvexQfziraPYvLMLw9EEzmuswj3XrMTH1zU5aqGwkCLxJJ555xgefK0DB8OjCJZ58YmLm/EnrS1Ys6S61OVRkc0U1PksJm4AcEBVO4w7+zmAjQByBjWdGT7oHcazu7rxXHsXOo6PIuB14da1TfiTS1v4cj4PAa8bn9mwDHe0tuD1g/34RdtR/Gz7ETzy+mGc01CJ2y5agtsuarLVUZpkjnyCuhnA0UkfdwK4bPpOIrIJwCYAWLZsWUGKI2tJJFN4+8gAXtrbi5fe60XH8VGIAFesrMM916zEbeuaLH0SJatyuQRXr67H1avrMTAWw7PtXXi2vRv3v7gf97+4H6sXV+LGNQ248fwGrG+pgduCZxMkc+Uz9fEpADer6j3Gx3cCuExVv5Lrazj14QyqioPhUbzR0Y/X9ofxxsF+DEcT8LoFl6+sw01rGnDzBY22OHmSHXUPjuP53T14aW8vth86gWRKESzz4qpVdbh6VQhXnF2H5XXlfOXiEAud+jgGoGXSx0uNz5HDxBIpvNc9hHeOnETb4ZPYfugEjo9EAQDNNWW4bV0TrlkdwjWr61HFkbPpmoJluPvqFbj76hUYHIvjDx+E8dr+MLYeOI7f7O4BADRU+7FhRR0uXV6Li1tqcV5TlSUuWUaFlc+I2oP0YuINSAf0WwA+q6o5LxjHEbX1xZMpHOgbwbvHBvHusUHsOjaIPV1DE0fPLQkGcNnKOmxYsQiXr+TIzUoyr3S2dfTjzUMnsP1QP3qH0v9QA14X1i4JYm1zEBctTb9fWV8BD8Pb8grRnncrgO8j3Z73sKr+zUz7M6itI5VSHBsYxwd9w9jfO4L9vcN4v3sYB/pGEEumQ7nC58YFS9J/2JecVYuLl9WgKVhW4sopX6rpn/E7Rwbw9pGT2N2Z/qc7Hk8CAHweF85pqMR5jdU4p6ESqxuqsHpxJZYEyyx59ZwzFQ94cThVRXgkiiP9YzjcP4bDx0dxqH8UHeFRdIRHEJ10jomGaj/ObazG+U1VWNNUjQuWVGNFfSUXqBwmmVIcDI9gT9cg9nYP472uIezrHUZ4ODqxT5nXjZWhCqwMVWJFXTmW11fgrLoKLFtUjvpKH19BFRmD2uYSyRT6hqPoHhzHsYEIugfG0XlyHJ0nx9B5chxHT44hEj8Vxm6XoKW2DCtDlTjb+ENcvbgSqxdXIVjOueUz2cnRGPb3DuNgeBQHwyM4GB5BR3gUnSfHkJoUBWVeN1oWlaGlthxLa8vQXFuGJTVlaAqWYUlNAKFKP6dTCmyhi4lkkmgiif6RGI6PRBEePvXWOxxB71AUfUMR9AxFEB6OTvkjAoBgmRdLa8uwor4C150TwrK6crQsKsfyugosrS3jghJlVVvhw2Ur63DZyropn48lUjh6cgwf9o/iSP8YjpwYx5ETYzg2MI43D53AcDQxZX+XpE/p2hAMoLHan9423oeq/AhV+VFX6UNdhf+MPWVAITGoCySRTGFwPI7B8ThOjsUxMBbDwFgcJ8diODEaw8mxGPpH0tsnRtPhPBRJZL2v2nIvFlcFsLjaj3Mbq9BYHUBjsAxNNQE015ShKRhg1wUVlM/jwtmhSpyd4+CawfE4ugfH0T0QwbGBcfQORdA9GEHPYASHjo9i+6ETGBiLZ/3aYJnXCG0fFlX4sKjCj0UVXtSW+9JvFV7UGNvBMi+qAx6O1qdhUCMdsqOxJEajCYxk3iLp98OROIYjCQxFJm2PxzEUiWNwPL09OB7HSDR76AKAxyWorfChttyLugo/1iypRl2FD/WVftRX+VFf6Z8YhdRX+njINVlOsMyLYJkX5zXmPrQ9Ek9OvDrsG45OvFo8PhJF/2gM/SNRdIRHsePDAZwciyE5/WXiJFV+D6qNxwyWeVFd5kF1wIuqgBdVAQ+qAumPK43tSn/6rcLvQWXAgwqfx1HrLpYPalVFIqWIxJOIxFOIxJOIJk5tR+IpjMeTGI8nEYklJ7bHYkmMxxLG+yRGJ22PRNPbo9EERmOJKfO7M6n0exAsO/WL0lwTwPlNVQiWeVFT5kOwzIOach+C5ZnRQnqkUB3wcGGGHC/gdWNpbTmW1s5+OTVVxdB4AifH0q82B8biGBg33o+lBz9DkfjEQOjD/jFjgJSYcVA0tR4XKv0elPvSAV7hc6PM50aFz4NyYzv93vjYm34LGNvlPjcCXhf8HjcC3vR2wOuG35N+X8zpRUsF9W3/9BqGIwlE4ylEE0lEE+kwnuEf74wyT/bED8ef/ri2ohwVPjfKM/+FfR5U+N3p/8qB9A+1yu9B1eT/2D4PW5mICkREECz3IljuxXJUzOlrkymdeLWbefU7HJn6ang0ljBeIacHbJlXzMORBHqHIhiNJhHJDOiMNsa5crtkIrT9Hhf8Hlf6Umz3XTGv+5uJpYJ69eIqqCr8Hjd8xjc++T9Y5r+b3+tK/+cz3tL/+dK3l/s8xm0ujmKJHMjtkokpkUJIpRTRRApjxqvuaCKJ8Vj6lXok82o9nkQ0nkIkMXU78z6WSCGaSKHMpNPUWiqo/+GO9aUugYjOMC6XoMx45V03++4lwaVVIiKLY1ATEVkcg5qIyOIY1EREFsegJiKyOAY1EZHFMaiJiCyOQU1EZHGmnI9aRMIAPpznl9cDOF7AcgqFdc0N65ob1jU3TqzrLFUNZbvBlKBeCBFpy3Xy7FJiXXPDuuaGdc3NmVYXpz6IiCyOQU1EZHFWDOoHSl1ADqxrbljX3LCuuTmj6rLcHDUREU1lxRE1ERFNwqAmIrK4kge1iHxPRN4XkV0i8oyI1OTY72YR2SciB0TkW0Wo6z+JyB4RSYlIznYbETksIrtFZKeItFmormI/X4tE5EUR+cB4X5tjv6TxXO0Ukc0m1jPj9y8ifhH5hXH7dhFZblYtc6zrLhEJT3qO7ilCTQ+LSJ+IvJvjdhGRfzRq3iUil5hdU551/bGIDE56rv5nkepqEZEtIvKe8bf4tSz7FPY5U9WSvgH4KACPsf13AP4uyz5uAAcBrATgA9AOYI3JdZ0P4FwArwBonWG/wwDqi/h8zVpXiZ6v/wPgW8b2t7L9HI3bRorwHM36/QP4LwD+xdj+NIBfWKSuuwD8sFi/T8ZjXgvgEgDv5rj9VgDPAxAAlwPYbpG6/hjAc8V8rozHbQJwibFdBWB/lp9jQZ+zko+oVfUFVc1cVngbgKVZdtsA4ICqdqhqDMDPAWw0ua69qrrPzMeYjzzrKvrzZdz/o8b2owD+g8mPN5N8vv/J9T4J4AYx/yKbpfi5zEpVXwVwYoZdNgJ4TNO2AagRkSYL1FUSqtqtqm8b28MA9gJonrZbQZ+zkgf1NHcj/V9oumYARyd93InTn5hSUQAviMgOEdlU6mIMpXi+GlS129juAdCQY7+AiLSJyDYRMSvM8/n+J/YxBgqDgOmXzMv35/JJ4+XykyLSYnJN+bDy398VItIuIs+LyAXFfnBjyuxiANun3VTQ56woF7cVkZcANGa56duq+mtjn28DSAB4vBg15VtXHq5W1WMishjAiyLyvjESKHVdBTdTXZM/UFUVkVx9n2cZz9dKAC+LyG5VPVjoWm3sWQBPqGpURO5FetR/fYlrsqq3kf59GhGRWwH8CsDqYj24iFQCeArA11V1yMzHKkpQq+qNM90uIncBuA3ADWpM8ExzDMDkkcVS43Om1pXnfRwz3veJyDNIv7xdUFAXoK6iP18i0isiTarabbzE68txH5nnq0NEXkF6NFLooM7n+8/s0ykiHgBBAP0FrmPOdanq5BoeRHruv9RM+X1aqMnhqKq/EZEfiUi9qpp+siYR8SId0o+r6tNZdinoc1byqQ8RuRnANwHcrqpjOXZ7C8BqEVkhIj6kF39M6xjIl4hUiEhVZhvphdGsK9RFVornazOAzxvbnx27fEwAAAE7SURBVAdw2shfRGpFxG9s1wO4CsB7JtSSz/c/ud5PAXg5xyChqHVNm8e8Hen5z1LbDOA/G50MlwMYnDTNVTIi0phZVxCRDUjnmdn/bGE85kMA9qrq/Tl2K+xzVuwV0ywrqAeQnsvZabxlVuKXAPjNtFXU/UiPvr5dhLo+gfS8UhRAL4DfTa8L6dX7duNtj1XqKtHzVQfg9wA+APASgEXG51sBPGhsXwlgt/F87QbwRRPrOe37B/BdpAcEABAA8Evj9+9NACvNfo7yrOt/Gb9L7QC2ADivCDU9AaAbQNz43foigPsA3GfcLgD+2ah5N2bogipyXV+Z9FxtA3Blkeq6Gum1qV2TcutWM58zHkJORGRxJZ/6ICKimTGoiYgsjkFNRGRxDGoiIotjUBMRWRyDmojI4hjUREQW9/8BONFa/NgvuysAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "grad_f = grad(f)\n", "\n", "x = np.linspace(-2, 2, 100)\n", "plt.plot(x, [grad_f(x_) for x_ in x])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To use this in our optimization routines:\n", "\n", " - we can automatically compute gradients that the optimizer can use. \n", " - Hessians can be exactly calculated (we will do this later)\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "T = (np.random.exponential(size=1000)/1.5) ** 2.3\n", "E = np.random.binomial(1, 0.95, size=1000)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def cumulative_hazard(params, t):\n", " lambda_, rho_ = params\n", " return (t / lambda_) ** rho_\n", "\n", "def log_hazard(params, t):\n", " lambda_, rho_ = params\n", " return np.log(rho_) - np.log(lambda_) + (rho_ - 1) * (np.log(t) - np.log(lambda_))\n", "\n", "def log_likelihood(params, t, e):\n", " return np.sum(e * log_hazard(params, t)) - np.sum(cumulative_hazard(params, t))\n", "\n", "def negative_log_likelihood(params, t, e):\n", " return -log_likelihood(params, t, e)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So `grad(negative_log_likelihood)` will find the gradient of `negative_log_likelihood` with respect to the first parameter, `params`. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-199.44295882 2906.3619735 ]\n" ] } ], "source": [ "grad_negative_log_likelihood = grad(negative_log_likelihood)\n", "\n", "print(grad_negative_log_likelihood(np.array([1., 1.]), T, E))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " fun: 243.1558229286153\n", " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([ 0.00101051, -0.00293405])\n", " message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", " nfev: 11\n", " nit: 7\n", " status: 0\n", " success: True\n", " x: array([0.4693 , 0.4273751])\n" ] } ], "source": [ "from scipy.optimize import minimize\n", "\n", "results = minimize(negative_log_likelihood, \n", " x0 = np.array([1.0, 1.0]),\n", " method=None, \n", " args=(T, E),\n", " jac=grad_negative_log_likelihood,\n", " bounds=((0.00001, None), (0.00001, None)))\n", "\n", "print(results)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from autograd import value_and_grad\n", "\n", "results = minimize(\n", " value_and_grad(negative_log_likelihood), \n", " x0 = np.array([1.0, 1.0]),\n", " method=None, \n", " args=(T, E),\n", " jac=True, # notice this set to True now.\n", " bounds=((0.00001, None), (0.00001, None)))\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " fun: 243.1558229286153\n", " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([ 0.00101051, -0.00293405])\n", " message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", " nfev: 11\n", " nit: 7\n", " status: 0\n", " success: True\n", " x: array([0.4693 , 0.4273751])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's continue this analytical-train 🚂 to Part 4! " ] } ], "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": 2 }