Probability density functions in Pints can be defined via models and problems, but they can also be defined directly.
In this example, we implement the Rosenbrock function and run an optimisation using it.
The rosenbrock function is a two dimensional defined as
f(x,y) = -((a - x)^2 + b(y - x^2)^2)
where a
and b
are constants and x
and y
are variable. In analogy with typical Pints models x
and y
are our parameters.
First, take a look at the LogPDF interface. It tells us two things:
n_parameters
that tells pints the dimension of the parameter space.__call__
.The input to this method should be a vector, so we should rewrite it as
f(p) = -((a - p[0])^2 + b(p[1] - p[0]^2)^2)
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.
So we should create an object that evaluates
-log(f(p))
We now have all we need to implement a Rosenbrock
class:
import numpy as np
import pints
class Rosenbrock(pints.LogPDF):
def __init__(self, a=1, b=100):
self._a = a
self._b = b
def __call__(self, x):
return - np.log((self._a - x[0])**2 + self._b * (x[1] - x[0]**2)**2)
def n_parameters(self):
return 2
We can test our class by creating an object and calling it with a few parameters:
r = Rosenbrock()
print(r([0, 0]))
print(r([0.1, 0.1]))
print(r([0.4, 0.2]))
-0.0 -0.482426149244 0.653926467407
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:
r([1, 1])
/usr/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log # Remove the CWD from sys.path while we load stuff.
inf
We get an error here, because the notebook doesn't like it, but it returns the correct value!
Now let's try an optimisation:
# Define some boundaries
boundaries = pints.RectangularBoundaries([-5, -5], [5, 5])
# Pick an initial point
x0 = [2, 2]
# And run!
xbest, fbest = pints.optimise(r, x0, boundaries=boundaries)
Maximising LogPDF using Covariance Matrix Adaptation Evolution Strategy (CMA-ES) Running in sequential mode. Population size: 6 Iter. Eval. Best Time m:s 0 6 -7.122578 0:00.1 1 12 -2.87918 0:00.1 2 18 -0.755 0:00.1 3 24 -0.755 0:00.1 20 126 2.270349 0:00.1 40 246 6.824817 0:00.1 60 366 17.84721 0:00.1 80 486 25.05412 0:00.2 100 606 35.67203 0:00.2 120 726 43.65702 0:00.2 140 846 55.17039 0:00.3 160 966 65.60115 0:00.3 180 1086 72.08731 0:00.3 200 1206 inf 0:00.3 220 1326 inf 0:00.4 240 1446 inf 0:00.4
/usr/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log # Remove the CWD from sys.path while we load stuff.
260 1566 inf 0:00.4 280 1686 inf 0:00.5 300 1806 inf 0:00.5 320 1926 inf 0:00.5 340 2046 inf 0:00.5 360 2166 inf 0:00.6 380 2286 inf 0:00.6 397 2382 inf 0:00.6 Halting: No significant change for 200 iterations.
Finally, print the returned point. If it worked, we should be at [1, 1]
:
print(xbest)
[ 1. 1.]