Probability distributions

Families and parameterizations

Replace {family} in the operation names below with one of these names. PHP numeric arguments accept integers or floats; binomial size specifically requires an integer. All parameters must be finite and satisfy their domains.

Gamma’s second parameter is scale, exponential’s is rate. Negative-binomial size is a positive real dispersion/shape parameter, not a binomial trial count. The two NB forms are equivalent when prob = size/(size + mu). Poisson lambda zero, binomial size zero/prob zero/prob one, and NB prob one/mu zero represent point masses. Parameter conversion in the mean-form NB avoids overflowing size + mu.

Scalar and batch operations

Continuous families expose pdf/logpdf; discrete families expose pmf/logpmf. Every operation has an _array function taking a PHP array as its first argument and the same parameters/defaults, for example:

$probability = Phit\dist_normal_cdf(x: 1.96); // Approximately 0.975.
$tails = Phit\dist_t_sf_array(values: [1, 2, 3], df: 12);
$logMass = Phit\dist_poisson_logpmf_array(counts: [0, 2, 5], lambda: 3.2);
$interval = Phit\dist_gamma_quantile_array(probabilities: [.025, .975], shape: 2, scale: 3);

Scalar continuous evaluation uses argument x; discrete uses k. Batch continuous evaluation uses values; discrete uses counts; all batch quantiles use probabilities. Parameter names are part of the contract because PHP supports named arguments. The generated stubs/phit.php lists every exact signature. Loading the extension makes all functions available globally.

Batch evaluation constructs the distribution once, then streams input values through the same core evaluation policy. Output is a packed list in input iteration order; keys are not retained. Empty input returns [], but parameters are still validated. An invalid element rejects the entire call. No missing rows are dropped. The benchmarks quantify when batching saves PHP boundary overhead and when numerical evaluation dominates.

Tails and likelihoods

Use survival functions for upper tails rather than subtracting a rounded CDF from one. Use log functions directly rather than taking the log of a probability that may already have underflowed.

Phit\dist_normal_sf(10);     // About 7.62e-24; 1 - Phit\dist_normal_cdf(10) is zero.
Phit\dist_normal_logsf(40);  // About -804.61, even though the ordinary SF is zero.
$logLikelihood = array_sum(Phit\dist_poisson_logpmf_array([2, 5, 3, 0], lambda: 2.5));

A density is not a probability at a point and may exceed one. Sum log densities or log masses for independent observations; likelihood underflow is avoided this way. Log probabilities at impossible events are -INF. PHP JSON encoding with JSON_THROW_ON_ERROR rejects nonfinite results; callers must choose their own serialization policy rather than silently replacing numerical values.

The runnable examples/distributions.php demonstrates a two-sided t test, a mean confidence interval, a Poisson likelihood, a central probability, extreme upper tails and equivalent NB parameterizations. The sample data are hypothetical and the computations illustrate the API.

Boundaries and errors

Evaluation points may be ±INF, but NAN raises 1003. Null, numeric strings, booleans and nested array elements are not numeric inputs; malformed batch elements raise 1001. Scalar argument-type failures follow Zend’s runtime conventions. There is no distribution missing option. Nonfinite parameters raise 1003; invalid domains and quantile probabilities outside [0,1] raise 1001. Binomial size above 2^53 raises 1007. A NAN numerical result raises 1006; it is not returned as a value. Parameters so extreme that the incomplete beta or gamma cannot be evaluated also raise 1006 promptly: a beta-type sum \(a + b\) beyond the double range, a beta-type parameter above \(10^{300}\) evaluated within \(10^{-300}\) of an end of the support, or a gamma-type shape or point above \(10^{300}\) (for example dist_f_cdf(2, 3, PHP_FLOAT_MAX) or dist_gamma_cdf(1e30, 1e30)). Earlier builds could hang or panic on such input.

Outside support, density/mass is zero and log density/mass is -INF. Boundary densities can be infinite, for example a beta shape below one at a support endpoint. Discrete PMFs at non-integer counts are zero, while CDF/SF use the floor of the evaluation point. SF always means strictly greater than that count.

Continuous quantile endpoints are the support bounds, including infinities. Discrete quantiles return the smallest support value whose CDF reaches the requested probability, using a documented 64 * EPSILON relative comparison allowance on the smaller tail to accommodate rounding at CDF steps. Endpoint p=0 follows the conventional lower support bound; p=1 follows the upper support bound for ordinary distributions. A point mass returns its actual mass location at p=1. This differs from R/SciPy in some degenerate cases. Returned count quantiles are PHP floats; precision is limited by binary64’s exact integer range.

Numerical implementation and dependency audit

The raw statrs audit evaluates the saved R cases with statrs 0.19, using log(CDF)/log(SF) because that crate has no log-tail API. In these cases it loses small tails, has roughly 4.8e-4 relative binomial CDF error at size one million, and suffers inaccurate extreme quantiles. Those findings motivated a separate distribution module while retaining existing estimator numerical paths.

The runtime pins MIT-licensed cdflib = 0.4.3; GPL nmath ports were excluded by the dependency audit. statrs remains the estimator dependency and a source of special functions/seeds where appropriate. The implementation forms both tails independently and uses custom log-scale expansions below 1e-280. Continuous quantiles use log-tail Newton steps protected by brackets. No R/Python/BLAS runtime dependency is introduced. Numerical failure is explicit; this is not proof of accuracy for every finite argument or arbitrarily extreme shape parameter.

Reference evidence and release scope

The executed comparison records R, SciPy and 60-digit mpmath references, documented boundary differences and quantity-specific tolerances. Scalar/batch identities and validation run in the offline release gate. Performance evidence covers 1,000, 100,000 and one million values, including peak process memory and regression preservation.

Random sampling uses explicitly seeded dist_{family}_rng(n, seed, ...parameters) calls with local ChaCha20 state. The sampling contract documents reproducibility, validation and independent reference vectors. The readiness assessment records this implementation, the estimator reuse experiment and remaining publication/platform gates.