Likelihood fitting with zmctp

Rasheedat Oladoja and Kabir Opeyemi Olorede

The package implements the complex triparametric Pearson (CTP) parent and a zero-inflated mixture. The parent follows Rodríguez-Avi et al. (2003), Statistical Papers 44, 67–88, doi:10.1007/s00362-002-0134-7. The zero-inflated extension is described in Chapter 3 of Oladoja’s 2021 Kwara State University thesis. The model’s original formulation and the numerical implementation are distinct contributions.

Distribution functions

The parent has

\[ f(k)=f(0)\frac{(a+ib)_k(a-ib)_k}{(\gamma)_k k!},\qquad f(0)=\frac{\Gamma(\gamma-a-ib)\Gamma(\gamma-a+ib)} {\Gamma(\gamma)\Gamma(\gamma-2a)}. \]

Arguments are scalar a, b >= 0 and gama > max(0, 2*a). The spelling gama is retained for compatibility. The normalising constant is computed independently of the counts requested. Negative or noninteger counts have probability zero. Missing vector entries propagate.

dctp(0:6,a=1,b=0.5,gama=6)
#> [1] 0.757021809 0.157712877 0.047877123 0.018452641 0.008329317 0.004206305
#> [7] 0.002310281
pctp(c(-Inf,0,2,Inf),1,0.5,6)
#> [1] 0.0000000 0.7570218 0.9626118 1.0000000
qctp(c(0,0.5,0.9,1),1,0.5,6)
#> [1]   0   0   1 Inf

Mean and variance need gama > 2*a+1 and gama > 2*a+2 respectively; normalisability alone does not imply finite moments. Terminating parents at b = 0 and nonpositive integer a are handled separately. A small positive b is not, by itself, evidence of an infinite variance.

Zero inflation

The mixture is

\[ P(Y=0)=\omega+(1-\omega)f(0),\qquad P(Y=k)=(1-\omega)f(k),\quad k>0. \]

For the same nondegenerate parent, every positive omega increases the probability of zero. There is no change from deflation to inflation at 0.5. At omega = 0 the distribution is CTP, and at omega = 1 it is concentrated at zero. Negative mixing probabilities are not implemented.

c(parent=dctp(0,1,0.5,6),mixture=dzictp(0,1,0.5,6,0.1))
#>    parent   mixture 
#> 0.7570218 0.7813196
mean_zictp(1,0.5,6,0.3)
#> [1] 0.2916667
var_zictp(1,0.5,6,0.3)
#> [1] 0.8263889
set.seed(123)
rzictp(12,1,0.5,6,0.3)
#>  [1] 0 0 0 1 2 0 0 1 0 0 2 0

Likelihood fitting

Fitting is restricted to b > 0 and gama > max(0, 2*a+2). It retains transformed coordinates, uses multiple starting points and polishes the best solutions. The mixture probability is profiled analytically, including the exact omega = 0 boundary. All-zero samples do not identify the parent parameters and produce an informative error.

The following is a constructed frequency example; it is not one of the published application datasets.

x <- rep(0:4,c(40,12,6,2,1))
fit <- zictp.fit(x)
fit$estimates
#>             a             b          gama         omega 
#>  2.006727e+02 9.797449e-196  3.880699e+04  4.684350e-01
fit$logLik
#> [1] -61.53984
fit$diagnostics
#> $uncertainty
#> [1] "Near boundary; Wald uncertainty withheld"
#> 
#> $variance_margin
#> [1] 38403.65
#> 
#> $b_near_zero
#> [1] TRUE
#> 
#> $gama_near_zero
#> [1] FALSE
#> 
#> $omega_boundary
#> [1] FALSE
#> 
#> $parent_logLik
#> [1] -62.15639
#> 
#> $parent_converged
#> [1] TRUE
#> 
#> $fit_domain
#> [1] "b > 0; gama > max(0, 2*a+2)"
#> 
#> $goodness_of_fit
#> [1] "Descriptive discrepancies; no calibrated p-values"
plot(fit,type="frequency")

A successful optimiser code is one diagnostic. Inspect the variance margin, small b, small gama, and the uncertainty status. Local standard errors are omitted when boundary checks fail or the Hessian is not sufficiently positive definite. These checks do not establish finite-sample validity of a Wald approximation.

fitted_freq includes zero-frequency counts and a final upper-tail bin. Its last x value denotes the start of that bin. pearson_chisq and the historically named wald_chisq are descriptive discrepancies. The package does not attach calibrated goodness-of-fit p-values to them.

Quantiles and runtime

max_q is an initial endpoint. A quantile search expands it until all requested probabilities are bracketed, subject to max_count (default 1,000,000). A failed search raises an error; it does not return a capped count. Infinite-support quantiles at probability one are infinite.

# This parent has F(k) = 1 - 1/(k+2).
qctp(0.9997,a=1,b=0,gama=3,max_q=1000)
#> [1] 3332

The default fitting argument method remains available and supplements the multistart search. omega_start is accepted for compatibility but does not affect the analytic profile. raw_optim$par now stores a, log_b, and log_margin; use estimates for natural parameters.

The research supplement contains the bootstrap studies, checkpoints and inference code. Bootstrap inference is not a public function in 0.1.4. The earlier study and the new package have separate code histories; the package checks do not constitute a new simulation calibration study.