3

I'm trying to produce contour plots for parameters with physical limits using the Minuit2 minimizer which is a part of the ROOT data analysis framework. Unfortunately, Minuit2 seems intent on drifting the parameters into regions outside of their limits when I try to produce contour plots:

>>> from minuit2 import Minuit2
>>> def f(x,y):
...     if x < 0 or y < 0:
...             print 'x = %.2f, y = %.2f' % (x,y)
...             raise Exception
...     return x**2 + y**2
... 
>>> m = Minuit2(f)
>>> m.limits['x'] = 0, 10
>>> m.limits['y'] = 0, 10
>>> m.migrad()
>>> xy = m.contour('x','y',3)
Info in <Minuit2>: MnMinos UP value has changed, need to update FunctionMinimum class
x = -9.95, y = 0.00
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 4, in f
Exception

Has anybody else dealt with this or a similar problem? Are there any workarounds?

I've already asked this question on the ROOT forums, but I thought there might also be some stack overflow users who have dealt with this or a similar issue.

user545424
  • 15,713
  • 11
  • 56
  • 70

1 Answers1

0

Try your example without raising an exception

def f(x,y):
    return x ** 2 + y ** 2

and you will get reasonable xy contour points (i.e. within 1e-3 of the true contour).

Note that the parameter sigmas=3 in your contour call m.contour('x', 'y', 3) means that the contour for sigmas ** 2 == 9 will be computed and that contour points along the parameter limits are computed. As far as I can see this is not mentioned in the contour() pyminuit documentation). In your example the contour starts at (0, 0), goes up to (3, 0), along the circle to (0, 3), and back to (0, 0).

A common method is to implement parameter limits (arbitrary shapes, not only min / max) in your cost function by returning very high values for excluded parameters:

def f(x,y):
    if x < 0 or y < 0:
        return 1e10
    return x ** 2 + y ** 2

This does throw the optimizer out of the forbidden regions, but it does not prevent it to probe them sometimes (i.e. evaluate f there).

I don't know why contour() should strictly respect the limits you set via

m.limits['x'] = 0, 10
m.limits['y'] = 0, 10

Here's a short description of the contour algorithm used by Minuit (and Minuit2) and here is the documentation for the Minuit2 code in ROOT, I did not manage to find the actual C file showing the implementation.

Christoph
  • 2,790
  • 2
  • 18
  • 23
  • Simply not raising an exception is not possible. The parameter limits represent physical bounds and it is not possible to evaluate the chi^2 outside of these limits. I've tried simply returning a large number like you suggest, but Minuit dies in more complicated scenarios than the simple example I gave; for example, try `f(x,y,z)` and only setting limits on `x` and `y`. I don't understand why Minuit *shouldn't* respect the limits I've set. They are physical limits outside of which my function is undefined. – user545424 Jul 27 '12 at 17:08
  • 1
    When running `Minuit.migrad()` with parameter limits, Minuit does the parameter transformation described in Section 1.3.1 of the [Minuit user guide](http://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf). I assume `Minuit.contour()` is run without parameter transformation (to be you'd have to look at the code or find a description of the CONTOUR implementation) and you have no way to get the contour from Minuit in complicated cases. – Christoph Aug 01 '12 at 10:25
  • Maybe you can use some other package like e.g. [lmfit-py](http://newville.github.com/lmfit-py] or compute the contour yourself by evaluating your cost function on a grid of `(x,y)` points, or running `Minuit.migrad()` on a grid of `(x,y)` points fixing `x` and `y` if you have more parameters and then use e.g. matplotlib to compute contours? Of course this is slow, but since you do more of the actual work yourself, you have more control and can make it work in complicated cases. I could provide a code snippet if this suggestion sounds cryptic. – Christoph Aug 01 '12 at 10:29