comp_osa.hpp File Reference

One-step-ahead (OSA) residual helpers for composition data. More...

Go to the source code of this file.

Functions
template<class Type>
vector< int >	osa_order (vector< Type > k)
template<class Type>
Type	osa_squeeze (Type u)
template<class Type>
Type	osa_squeeze_cdf (Type u)
	Map a CDF onto the open interval before its logarithm is taken.
template<class Type>
Type	osa_pbinom (Type x, Type n, Type p)
	Conditional binomial CDF P(X <= x), X ~ Binomial(n, p), at fractional x.
template<class Type>
Type	osa_norm_cdf_terms (Type x, Type mu, Type sd, Type cdf_lower, Type cdf_upper)
	One-step-ahead CDF terms for a normal observation.
template<class Type>
Type	osa_trunc_norm_cdf_terms (Type x, Type mu, Type sd, Type cdf_lower, Type cdf_upper)
	One-step-ahead CDF terms for a normal left-truncated at zero.
template<class Type>
Type	dmultinom_osa (vector< Type > x, vector< Type > p, data_indicator< vector< Type >, Type > keep, int do_osa, int give_log, int do_cdf=0)
	Keep-aware multinomial log-density via conditional binomials.
template<class Type>
Type	ddirmultinom_osa (vector< Type > obs, vector< Type > alpha, data_indicator< vector< Type >, Type > keep, int do_osa, int give_log)
	Keep-aware Dirichlet-multinomial log-density via conditional beta-binomials.

Detailed Description

One-step-ahead (OSA) residual helpers for composition data.

TMB's oneStepPredict() needs every observation in a flat obsvec, its likelihood term multiplied by its own keep element, which multivariate discrete composition data cannot supply directly. Following Trijoulet et al. (2023) and WHAM (src/age_comp_osa.hpp), a composition of A bins is decomposed into A-1 univariate conditionals – binomials for the multinomial, beta-binomials for the Dirichlet-multinomial – each gated by its keep element. The last bin is fixed by the sum-to-N constraint, contributes nothing, and is reported as NA on the R side. With every keep at 1 the conditional log-densities sum to the joint log-density, so OSA changes how the likelihood is decomposed, not its value. The _osa suffix and the data_indicator argument keep these clear of TMB's dmultinom() and of the ddirmultinom() in helper_functions.hpp that ordinary fitting uses.

The file also holds the conditional CDFs oneStepPredict(method = "cdf") reads – continuous binomial for a composition bin, normal / left-truncated normal for an aggregate series. That residual is Phi^-1(F(x)), formed without a conditional mean, so it cannot leave the support of a bounded observation.

Shared contract of the CDF terms below – the two _cdf_terms helpers and dmultinom_osa()'s do_cdf path. Each contributes to the log-likelihood keep.cdf_lower times log F plus keep.cdf_upper times log of its complement, which is 0 whenever the cdf gates are 0: every evaluation except the two oneStepPredict(method = "cdf") makes per observation. The terms are added only when the caller asks (do_cdf, from osa_mode == 2), so ordinary fitting and the Gaussian OSA methods are unchanged and pay nothing for them. Every F is mapped onto the open interval by osa_squeeze_cdf() first, because an exactly 0 or 1 CDF puts a -Inf on the tape and a -Inf times a zero gate is a NaN.

Definition in file comp_osa.hpp.

Function Documentation

ddirmultinom_osa()

template<class Type>

Type ddirmultinom_osa	(	vector< Type >	obs,
		vector< Type >	alpha,
		data_indicator< vector< Type >, Type >	keep,
		int	do_osa,
		int	give_log )

Keep-aware Dirichlet-multinomial log-density via conditional beta-binomials.

No do_cdf counterpart: a beta-binomial has no elementary CDF, and the reference route – summing the pmf over 0..floor(x) (contrib/OSA_multivariate_dists-main/distr.hpp, pbetabinom) – is a step function of a fractional count costing O(x) beta functions per bin, thousands of times the multinomial path at CEATTLE sample sizes. So a Dirichlet-multinomial fleet cannot be residualized with oneStepPredict(method = "cdf"); osa_residuals() sends it to a Gaussian method instead.

Parameters

obs	Observed counts (length A).
alpha	Dirichlet-multinomial concentration parameters (length A).
keep	Data indicator for these A bins.
do_osa	1 = conditional beta-binomial decomposition; 0 = plain D-M.
give_log	1 = return the log-density.

Definition at line 252 of file comp_osa.hpp.

References ddirmultinom(), and osa_order().

Referenced by comp_row_nll().

dmultinom_osa()

template<class Type>

Type dmultinom_osa	(	vector< Type >	x,
		vector< Type >	p,
		data_indicator< vector< Type >, Type >	keep,
		int	do_osa,
		int	give_log,
		int	do_cdf = 0 )

Keep-aware multinomial log-density via conditional binomials.

Parameters

x	Observed counts (length A; may be non-integer, e.g. Neff * prop).
p	Predicted proportions (length A; normalized internally).
keep	Data indicator for these A bins (from keep.segment(...)).
do_osa	1 = conditional-binomial decomposition; 0 = plain multinomial.
give_log	1 = return the log-density.
do_cdf	1 = also add the CDF terms (osa_mode == 2); 0 skips the pbeta.

Definition at line 159 of file comp_osa.hpp.

References osa_order(), osa_pbinom(), osa_squeeze(), and osa_squeeze_cdf().

Referenced by comp_row_nll().

osa_norm_cdf_terms()

template<class Type>

Type osa_norm_cdf_terms	(	Type	x,
		Type	mu,
		Type	sd,
		Type	cdf_lower,
		Type	cdf_upper )

One-step-ahead CDF terms for a normal observation.

Every aggregate series in CEATTLE is normal on the scale obsvec stores it in – log for lognormal catch and index, natural scale for a "Normal" index, the whitened L^-1 scale for "MVN" – so its conditional CDF is pnorm of the standardized residual.

Parameters

x	Observation, on the scale it is stored in.
mu	Predicted value on that same scale.
sd	Standard deviation on that same scale.
cdf_lower,cdf_upper	The observation's keep.cdf_lower/keep.cdf_upper.

Definition at line 117 of file comp_osa.hpp.

References osa_squeeze_cdf().

osa_order()

template<class Type>

vector< int > osa_order

(

vector< Type >

k

)

Definition at line 36 of file comp_osa.hpp.

Referenced by ddirmultinom_osa(), and dmultinom_osa().

osa_pbinom()

template<class Type>

Type osa_pbinom	(	Type	x,
		Type	n,
		Type	p )

Conditional binomial CDF P(X <= x), X ~ Binomial(n, p), at fractional x.

Written through the regularized incomplete beta, 1 - I_p(x + 1, n - x), which equals the binomial CDF at whole counts and stays smooth at the fractional ones composition data carry: a bin count is proportion * effective sample size, an integer only by coincidence. Same expression as the Trijoulet et al. (2023) reference implementation TMB ships (contrib/OSA_multivariate_dists-main/distr.hpp, dists::pbinom).

Parameters

x	Count in this bin.
n	Count still unplaced across this bin and every later one.
p	Probability of landing in this bin rather than a later one.

Definition at line 95 of file comp_osa.hpp.

Referenced by dmultinom_osa().

osa_squeeze()

template<class Type>

Type osa_squeeze

(

Type

u

)

Definition at line 55 of file comp_osa.hpp.

Referenced by dmultinom_osa().

osa_squeeze_cdf()

template<class Type>

Type osa_squeeze_cdf

(

Type

u

)

Map a CDF onto the open interval before its logarithm is taken.

Four machine epsilons wide where osa_squeeze() uses one: oneStepPredict recovers F from 1/(1 + exp(nlcdf.lower - nlcdf.upper)) in double precision, and at one epsilon that denominator sits on the round-half-to-even tie above 1 – where it rounds to 1 so does F, and the residual is infinite. Measured on BS2017SS: 6 of 4538 composition residuals infinite at one epsilon, none at four; the cost is where |residual| saturates – 8.04 against 8.13 – a tail no Q-Q plot resolves.

It maps rather than clamps, so a caller that forms F by subtracting two CDFs relies on it to absorb the cancellation: osa_trunc_norm_cdf_terms()'s numerator can round to about -2e-16 at an observation near zero, which this shift carries back above 0 before the log. Narrowing it would return NaN there.

Definition at line 76 of file comp_osa.hpp.

Referenced by dmultinom_osa(), osa_norm_cdf_terms(), and osa_trunc_norm_cdf_terms().

osa_trunc_norm_cdf_terms()

template<class Type>

Type osa_trunc_norm_cdf_terms	(	Type	x,
		Type	mu,
		Type	sd,
		Type	cdf_lower,
		Type	cdf_upper )

One-step-ahead CDF terms for a normal left-truncated at zero.

With z = (x - mu)/sd and m = mu/sd, the truncated CDF is F(x) = [Phi(z) - Phi(-m)] / Phi(m) and its complement is 1 - F(x) = Phi(-z) / Phi(m), because Phi(m) + Phi(-m) = 1. Both are written in that form rather than as 1 - F, so the upper tail keeps its precision. This is exact wherever the truncation carries mass, which is what the Gaussian OSA methods cannot see: the truncation enters the density only through log Phi(m), a function of the prediction and not of the observation.

Parameters

x	Observation (natural scale, > 0).
mu	Predicted index (natural scale).
sd	Observation standard deviation (absolute, natural scale).
cdf_lower,cdf_upper	The observation's keep.cdf_lower/keep.cdf_upper.

Definition at line 139 of file comp_osa.hpp.

References osa_squeeze_cdf().