{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Writing a LogPDF\n", "\n", "Probability density functions in Pints can be defined via [models](writing-a-model.ipynb) and [problems](sampling-first-example.ipynb), but they can also be defined directly.\n", "\n", "In this example, we implement the [Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function) and run an [optimisation](optimisation-first-example.ipynb) using it." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The rosenbrock function is a two dimensional defined as\n", "\n", "```\n", "f(x,y) = -((a - x)^2 + b(y - x^2)^2)\n", "```\n", "\n", "where ``a`` and ``b`` are constants and ``x`` and ``y`` are variable. In analogy with typical Pints models ``x`` and ``y`` are our _parameters_." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, take a look at the [LogPDF](http://pints.readthedocs.io/en/latest/log_pdfs.html#pints.LogPDF) interface. It tells us two things:\n", "\n", "1. We need to add a method `n_parameters` that tells pints the dimension of the parameter space.\n", "2. Objects of our class should be _callable_. In Python, we can do this using the [special method](https://docs.python.org/3/reference/datamodel.html#object.__call__) `__call__`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The input to this method should be a vector, so we should rewrite it as\n", "\n", "```\n", "f(p) = -((a - p[0])^2 + b(p[1] - p[0]^2)^2)\n", "```\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of calling this method should be the logarithm of a normalised log likelihood. That means we should (1) take the logarithm of ``f`` instead of returning it directly, and (2) invert the method, so that it has a clearly defined _maximum_ that we can search for.\n", "\n", "So we should create an object that evaluates\n", "\n", "```\n", "-log(f(p))\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have all we need to implement a ```Rosenbrock``` class:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pints\n", "\n", "class Rosenbrock(pints.LogPDF):\n", " def __init__(self, a=1, b=100):\n", " self._a = a\n", " self._b = b\n", "\n", " def __call__(self, x):\n", " return - np.log((self._a - x[0])**2 + self._b * (x[1] - x[0]**2)**2)\n", "\n", " def n_parameters(self):\n", " return 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can test our class by creating an object and calling it with a few parameters:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-0.0\n", "-0.482426149244\n", "0.653926467407\n" ] } ], "source": [ "r = Rosenbrock()\n", "print(r([0, 0]))\n", "print(r([0.1, 0.1]))\n", "print(r([0.4, 0.2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wikipedia tells for that for ``a = 1`` and ``b = 100`` the minimum value should be at ``[1, 1]``. We can test this by inspecting its value at that point:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log\n", " # Remove the CWD from sys.path while we load stuff.\n" ] }, { "data": { "text/plain": [ "inf" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r([1, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get an error here, because the notebook doesn't like it, but it returns the correct value!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's try an optimisation:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Maximising LogPDF\n", "using Covariance Matrix Adaptation Evolution Strategy (CMA-ES)\n", "Running in sequential mode.\n", "Population size: 6\n", "Iter. Eval. Best Time m:s\n", "0 6 -7.122578 0:00.1\n", "1 12 -2.87918 0:00.1\n", "2 18 -0.755 0:00.1\n", "3 24 -0.755 0:00.1\n", "20 126 2.270349 0:00.1\n", "40 246 6.824817 0:00.1\n", "60 366 17.84721 0:00.1\n", "80 486 25.05412 0:00.2\n", "100 606 35.67203 0:00.2\n", "120 726 43.65702 0:00.2\n", "140 846 55.17039 0:00.3\n", "160 966 65.60115 0:00.3\n", "180 1086 72.08731 0:00.3\n", "200 1206 inf 0:00.3\n", "220 1326 inf 0:00.4\n", "240 1446 inf 0:00.4\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log\n", " # Remove the CWD from sys.path while we load stuff.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "260 1566 inf 0:00.4\n", "280 1686 inf 0:00.5\n", "300 1806 inf 0:00.5\n", "320 1926 inf 0:00.5\n", "340 2046 inf 0:00.5\n", "360 2166 inf 0:00.6\n", "380 2286 inf 0:00.6\n", "397 2382 inf 0:00.6\n", "Halting: No significant change for 200 iterations.\n" ] } ], "source": [ "# Define some boundaries\n", "boundaries = pints.RectangularBoundaries([-5, -5], [5, 5])\n", "\n", "# Pick an initial point\n", "x0 = [2, 2]\n", "\n", "# And run!\n", "xbest, fbest = pints.optimise(r, x0, boundaries=boundaries)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, print the returned point. If it worked, we should be at `[1, 1]`:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1. 1.]\n" ] } ], "source": [ "print(xbest)" ] } ], "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }