bMax-free closed forms for the LCD / modified Cramer-von-Mises distance of an even-dimensional standard normal target More...

#include <cassert>
#include <cmath>
#include <cstddef>
#include "lcd_even_closed_form.h"

Functions
double	lcd_binf_b0_zero (size_t d, double e0, double ei0, double c)
	B_{0,d}(0, c), the lower-limit base case of the attraction term.

double	lcd_binf_b0_zero_deriv (size_t d, double e0, double ei0, double c)
	d/dc B_{0,d}(0, c)

double	lcd_binf_bkk_zero (size_t k, double c)
	B_{k,k}(0, c), the lower-limit attraction term.

double	lcd_binf_bkk_zero_deriv (size_t k, double c)
	d/dc B_{k,k}(0, c)

double	lcd_binf_s (size_t k, double c)
	S_k(c), the pole-cancelling companion of B_{k,k}(0, c)

double	lcd_binf_s_deriv (size_t k, double c)
	S_k'(c)

double	lcd_binf_phi (size_t k, double c)
	phi_k(c), the per-sample attraction potential of the bMax-free distance

double	lcd_binf_phi_deriv (size_t k, double c)
	phi_k'(c)

double	lcd_binf_constant (size_t k)
	K_k, the sample-independent constant of the bMax-free distance.

Variables
constexpr double	lcdBinfEulerGamma = 0.5772156649015328606
	Euler-Mascheroni constant.

Detailed Description

bMax-free closed forms for the LCD / modified Cramer-von-Mises distance of an even-dimensional standard normal target

lcd_even_closed_form.h evaluates the mCvM distance for a finite integration bound bMax. This header takes the limit bMax -> infinity analytically, so bMax disappears from the problem entirely.

With N = 2k, samples s_i, weights w_i summing to 1, c_i = ||s_i||^2 and T_ij = ||s_i - s_j||^2 the limit distance is

D = K_k - 2 * sum_i w_i * phi_k(c_i)
        + sum_{i,j, T_ij > 0} w_i w_j * (T_ij / 8) * ln(T_ij / 4)

Note: PRECONDITION: the limit exists only for a ZERO-MEAN sample set, sum_i w_i s_i = 0. For a non-zero mean the true distance diverges like ||mu||^2 / 4 * ln(bMax^2); the formulas here still return a finite number in that case, but that number is NOT the distance. gm_to_dirac_even_binf removes the hazard by optimising on the zero-mean manifold and projecting before every evaluation.; PRECONDITION: the target must be the STANDARD normal N(0, I) and N = 2k must be EVEN.

Warning: c == 0 IS A DOMAIN ERROR, not merely inaccurate: it reaches lcd_ei(0), whose GSL_EDOM result GSL's default error handler turns into an abort long before lcd_ei() can map it to 0.0. Every entry point here asserts c > 0, and gm_to_dirac_even_binf refuses a sample set that projects any component onto the origin.

Function Documentation

◆ lcd_binf_b0_zero()

double lcd_binf_b0_zero	(	size_t	d,
		double	e0,
		double	ei0,
		double	c
	)

inline

B_{0,d}(0, c), the lower-limit base case of the attraction term.

With e0 = exp(-c/2) and Ei0 = Ei(-c/2),

B_{0,0}(0,c) = e0/4 + (c/8) * Ei0

B_{0,1}(0,c) = -Ei0/4

B_{0,d}(0,c) = e0 * sum_{m=2..d} (d-2)! * 2^(d-m-1)
                                 / ( (m-2)! * c^(d-m+1) )   for d >= 2

The d >= 2 sum is accumulated downwards from m = d, where the term is 1 / (2c), using term_{m-1} = term_m * 2 (m - 2) / c. That is the same recurrence lcd_delta_b0() uses and it keeps both (d-2)! and c^(d-1) out of the intermediate values.

Parameters

d	index of the base case
e0	exp(-c / 2)
ei0	Ei(-c / 2)
c	squared norm of the sample, must be > 0

Returns: B_{0,d}(0, c)

◆ lcd_binf_b0_zero_deriv()

double lcd_binf_b0_zero_deriv	(	size_t	d,
		double	e0,
		double	ei0,
		double	c
	)

inline

d/dc B_{0,d}(0, c)

Using d/dc Ei(-c/2) = e0 / c,

dB_{0,0}/dc = Ei0 / 8

dB_{0,1}/dc = -e0 / (4c)

dB_{0,d}/dc = -e0/2 * T_d(c) + e0 * T_d'(c)                 for d >= 2

where T_d(c) is the d >= 2 sum of lcd_binf_b0_zero() without the e0 factor and T_d'(c) is its termwise derivative, i.e. term_m scaled by -(d - m + 1) / c. Both sums share the one downward recurrence.

Parameters

d	index of the base case
e0	exp(-c / 2)
ei0	Ei(-c / 2)
c	squared norm of the sample, must be > 0

Returns: d/dc B_{0,d}(0, c)

◆ lcd_binf_bkk_zero()

double lcd_binf_bkk_zero	(	size_t	k,
		double	c
	)

inline

B_{k,k}(0, c), the lower-limit attraction term.

B_{k,k}(0,c) = 2^-k * sum_{j=0..k} (-1)^j * binom(k,j) * B_{0,j}(0,c)

Parameters

k	half the dimension, N = 2k; must be >= 1
c	squared norm of the sample, must be > 0

Returns: B_{k,k}(0, c)

Note: standard-normal target and even N only

◆ lcd_binf_bkk_zero_deriv()

double lcd_binf_bkk_zero_deriv	(	size_t	k,
		double	c
	)

inline

d/dc B_{k,k}(0, c)

dB_{k,k}(0,c)/dc = 2^-k * sum_{j=0..k} (-1)^j * binom(k,j)

dB_{0,j}(0,c)/dc

Parameters

k	half the dimension, N = 2k; must be >= 1
c	squared norm of the sample, must be > 0

Returns: d/dc B_{k,k}(0, c)

Note: standard-normal target and even N only

◆ lcd_binf_constant()

double lcd_binf_constant ( size_t k )

inline

K_k, the sample-independent constant of the bMax-free distance.

K_k = k * sum_{i=2..k} 1/(2i) - 1/2 - k * gamma / 2

K_k does not depend on the samples, so it shifts the objective without moving the optimum. It is carried anyway so the reported value is the true distance.

Parameters

k	half the dimension, N = 2k; must be >= 1

Returns: K_k

◆ lcd_binf_phi()

double lcd_binf_phi	(	size_t	k,
		double	c
	)

inline

phi_k(c), the per-sample attraction potential of the bMax-free distance

phi_k(c) = (c + 2k)/8 * ln(c/4) + S_k(c) - 2^k * B_{k,k}(0,c)

The distance uses it as D = K_k - 2 * sum_i w_i * phi_k(c_i) + repulsion.

Parameters

k	half the dimension, N = 2k; must be >= 1
c	squared norm of the sample, must be > 0

Returns: phi_k(c)

◆ lcd_binf_phi_deriv()

double lcd_binf_phi_deriv	(	size_t	k,
		double	c
	)

inline

phi_k'(c)

phi_k'(c) = ln(c/4)/8 + (c + 2k)/(8c) + S_k'(c)

2^k * dB_{k,k}(0,c)/dc

Parameters

k	half the dimension, N = 2k; must be >= 1
c	squared norm of the sample, must be > 0

Returns: phi_k'(c)

◆ lcd_binf_s()

double lcd_binf_s	(	size_t	k,
		double	c
	)

inline

S_k(c), the pole-cancelling companion of B_{k,k}(0, c)

S_k(c) = sum_{j=2..k} (-1)^j * binom(k,j) * (j-2)! * 2^(j-3) / c^(j-1)

The sum is empty for k = 1. binom(k,j), (j-2)! * 2^(j-3) and c^(j-1) are each carried forward by their own recurrence, as lcd_delta_bkk_zero() does.

Parameters

k	half the dimension, N = 2k; must be >= 1
c	squared norm of the sample, must be > 0

Returns: S_k(c), 0 for k = 1

◆ lcd_binf_s_deriv()

double lcd_binf_s_deriv	(	size_t	k,
		double	c
	)

inline

S_k'(c)

S_k'(c) = sum_{j=2..k} (-1)^j * binom(k,j) * (j-2)! * 2^(j-3)

(1 - j) / c^j

Parameters

k	half the dimension, N = 2k; must be >= 1
c	squared norm of the sample, must be > 0

Returns: S_k'(c), 0 for k = 1

Functions

Variables

Detailed Description

Function Documentation

◆ lcd_binf_b0_zero()

◆ lcd_binf_b0_zero_deriv()

◆ lcd_binf_bkk_zero()

◆ lcd_binf_bkk_zero_deriv()

◆ lcd_binf_constant()

◆ lcd_binf_phi()

◆ lcd_binf_phi_deriv()

◆ lcd_binf_s()

◆ lcd_binf_s_deriv()