Add higher order derivatives

[JCGoran]

The DALI method (as described here) provides a straightforward generalization to Fisher matrices, and only requires minor tweaks to make it work:

• instead of having values in the constructor, have *values, i.e. let the user pass n x n, n x n x n, etc. tensors (or just supply the keywords matrix, flexion, quarxion to make it unambiguous)
• all math operators (with the exception of matrix multiplication) just act on all elements
• the Jacobian transformation requires a bit of care, but is probably straightforward

Now, what isn't super straightforward is the computation of the (un)marginalized constraints; since we are dealing with a non-Gaussian distribution, in particular (assuming zero mean):

$$ \mathcal{L}(\boldsymbol{\theta}) = \frac{1}{N} \mathrm{exp}(-\frac{1}{2} \boldsymbol{\theta}^T \mathsf{F} \boldsymbol{\theta} - H.O.) $$

we have two problems:

• we need to figure out the normalizing constant $N$; this can probably be accomplished by using an n-dimensional integration (for instance, scipy.integrate.nquad seems to be adequate for the job of up to, say, 7 parameters, otherwise some MC-based thing is preferred; this paper describes a nice, well-maintained alternative)
• we need to get the 2D contours somehow, since we can't just invert a Fisher matrix and call it a day. One idea is to integrate over n-2 dimensions, then define iso-probability contours (once we have the normalization, of course), and then find points in the x-y plane which correspond to those contours. Scipy seems to have some root-finding routines for multidimensional problems, but it's not immediately clear how to code everything up.