map_gamma.h File Reference

Autocorrelation decay rate of a MAP: the gamma of the geometric model rho(k) = rho0 * gamma^k with rho0 = (1 - 1/scv)/2. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/mam/map_moment.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/levmar.h"
#include "line/util/matrix.h"

Include dependency graph for map_gamma.h:

Go to the source code of this file.

Classes
struct	line::mam::MapGammaResult< T >
	Result of map_gamma, mirroring the MATLAB [GAMMA, RHO0] pair. More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
MapGammaResult< T >	line::mam::map_gamma_full (const Map< T > &m, long limit=1000)
template<class T>
T	line::mam::map_gamma (const Map< T > &m, long limit=1000)
	Autocorrelation decay rate of a MAP (map_gamma.m).

Detailed Description

Autocorrelation decay rate of a MAP: the gamma of the geometric model rho(k) = rho0 * gamma^k with rho0 = (1 - 1/scv)/2.

Templated port of matlab/lib/kpctoolbox/map/map_gamma.m, cross-checked against jar/src/main/java/jline/api/mam/Map_gamma.java.

Order one is Poisson and has decay rate zero. Order two has a genuinely geometric acf, so the rate is the exact ratio acf(2)/acf(1); a vanishing acf(1) means the MAP degenerates to a phase-type renewal process and the rate is again zero. Only above order two is a fit needed, and there both references evaluate the acf on the ten lags 1, 1+limit/10, ..., and regress.

WHY THE FIT IS ROBUST AND NOT PLAIN LEAST SQUARES: MATLAB uses nlinfit with RobustWgtFun 'fair'. On a higher-order MAP the acf is a sum of geometric terms, so the short lags sit far off any single geometric and would otherwise dominate the fit and bias gamma low. The reproduction here is the JAR's: an ordinary fit, the leverage of its Jacobian held fixed, then iteratively reweighted fits with w = 1/(1 + |r_adj|/(1.4 sigma)) and sigma the MAD estimate floored against the spread of the response. Dropping the robustness is not a simplification, it changes the answer.

ARITHMETIC: unlike most of api/mam this is a floating-point algorithm. The lags enter as real exponents and the MAD needs an order statistic, so the regression runs in double whatever T is; only the acf, mean and second moment that feed it are computed in T. Orders one and two return an exact T.

Definition in file map_gamma.h.