map_bernstein.h File Reference

Acyclic phase-type approximation of an arbitrary density by Bernstein exponentials. More...

#include <cmath>
#include <cstddef>
#include <functional>
#include "line/api/mam/map_moment.h"
#include "line/api/mam/map_transform.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for map_bernstein.h:

Go to the source code of this file.

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
Map< T >	line::mam::map_bernstein (const std::function< double(double)> &f, unsigned order=20)
	Acyclic phase-type approximation of an arbitrary density by Bernstein exponentials.
template<class T>
Map< T >	line::mam::aph_bernstein (const std::function< double(double)> &f, unsigned order=20)
	The reference's (D0, D1) spelling of the same fit, kept because the JAR exposes it under this name (Aph_bernstein) and callers ported from it expect the pair rather than a Map.

Detailed Description

Acyclic phase-type approximation of an arbitrary density by Bernstein exponentials.

Templated port of matlab/lib/kpctoolbox/map/map_bernstein.m and jar/src/main/java/jline/api/mam/Aph_bernstein.java, after Horvath and Vicario, "Construction of Phase Type Distributions by Bernstein Exponentials", EPEW 2023, 201-215.

The construction evaluates the target density at the n Bernstein nodes x_i = -log(i/n) and uses the values as the entry law of the Erlang cascade T = diag(-1..-n) + superdiag(1..n-1), whose i-th phase has an Erlang(i,i) residual time. The result is a renewal process, D1 = -T 1 alpha, so the fit is a PH and not a general MAP; the caller rescales it to the target mean with map_scale, since the construction pins the time unit at one.

WHICH REFERENCE. The MATLAB and JAR versions differ and MATLAB is the ground truth here: it skips a node where the density is not finite and positive, renormalizes alpha (which the bare 1/(i c) weights only satisfy up to the skipped nodes), and falls back to an Erlang-n of unit mean when the normalization constant is not usable at all. The JAR divides by c unconditionally, so a density that underflows at some node – a Pareto below its scale, a Uniform outside its support, any bounded law – gives it a NaN generator instead of a fit. That path is exactly the one sn_nonmarkov_toph takes for Uniform and Pareto service.

ARITHMETIC: transcendental, since the nodes are logarithms of i/n.

Definition in file map_bernstein.h.