laplaceapprox.h File Reference

Laplace approximation of a multidimensional integral around a given point. More...

#include <cmath>
#include <cstddef>
#include <functional>
#include <vector>
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for laplaceapprox.h:

Go to the source code of this file.

Classes
struct	line::pfqn::LaplaceResult< T >
	Return value of laplaceapprox, mirroring [I, H, logI]. More...

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
LaplaceResult< T >	line::pfqn::laplaceapprox (const std::function< T(const std::vector< T > &)> &h, const std::vector< T > &x0)
	Laplace approximation of a multidimensional integral around a given point.

Detailed Description

Laplace approximation of a multidimensional integral around a given point.

Templated port of matlab/src/api/pfqn/laplaceapprox.m together with the two thirdparty helpers it calls, matlab/lib/thirdparty/num_hess.m and num_grad.m (central differences at step h, the Hessian built as a difference of gradients). Those helpers have no pfqn counterpart of their own, so they are ported here as num_grad / num_hess in the pfqn::detail namespace rather than given a header of their own.

I = h(x0) sqrt((2 pi)^d / det(-H)), H = Hessian of log h at x0.

MATLAB's retry ladder on a negative det(-H) – widen the differentiation step from 1e-5 to 1e-4 to 1e-3, then warn and carry on – is reproduced exactly, including the fact that the final value is still returned when the determinant stays negative. Callers (pfqn_nrl, pfqn_nrp) take the real part of the logarithm afterwards, which is what makes that path survivable.

logI IS the logarithm of I. laplaceapprox.m returns logI = log(h(x0)) + (d/2) log(2 pi) - (1/2) log(det(-H)), i.e. exactly log(I): the curvature term carries HALF the log-determinant, because I carries sqrt(det(-H)) in its denominator. That is worth stating explicitly because the two are easy to desynchronize – the reference itself carried a full log(detnH) at one point, which shifted every pfqn_nrl and pfqn_nrp value by (1/2) log det(-H), and pfqn_nrl / pfqn_nrp read logI, not I.

ARITHMETIC. Both the (2 pi)^{d/2} factor and the finite-difference Hessian are inexact, so the routine is gated on num_traits<T>::has_transcendental.

Definition in file laplaceapprox.h.