Classes
struct	Radial
	log J(c) and the moments of the tilted law of the radius. More...
struct	Mode
	Mode, curvature and log-integrand at the mode of the simplex factor. More...

Functions
template<class T>
void	sym_eig (Matrix< T > A, std::vector< T > &d, Matrix< T > &V)
	Eigenvalues and eigenvectors of a symmetric matrix by cyclic Jacobi, ascending.
template<class T>
void	golub_welsch (const std::vector< T > &off, const T &mu0, std::vector< T > &x, std::vector< T > &w)
	Nodes and weights of the Gauss rule with the given Jacobi off-diagonal.
template<class T>
void	gauss_legendre (std::size_t n, std::vector< T > &x, std::vector< T > &w)
	N-point Gauss-Legendre rule on [-1,1].
template<class T>
void	gauss_hermite (std::size_t q, std::vector< T > &z, std::vector< T > &w)
	Q-point Gauss-Hermite rule of the probabilists' weight exp(-z^2/2).
template<class T>
T	radial_logf (const T &t, const std::vector< T > &c, const std::vector< T > &N, const std::vector< T > &Z, std::size_t M)
	Log-integrand of the radial integral in t = log v, Jacobian included.
template<class T>
Radial< T >	radial (const std::vector< T > &c, const std::vector< T > &N, const std::vector< T > &Z, std::size_t M)
	log J(c) = log int_0^inf exp(-v) v^(M-1) prod_r (Z_r + v c_r)^N_r dv, plus the moments of the tilted law of v that the simplex derivatives need: G_r = E[T_r], vbar = E[v] and Lam = cov(T) - diag(E[T^2]/N) = grad^2_c log J, with T_r(v) = N_r v/(Z_r + v c_r).
template<class T>
Mode< T >	simplex_mode (const Matrix< T > &L, const std::vector< T > &N, const std::vector< T > &Z)
	Mode and curvature of h(w) = log J(L'x(w)) + sum_i log x_i with J the exact radial integral.
template<class T>
std::vector< T >	softmax_gauge (const std::vector< T > &w)
	softmax of [w; 0], the logistic parametrisation of the simplex with gauge w_M = 0.
template<class T, class F>
T	aghq_rule (const F &h, const std::vector< T > &w0, const T &h0, const Matrix< T > &A, std::size_t q, std::size_t d)
	Log of the tensor Gauss-Hermite sum, accumulated with a running maximum; the det(A)^(-1/2) of the rule is applied by the caller.

Function Documentation

◆ aghq_rule()

template<class T, class F>

T line::pfqn::simplex::aghq_rule	(	const F &	h,
		const std::vector< T > &	w0,
		const T &	h0,
		const Matrix< T > &	A,
		std::size_t	q,
		std::size_t	d )

Log of the tensor Gauss-Hermite sum, accumulated with a running maximum; the det(A)^(-1/2) of the rule is applied by the caller.

A tensor rule is NOT invariant to the choice of A^(-1/2): the principal-axis frame is used, as in the reference results.

Definition at line 450 of file pfqn_simplex.h.

References aghq_rule(), gauss_hermite(), line::InputError::InputError(), and sym_eig().

Referenced by aghq_rule(), and line::pfqn::pfqn_aghq().

◆ gauss_hermite()

template<class T>

void line::pfqn::simplex::gauss_hermite	(	std::size_t	q,
		std::vector< T > &	z,
		std::vector< T > &	w )

Q-point Gauss-Hermite rule of the probabilists' weight exp(-z^2/2).

Definition at line 144 of file pfqn_simplex.h.

References gauss_hermite(), and golub_welsch().

Referenced by aghq_rule(), and gauss_hermite().

◆ gauss_legendre()

template<class T>

void line::pfqn::simplex::gauss_legendre	(	std::size_t	n,
		std::vector< T > &	x,
		std::vector< T > &	w )

N-point Gauss-Legendre rule on [-1,1].

Definition at line 131 of file pfqn_simplex.h.

References gauss_legendre(), and golub_welsch().

Referenced by gauss_legendre(), and radial().

◆ golub_welsch()

template<class T>

void line::pfqn::simplex::golub_welsch	(	const std::vector< T > &	off,
		const T &	mu0,
		std::vector< T > &	x,
		std::vector< T > &	w )

Nodes and weights of the Gauss rule with the given Jacobi off-diagonal.

Definition at line 112 of file pfqn_simplex.h.

References golub_welsch(), and sym_eig().

Referenced by gauss_hermite(), gauss_legendre(), and golub_welsch().

◆ radial()

template<class T>

Radial< T > line::pfqn::simplex::radial	(	const std::vector< T > &	c,
		const std::vector< T > &	N,
		const std::vector< T > &	Z,
		std::size_t	M )

log J(c) = log int_0^inf exp(-v) v^(M-1) prod_r (Z_r + v c_r)^N_r dv, plus the moments of the tilted law of v that the simplex derivatives need: G_r = E[T_r], vbar = E[v] and Lam = cov(T) - diag(E[T^2]/N) = grad^2_c log J, with T_r(v) = N_r v/(Z_r + v c_r).

Quadrature runs in t = log v, where the integrand is bounded at both ends, over two Gauss-Legendre panels meeting at the mode, each widened until the log-integrand has fallen 60 nats so the discarded tails are below 1e-26 in relative terms.

Definition at line 194 of file pfqn_simplex.h.

References line::pfqn::simplex::Radial< T >::G, gauss_legendre(), line::pfqn::simplex::Radial< T >::Lam, line::pfqn::simplex::Radial< T >::lJ, line::Matrix< T >::Matrix(), radial(), radial_logf(), and line::pfqn::simplex::Radial< T >::vbar.

Referenced by line::pfqn::pfqn_aghq(), radial(), and simplex_mode().

◆ radial_logf()

template<class T>

T line::pfqn::simplex::radial_logf	(	const T &	t,
		const std::vector< T > &	c,
		const std::vector< T > &	N,
		const std::vector< T > &	Z,
		std::size_t	M )

Log-integrand of the radial integral in t = log v, Jacobian included.

Definition at line 160 of file pfqn_simplex.h.

References radial_logf().

Referenced by radial(), and radial_logf().

◆ simplex_mode()

template<class T>

Mode< T > line::pfqn::simplex::simplex_mode	(	const Matrix< T > &	L,
		const std::vector< T > &	N,
		const std::vector< T > &	Z )

Mode and curvature of h(w) = log J(L'x(w)) + sum_i log x_i with J the exact radial integral.

The fixed point x = (1 + x.*(L*G))/vbar is the Z > 0 analogue of pfqn_le_fpi: integrating by parts gives sum_i x_i (L*G)_i = vbar - M, so the update is normalised by construction, and at Z = 0 it reduces to pfqn_le_fpi. The term in the second derivative of x(w) drops at the mode against sum_i x_i == 1.

Definition at line 343 of file pfqn_simplex.h.

References line::pfqn::simplex::Mode< T >::A, line::Matrix< T >::cols(), line::pfqn::simplex::Radial< T >::G, line::pfqn::simplex::Mode< T >::h0, line::pfqn::simplex::Radial< T >::Lam, line::pfqn::simplex::Mode< T >::ld, line::pfqn::simplex::Radial< T >::lJ, line::Matrix< T >::Matrix(), line::pfqn::pfqn_le_fpiZ(), radial(), line::Matrix< T >::rows(), simplex_mode(), line::pfqn::simplex::Radial< T >::vbar, and line::pfqn::simplex::Mode< T >::x.