diff -dPNur stats-1.0.2/config.m4 trunk/config.m4 --- stats-1.0.2/config.m4 2006-05-31 19:24:26.000000000 +0200 +++ trunk/config.m4 2012-10-29 17:22:36.000000000 +0100 @@ -1,9 +1,9 @@ -dnl $Id: config.m4,v 1.2 2004/01/09 19:22:36 andrey Exp $ +dnl $Id: config.m4 256977 2008-04-08 16:13:17Z sfox $ dnl config.m4 for extension stats PHP_ARG_ENABLE(stats, whether to enable stats support, [ --enable-stats Enable statistics support]) if test "$PHP_STATS" != "no"; then - PHP_NEW_EXTENSION(stats, statistics.c com.c dcdflib.c ipmpar.c linpack.c randlib.c, $ext_shared) + PHP_NEW_EXTENSION(stats, php_stats.c com.c dcdflib.c ipmpar.c linpack.c randlib.c , $ext_shared) fi diff -dPNur stats-1.0.2/config.w32 trunk/config.w32 --- stats-1.0.2/config.w32 2006-05-31 19:24:26.000000000 +0200 +++ trunk/config.w32 2012-10-29 17:22:36.000000000 +0100 @@ -1,9 +1,9 @@ -// $Id: config.w32,v 1.2 2006/05/30 17:53:18 andrey Exp $ +// $Id: config.w32 256977 2008-04-08 16:13:17Z sfox $ // vim:ft=javascript ARG_ENABLE("stats", "statistics support", "no"); if (PHP_STATS != "no") { - EXTENSION('stats', 'statistics.c com.c dcdflib.c ipmpar.c linpack.c randlib.c fd_e_lgamma_r.c fd_e_log.c fd_k_cos.c fd_k_sin.c fd_w_lgamma.c'); + EXTENSION('stats', 'php_stats.c com.c dcdflib.c ipmpar.c linpack.c randlib.c fd_e_lgamma_r.c fd_e_log.c fd_k_cos.c fd_k_sin.c fd_w_lgamma.c'); } diff -dPNur stats-1.0.2/package.xml trunk/package.xml --- stats-1.0.2/package.xml 1970-01-01 01:00:00.000000000 +0100 +++ trunk/package.xml 2012-10-29 17:22:36.000000000 +0100 @@ -0,0 +1,104 @@ + + + stats + pecl.php.net + Extension with routines for statistical computation. + Extension that provides few dozens routines for statistical computation. + + + Andrey Hristov + andrey + andrey_php_net + yes + + 2006-05-31 + + + 1.0.2 + 1.0.2 + + + stable + stable + + PHP License + Added fd_lgamma from FDLibM, so now all functions are also available on Windows. + Previously some were disabled because Microsoft does not ship lgamma with their C library. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 4.0.0 + + + 1.4.0 + + + + stats + + + + + 1.0.0 + 1.0.0 + + + stable + stable + + 2005-11-24 + PHP License + Initial release. + + + + + 1.0.1 + 1.0.1 + + + stable + stable + + 2005-3-25 + PHP License + Fixed shared build. + + + + diff -dPNur stats-1.0.2/php_statistics.h trunk/php_statistics.h --- stats-1.0.2/php_statistics.h 2006-05-31 19:24:26.000000000 +0200 +++ trunk/php_statistics.h 1970-01-01 01:00:00.000000000 +0100 @@ -1,127 +0,0 @@ -/* - +----------------------------------------------------------------------+ - | PHP Version 5 | - +----------------------------------------------------------------------+ - | Copyright (c) 1997-2004 The PHP Group | - +----------------------------------------------------------------------+ - | This source file is subject to version 3.0 of the PHP license, | - | that is bundled with this package in the file LICENSE, and is | - | available through the world-wide-web at the following url: | - | http://www.php.net/license/3_0.txt. | - | If you did not receive a copy of the PHP license and are unable to | - | obtain it through the world-wide-web, please send a note to | - | license@php.net so we can mail you a copy immediately. | - +----------------------------------------------------------------------+ - | Author: Andrey Hristov | - +----------------------------------------------------------------------+ -*/ - -/* $Id: php_statistics.h,v 1.5 2005/05/13 12:17:07 andrey Exp $ */ - -#ifndef PHP_STATISTICS_H -#define PHP_STATISTICS_H - -extern zend_module_entry stats_module_entry; -#define phpext_stats_ptr &stats_module_entry - -#ifdef PHP_WIN32 -#define PHP_STATS_API __declspec(dllexport) -#else -#define PHP_STATS_API -#endif - - -PHP_MINFO_FUNCTION(stats); - -PHP_FUNCTION(stats_bin_counts); -PHP_FUNCTION(stats_cdf_t); -PHP_FUNCTION(stats_cdf_normal); -PHP_FUNCTION(stats_cdf_gamma); -PHP_FUNCTION(stats_cdf_chisquare); -PHP_FUNCTION(stats_cdf_beta); -PHP_FUNCTION(stats_cdf_binomial); -PHP_FUNCTION(stats_cdf_noncentral_chisquare); -PHP_FUNCTION(stats_cdf_f); -PHP_FUNCTION(stats_cdf_noncentral_f); -PHP_FUNCTION(stats_cdf_noncentral_t); -PHP_FUNCTION(stats_cdf_negative_binomial); -PHP_FUNCTION(stats_cdf_poisson); -PHP_FUNCTION(stats_cdf_laplace); -PHP_FUNCTION(stats_cdf_cauchy); -PHP_FUNCTION(stats_cdf_logistic); -PHP_FUNCTION(stats_cdf_weibull); -PHP_FUNCTION(stats_cdf_uniform); -PHP_FUNCTION(stats_cdf_exponential); -PHP_FUNCTION(stats_rand_setall); -PHP_FUNCTION(stats_rand_getsd); -PHP_FUNCTION(stats_rand_gen_iuniform); -PHP_FUNCTION(stats_rand_gen_funiform); -PHP_FUNCTION(stats_rand_ignlgi); -PHP_FUNCTION(stats_rand_ranf); -PHP_FUNCTION(stats_rand_gen_beta); -PHP_FUNCTION(stats_rand_gen_chisquare); -PHP_FUNCTION(stats_rand_gen_exponential); -PHP_FUNCTION(stats_rand_gen_f); -PHP_FUNCTION(stats_rand_gen_gamma); -PHP_FUNCTION(stats_rand_gen_noncentral_chisquare); -PHP_FUNCTION(stats_rand_gen_noncenral_f); -PHP_FUNCTION(stats_rand_gen_normal); -PHP_FUNCTION(stats_rand_phrase_to_seeds); -PHP_FUNCTION(stats_rand_ibinomial); -PHP_FUNCTION(stats_rand_ibinomial_negative); -PHP_FUNCTION(stats_rand_gen_ipoisson); -PHP_FUNCTION(stats_rand_gen_noncentral_t); -PHP_FUNCTION(stats_rand_gen_t); -PHP_FUNCTION(stats_dens_normal); -PHP_FUNCTION(stats_dens_cauchy); -PHP_FUNCTION(stats_dens_laplace); -PHP_FUNCTION(stats_dens_logistic); -PHP_FUNCTION(stats_dens_beta); -PHP_FUNCTION(stats_dens_weibull); -PHP_FUNCTION(stats_dens_uniform); -PHP_FUNCTION(stats_dens_chisquare); -PHP_FUNCTION(stats_dens_t); -PHP_FUNCTION(stats_dens_gamma); -PHP_FUNCTION(stats_dens_exponential); -PHP_FUNCTION(stats_dens_f); -PHP_FUNCTION(stats_dens_pmf_binomial); -PHP_FUNCTION(stats_dens_pmf_poisson); -PHP_FUNCTION(stats_dens_pmf_negative_binomial); -PHP_FUNCTION(stats_dens_pmf_hypergeometric); -PHP_FUNCTION(stats_stat_powersum); -PHP_FUNCTION(stats_stat_innerproduct); -PHP_FUNCTION(stats_stat_independent_t); -PHP_FUNCTION(stats_stat_paired_t); -PHP_FUNCTION(stats_stat_percentile); -PHP_FUNCTION(stats_stat_correlation); -PHP_FUNCTION(stats_stat_binomial_coef); -PHP_FUNCTION(stats_stat_factorial); -PHP_FUNCTION(stats_absolute_deviation); -PHP_FUNCTION(stats_standard_deviation); -PHP_FUNCTION(stats_variance); -PHP_FUNCTION(stats_harmonic_mean); -PHP_FUNCTION(stats_skew); -PHP_FUNCTION(stats_kurtosis); -PHP_FUNCTION(stats_covariance); - - -#ifdef ZTS -#define STATS_D zend_stats_globals *stats_globals -#define STATS_G(v) (stats_globals->v) -#define STATS_FETCH() zend_stats_globals *stats_globals = ts_resource(stats_globals_id) -#else -#define STATS_D -#define STATS_G(v) (stats_globals.v) -#define STATS_FETCH() -#endif - -#endif /* PHP_STATISTICS_H */ - - -/* - * Local variables: - * tab-width: 4 - * c-basic-offset: 4 - * indent-tabs-mode: t - * End: - */ diff -dPNur stats-1.0.2/php_stats.c trunk/php_stats.c --- stats-1.0.2/php_stats.c 1970-01-01 01:00:00.000000000 +0100 +++ trunk/php_stats.c 2012-10-29 17:22:36.000000000 +0100 @@ -0,0 +1,3729 @@ +/* + +----------------------------------------------------------------------+ + | PHP Version 5 | + +----------------------------------------------------------------------+ + | Copyright (c) 1997-2004 The PHP Group | + +----------------------------------------------------------------------+ + | This source file is subject to version 3.0 of the PHP license, | + | that is bundled with this package in the file LICENSE, and is | + | available through the world-wide-web at the following url: | + | http://www.php.net/license/3_0.txt. | + | If you did not receive a copy of the PHP license and are unable to | + | obtain it through the world-wide-web, please send a note to | + | license@php.net so we can mail you a copy immediately. | + +----------------------------------------------------------------------+ + | Author: Andrey Hristov | + +----------------------------------------------------------------------+ +*/ + +/* $Id: php_stats.c 313529 2011-07-21 12:14:10Z derick $ */ + + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "php.h" +#include "php_stats.h" +#include "ext/standard/info.h" +#include "ext/standard/head.h" +#include +#include +#include +#include +#include "randlib.h" +#include "cdflib.h" + +#define STATS_PI 3.14159265358979323846 + + +#ifdef PHP_WIN32 +extern double fd_lgamma(double x); +#define lgamma fd_lgamma +#endif + +static double logistic_quantile(double p); +static double logistic_cdf(double x); +static double cauchy_quantile(double p); +static double cauchy_cdf(double x); +static double laplace_quantile(double p); +static double laplace_cdf(double x); +static double exponential_quantile(double p); +static double exponential_cdf(double x); +static double binom(double x, double n); + +zend_function_entry statistics_functions[] = { + PHP_FE(stats_cdf_t, NULL) + PHP_FE(stats_cdf_normal, NULL) + PHP_FE(stats_cdf_gamma, NULL) + PHP_FE(stats_cdf_chisquare, NULL) + PHP_FE(stats_cdf_beta, NULL) + PHP_FE(stats_cdf_binomial, NULL) + PHP_FE(stats_cdf_noncentral_chisquare,NULL) + PHP_FE(stats_cdf_f, NULL) + PHP_FE(stats_cdf_noncentral_f, NULL) + PHP_FE(stats_cdf_noncentral_t, NULL) + PHP_FE(stats_cdf_negative_binomial, NULL) + PHP_FE(stats_cdf_poisson, NULL) + PHP_FE(stats_cdf_laplace, NULL) + PHP_FE(stats_cdf_cauchy, NULL) + PHP_FE(stats_cdf_logistic, NULL) + PHP_FE(stats_cdf_weibull, NULL) + PHP_FE(stats_cdf_uniform, NULL) + PHP_FE(stats_cdf_exponential, NULL) + PHP_FE(stats_rand_setall, NULL) + PHP_FE(stats_rand_getsd, NULL) + PHP_FE(stats_rand_gen_iuniform, NULL) + PHP_FE(stats_rand_gen_funiform, NULL) + PHP_FE(stats_rand_ignlgi, NULL) + PHP_FE(stats_rand_ranf, NULL) + PHP_FE(stats_rand_gen_beta, NULL) + PHP_FE(stats_rand_gen_chisquare, NULL) + PHP_FE(stats_rand_gen_exponential, NULL) + PHP_FE(stats_rand_gen_f, NULL) + PHP_FE(stats_rand_gen_gamma, NULL) + PHP_FE(stats_rand_gen_noncentral_chisquare,NULL) + PHP_FE(stats_rand_gen_noncenral_f, NULL) + PHP_FE(stats_rand_gen_normal, NULL) + PHP_FE(stats_rand_phrase_to_seeds, NULL) + PHP_FE(stats_rand_ibinomial, NULL) + PHP_FE(stats_rand_ibinomial_negative,NULL) + PHP_FE(stats_rand_gen_ipoisson, NULL) + PHP_FE(stats_rand_gen_noncentral_t, NULL) + PHP_FE(stats_rand_gen_t, NULL) + PHP_FE(stats_dens_normal, NULL) + PHP_FE(stats_dens_cauchy, NULL) + PHP_FE(stats_dens_laplace, NULL) + PHP_FE(stats_dens_logistic, NULL) + PHP_FE(stats_dens_beta, NULL) + PHP_FE(stats_dens_weibull, NULL) + PHP_FE(stats_dens_uniform, NULL) + PHP_FE(stats_dens_chisquare, NULL) + PHP_FE(stats_dens_t, NULL) + PHP_FE(stats_dens_gamma, NULL) + PHP_FE(stats_dens_exponential, NULL) + PHP_FE(stats_dens_f, NULL) + PHP_FE(stats_dens_pmf_binomial, NULL) + PHP_FE(stats_dens_pmf_poisson, NULL) + PHP_FE(stats_dens_pmf_negative_binomial,NULL) + PHP_FE(stats_dens_pmf_hypergeometric, NULL) + PHP_FE(stats_stat_powersum, NULL) + PHP_FE(stats_stat_innerproduct, NULL) + PHP_FE(stats_stat_independent_t, NULL) + PHP_FE(stats_stat_paired_t, NULL) + PHP_FE(stats_stat_percentile, NULL) + PHP_FE(stats_stat_correlation, NULL) + PHP_FE(stats_stat_binomial_coef, NULL) + PHP_FE(stats_stat_factorial, NULL) + PHP_FE(stats_standard_deviation, NULL) + PHP_FE(stats_absolute_deviation, NULL) + PHP_FE(stats_variance, NULL) + PHP_FE(stats_harmonic_mean, NULL) + PHP_FE(stats_skew, NULL) + PHP_FE(stats_kurtosis, NULL) + PHP_FE(stats_covariance, NULL) + {NULL, NULL, NULL} +}; + +zend_module_entry stats_module_entry = { + STANDARD_MODULE_HEADER, + "stats", + statistics_functions, + NULL, + NULL, + NULL, + NULL, + PHP_MINFO(stats), + PHP_STATS_VERSION, + STANDARD_MODULE_PROPERTIES, +}; + +#ifdef COMPILE_DL_STATS +ZEND_GET_MODULE(stats) +#endif + + +PHP_MINFO_FUNCTION(stats) +{ + php_info_print_table_start(); + php_info_print_table_header(2, "Statistics Support", "enabled"); + php_info_print_table_row(2, "Version", PHP_STATS_VERSION); + php_info_print_table_end(); +} + + + +/* Numbers are always smaller than strings int this function as it + * anyway doesn't make much sense to compare two different data types. + * This keeps it consistant and simple. + * + * This is not correct any more, depends on what compare_func is set to. + */ +static int stats_array_data_compare(const void *a, const void *b TSRMLS_DC) +{ + Bucket *f; + Bucket *s; + zval result; + zval *first; + zval *second; + + f = *((Bucket **) a); + s = *((Bucket **) b); + + first = *((zval **) f->pData); + second = *((zval **) s->pData); + + if (numeric_compare_function(&result, first, second TSRMLS_CC) == FAILURE) { + return 0; + } + + if (Z_TYPE(result) == IS_DOUBLE) { + if (Z_DVAL(result) < 0) { + return -1; + } else if (Z_DVAL(result) > 0) { + return 1; + } else { + return 0; + } + } + + convert_to_long(&result); + + if (Z_LVAL(result) < 0) { + return -1; + } else if (Z_LVAL(result) > 0) { + return 1; + } + + return 0; +} + + + +/**************************************/ +/* Cumulative Distributions Functions */ +/**************************************/ + +/****************************************************** + Cumulative Distribution Function + T distribution + + Function + + + Calculates any one parameter of the t distribution given + values for the others. + + Arguments + + WHICH --> Integer indicating which argument + values is to be calculated from the others. + Legal range: 1..3 + iwhich = 1 : Calculate P and Q from T and DF + iwhich = 2 : Calculate T from P,Q and DF + iwhich = 3 : Calculate DF from P,Q and T + + P <--> The integral from -infinity to t of the t-density. + Input range: (0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + T <--> Upper limit of integration of the t-density. + Input range: ( -infinity, +infinity). + Search range: [ -1E100, 1E100 ] + + DF <--> Degrees of freedom of the t-distribution. + Input range: (0 , +infinity). + Search range: [1e-100, 1E10] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.5.27 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to reduce the computation + of the cumulative distribution function to that of an incomplete + beta. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. +******************************************************/ + +/* {{{ proto float stats_cdf_t(float par1, float par2, int which) + Calculates any one parameter of the T distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_t) +{ + double arg1; + double arg2; + double df; + double bound; + double p; + double q; + double t; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddl", &arg1, &arg2, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 3) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); + RETURN_FALSE; + } + + if (which < 3 ) { + df = arg2; + } else { + t = arg2; + } + if (which == 1) { + t = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdft((int *)&which, &p, &q, &t, &df, &status, &bound); + + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(t); + case 3: RETURN_DOUBLE(df); + } + RETURN_FALSE; /* should never be reached */ +} +/* }}} */ + +/********************************************************************* + Cumulative Distribution Function NORmal distribution + + Calculates any one parameter of the normal + distribution given values for the others. + + + Arguments + + + WHICH --> Integer indicating which of the next parameter + values is to be calculated using values of the others. + Legal range: 1..4 + iwhich = 1 : Calculate P and Q from X,MEAN and SD + iwhich = 2 : Calculate X from P,Q,MEAN and SD + iwhich = 3 : Calculate MEAN from P,Q,X and SD + iwhich = 4 : Calculate SD from P,Q,X and MEAN + + P <--> The integral from -infinity to X of the normal density. + Input range: (0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + X < --> Upper limit of integration of the normal-density. + Input range: ( -infinity, +infinity) + + MEAN <--> The mean of the normal density. + Input range: (-infinity, +infinity) + + SD <--> Standard Deviation of the normal density. + Input range: (0, +infinity). + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + A slightly modified version of ANORM from + + Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN + Package of Special Function Routines and Test Drivers" + acm Transactions on Mathematical Software. 19, 22-32. + + is used to calulate the cumulative standard normal distribution. + + The rational functions from pages 90-95 of Kennedy and Gentle, + Statistical Computing, Marcel Dekker, NY, 1980 are used as + starting values to Newton's Iterations which compute the inverse + standard normal. Therefore no searches are necessary for any + parameter. + + For X < -15, the asymptotic expansion for the normal is used as + the starting value in finding the inverse standard normal. + This is formula 26.2.12 of Abramowitz and Stegun. + + Note + + The normal density is proportional to + exp( - 0.5 * (( X - MEAN)/SD)**2) +***********************************************************************/ + +/* {{{ proto float stats_stat_gennch(float par1, float par2, float par3, int which) + Calculates any one parameter of the normal distribution given values for thee others. */ +PHP_FUNCTION(stats_cdf_normal) +{ + double arg1; + double arg2; + double arg3; + double sd; + double bound; + double p; + double q; + double x; + double mean; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + sd = arg3; + } else { + mean = arg3; + } + + if (which < 3) { + mean = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfnor((int *)&which, &p, &q, &x, &mean, &sd, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(mean); + case 4: RETURN_DOUBLE(sd); + } + RETURN_FALSE; /* should never be reached */ +} +/* }}} */ + + +/********************************************************************* + Cumulative Distribution Function + GAMma Distribution + + + Function + + + Calculates any one parameter of the gamma + distribution given values for the others. + + + Arguments + + + WHICH --> Integer indicating which of the next four argument + values is to be calculated from the others. + Legal range: 1..4 + iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE + iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE + iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE + iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE + + P <--> The integral from 0 to X of the gamma density. + Input range: [0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + X <--> The upper limit of integration of the gamma density. + Input range: [0, +infinity). + Search range: [0,1E100] + + SHAPE <--> The shape parameter of the gamma density. + Input range: (0, +infinity). + Search range: [1E-100,1E100] + + SCALE <--> The scale parameter of the gamma density. + Input range: (0, +infinity). + Search range: (1E-100,1E100] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + 10 if the gamma or inverse gamma routine cannot + compute the answer. Usually happens only for + X and SHAPE very large (gt 1E10 or more) + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + + Method + + + Cumulative distribution function (P) is calculated directly by + the code associated with: + + DiDinato, A. R. and Morris, A. H. Computation of the incomplete + gamma function ratios and their inverse. ACM Trans. Math. + Softw. 12 (1986), 377-393. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + + + Note + + + + The gamma density is proportional to + T**(SHAPE - 1) * EXP(- SCALE * T) +**************************************************************************/ +/* {{{ proto float stats_cdf_gamma(float par1, float par2, float par3, int which) + Calculates any one parameter of the gamma distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_gamma) +{ + double arg1; + double arg2; + double arg3; + double bound; + double p; + double q; + double x; + double shape; + double scale; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + scale = arg3; + } else { + shape = arg3; + } + + if (which < 3) { + shape = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfgam((int *)&which, &p, &q, &x, &shape, &scale, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(shape); + case 4: RETURN_DOUBLE(scale); + } + RETURN_FALSE; /* should never be reached */ +} +/* }}} */ + +/***************************************************************** + Cumulative Distribution Function + CHI-Square distribution + + Function + + Calculates any one parameter of the chi-square + distribution given values for the others. + + Arguments + + + WHICH --> Integer indicating which of the next three argument + values is to be calculated from the others. + Legal range: 1..3 + iwhich = 1 : Calculate P and Q from X and DF + iwhich = 2 : Calculate X from P,Q and DF + iwhich = 3 : Calculate DF from P,Q and X + + P <--> The integral from 0 to X of the chi-square + distribution. + Input range: [0, 1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + X <--> Upper limit of integration of the non-central + chi-square distribution. + Input range: [0, +infinity). + Search range: [0,1E100] + + DF <--> Degrees of freedom of the + chi-square distribution. + Input range: (0, +infinity). + Search range: [ 1E-100, 1E100] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + 10 indicates error returned from cumgam. See + references in cdfgam + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + + Method + + + Formula 26.4.19 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to reduce the chisqure + distribution to the incomplete distribution. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. +*****************************************************************/ +/* {{{ proto float stats_cdf_chisquare(float par1, float par2, float par3, int which) + Calculates any one parameter of the chi-square distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_chisquare) +{ + double arg1; + double arg2; + double bound; + double p; + double q; + double x; + double df; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddl", &arg1, &arg2, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 3) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); + RETURN_FALSE; + } + + if (which < 3 ) { + df = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfchi((int *)&which, &p, &q, &x, &df, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(df); + } + RETURN_FALSE; /* should never be here */ +} +/* }}} */ + +/******************************************************************* + Cumulative Distribution Function + BETa Distribution + + Function + + Calculates any one parameter of the beta distribution given + values for the others. + + Arguments + + WHICH --> Integer indicating which of the next four argument + values is to be calculated from the others. + Legal range: 1..4 + iwhich = 1 : Calculate P and Q from X,Y,A and B + iwhich = 2 : Calculate X and Y from P,Q,A and B + iwhich = 3 : Calculate A from P,Q,X,Y and B + iwhich = 4 : Calculate B from P,Q,X,Y and A + + P <--> The integral from 0 to X of the chi-square + distribution. + Input range: [0, 1]. + + Q <--> 1-P. + Input range: [0, 1]. + P + Q = 1.0. + + X <--> Upper limit of integration of beta density. + Input range: [0,1]. + Search range: [0,1] + + Y <--> 1-X. + Input range: [0,1]. + Search range: [0,1] + X + Y = 1.0. + + A <--> The first parameter of the beta density. + Input range: (0, +infinity). + Search range: [1D-100,1D100] + + B <--> The second parameter of the beta density. + Input range: (0, +infinity). + Search range: [1D-100,1D100] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + 4 if X + Y .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Cumulative distribution function (P) is calculated directly by + code associated with the following reference. + + DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant + Digit Computation of the Incomplete Beta Function Ratios. ACM + Trans. Math. Softw. 18 (1993), 360-373. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + + Note + + The beta density is proportional to + t^(A-1) * (1-t)^(B-1) + +*******************************************************************/ + +/* {{{ proto float stats_cdf_beta(float par1, float par2, float par3, int which) + Calculates any one parameter of the beta distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_beta) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double x; + double bound; + double y; + double a; + double b; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + + if (which < 4) { + b = arg3; + } else { + a = arg3; + } + + if (which < 3) { + a = arg2; + } else { + x = arg2; + y = 1.0 - x; + } + + if (which == 1) { + x = arg1; + y = 1.0 - x; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfbet((int *)&which, &p, &q, &x, &y, &a, &b, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(a); + case 4: RETURN_DOUBLE(b); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/********************************************************************* + Cumulative Distribution Function + BINomial distribution + + Function + + Calculates any one parameter of the binomial + distribution given values for the others. + + Arguments + + WHICH --> Integer indicating which of the next four argument + values is to be calculated from the others. + Legal range: 1..4 + iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR + iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR + iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR + iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN + + P <--> The cumulation from 0 to S of the binomial distribution. + (Probablility of S or fewer successes in XN trials each + with probability of success PR.) + Input range: [0,1]. + + Q <--> 1-P. + Input range: [0, 1]. + P + Q = 1.0. + + S <--> The number of successes observed. + Input range: [0, XN] + Search range: [0, XN] + + XN <--> The number of binomial trials. + Input range: (0, +infinity). + Search range: [1E-100, 1E100] + + PR <--> The probability of success in each binomial trial. + Input range: [0,1]. + Search range: [0,1] + + OMPR <--> 1-PR + Input range: [0,1]. + Search range: [0,1] + PR + OMPR = 1.0 + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + 4 if PR + OMPR .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.5.24 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to reduce the binomial + distribution to the cumulative incomplete beta distribution. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + +*********************************************************************/ + +/* {{{ proto float stats_cdf_binomial(float par1, float par2, float par3, int which) + Calculates any one parameter of the binomial distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_binomial) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double xn; + double bound; + double sn; + double pr; + double ompr; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + + if (which < 4) { + pr = arg3; + ompr = 1.0 - pr; + } else { + xn = arg3; + } + + if (which < 3) { + xn = arg2; + } else { + sn = arg2; + } + + if (which == 1) { + sn = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfbin((int *)&which, &p, &q, &sn, &xn, &pr, &ompr, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in binomialcdf"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(sn); + case 3: RETURN_DOUBLE(xn); + case 4: RETURN_DOUBLE(pr); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/***************************************************************** + Cumulative Distribution Function + Non-central Chi-Square + + Function + + Calculates any one parameter of the non-central chi-square + distribution given values for the others. + + Arguments + + WHICH --> Integer indicating which of the next three argument + values is to be calculated from the others. + Input range: 1..4 + iwhich = 1 : Calculate P and Q from X and DF + iwhich = 2 : Calculate X from P,DF and PNONC + iwhich = 3 : Calculate DF from P,X and PNONC + iwhich = 3 : Calculate PNONC from P,X and DF + + P <--> The integral from 0 to X of the non-central chi-square + distribution. + Input range: [0, 1-1E-16). + + Q <--> 1-P. + Q is not used by this subroutine and is only included + for similarity with other cdf* routines. + + X <--> Upper limit of integration of the non-central + chi-square distribution. + Input range: [0, +infinity). + Search range: [0,1E100] + + DF <--> Degrees of freedom of the non-central + chi-square distribution. + Input range: (0, +infinity). + Search range: [ 1E-100, 1E100] + + PNONC <--> Non-centrality parameter of the non-central + chi-square distribution. + Input range: [0, +infinity). + Search range: [0,1E4] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.4.25 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to compute the cumulative + distribution function. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + + WARNING + + The computation time required for this routine is proportional + to the noncentrality parameter (PNONC). Very large values of + this parameter can consume immense computer resources. This is + why the search range is bounded by 10,000. + +*****************************************************************/ + +/* {{{ proto float stats_cdf_noncentral_chisquare(float par1, float par2, float par3, int which) + Calculates any one parameter of the non-central chi-square distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_noncentral_chisquare) +{ + double arg1; + double arg2; + double arg3; + double p; + double bound; + double q; + double x; + double df; + double pnonc; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + pnonc = arg3; + } else { + df = arg3; + } + + if (which < 3) { + df = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfchn((int *)&which, &p, &q, &x, &df, &pnonc, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdfchn"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(df); + case 4: RETURN_DOUBLE(pnonc); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/************************************************************** + Cumulative Distribution Function F distribution + + Function + + Calculates any one parameter of the F distribution + given values for the others. + + Arguments + + WHICH --> Integer indicating which of the next four argument + values is to be calculated from the others. + Legal range: 1..4 + iwhich = 1 : Calculate P and Q from F,DFN and DFD + iwhich = 2 : Calculate F from P,Q,DFN and DFD + iwhich = 3 : Calculate DFN from P,Q,F and DFD + iwhich = 4 : Calculate DFD from P,Q,F and DFN + + P <--> The integral from 0 to F of the f-density. + Input range: [0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + F <--> Upper limit of integration of the f-density. + Input range: [0, +infinity). + Search range: [0,1E100] + + DFN < --> Degrees of freedom of the numerator sum of squares. + Input range: (0, +infinity). + Search range: [ 1E-100, 1E100] + + DFD < --> Degrees of freedom of the denominator sum of squares. + Input range: (0, +infinity). + Search range: [ 1E-100, 1E100] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.6.2 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to reduce the computation + of the cumulative distribution function for the F variate to + that of an incomplete beta. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + + WARNING + + The value of the cumulative F distribution is not necessarily + monotone in either degrees of freedom. There thus may be two + values that provide a given CDF value. This routine assumes + monotonicity and will find an arbitrary one of the two values. + +**********************************************************************/ +/* {{{ proto float stats_cdf_f(float par1, float par2, float par3, int which) + Calculates any one parameter of the F distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_f) +{ + double arg1; + double arg2; + double arg3; + double p; + double bound; + double dfn; + double q; + double f; + double dfd; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + if (which < 4) { + dfd = arg3; + } else { + dfn = arg3; + } + if (which < 3) { + dfn = arg2; + } else { + f = arg2; + } + if (which == 1) { + f = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdff((int *)&which, &p, &q, &f, &dfn, &dfd, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdff"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(f); + case 3: RETURN_DOUBLE(dfn); + case 4: RETURN_DOUBLE(dfd); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/*********************************************************************** + Cumulative Distribution Function + Non-central F distribution + + Function + + Calculates any one parameter of the Non-central F + distribution given values for the others. + + Arguments + + WHICH --> Integer indicating which of the next five argument + values is to be calculated from the others. + Legal range: 1..5 + iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC + iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC + iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC + iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC + iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD + + P <--> The integral from 0 to F of the non-central f-density. + Input range: [0,1-1E-16). + + Q <--> 1-P. + Q is not used by this subroutine and is only included + for similarity with other cdf* routines. + + F <--> Upper limit of integration of the non-central f-density. + Input range: [0, +infinity). + Search range: [0,1E100] + + DFN < --> Degrees of freedom of the numerator sum of squares. + Input range: (0, +infinity). + Search range: [ 1E-100, 1E100] + + DFD < --> Degrees of freedom of the denominator sum of squares. + Must be in range: (0, +infinity). + Input range: (0, +infinity). + Search range: [ 1E-100, 1E100] + + PNONC <-> The non-centrality parameter + Input range: [0,infinity) + Search range: [0,1E4] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.6.20 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to compute the cumulative + distribution function. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + + WARNING + + The computation time required for this routine is proportional + to the noncentrality parameter (PNONC). Very large values of + this parameter can consume immense computer resources. This is + why the search range is bounded by 10,000. + + WARNING + + The value of the cumulative noncentral F distribution is not + necessarily monotone in either degrees of freedom. There thus + may be two values that provide a given CDF value. This routine + assumes monotonicity and will find an arbitrary one of the two + values. + +***********************************************************************/ + +/* {{{ proto float stats_cdf_noncentral_f(float par1, float par2, float par3, float par4, int which) + Calculates any one parameter of the Non-central F distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_noncentral_f) +{ + double arg1; + double arg2; + double arg3; + double arg4; + double p; + double q; + double f; + double dfn; + double dfd; + double pnonc; + double bound; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddddl", &arg1, &arg2, &arg3, &arg4, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 5) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fifth parameter should be in the 1..5 range"); + RETURN_FALSE; + } + + if (which < 5) { + pnonc = arg4; + } else { + dfd = arg4; + } + + if (which < 4) { + dfd = arg3; + } else { + dfn = arg3; + } + + if (which < 3) { + dfn = arg2; + } else { + f = arg2; + } + + if (which == 1) { + f = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdffnc((int *)&which, &p, &q, &f, &dfn, &dfd, &pnonc, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdffnc"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(f); + case 3: RETURN_DOUBLE(dfn); + case 4: RETURN_DOUBLE(dfd); + case 5: RETURN_DOUBLE(pnonc); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/************************************************************************ + Cumulative Distribution Function Non-Central T distribution + + Function + + Calculates any one parameter of the noncentral t distribution give + values for the others. + + Arguments + + WHICH --> Integer indicating which argument + values is to be calculated from the others. + Legal range: 1..3 + iwhich = 1 : Calculate P and Q from T,DF,PNONC + iwhich = 2 : Calculate T from P,Q,DF,PNONC + iwhich = 3 : Calculate DF from P,Q,T + iwhich = 4 : Calculate PNONC from P,Q,DF,T + + P <--> The integral from -infinity to t of the noncentral t-den + Input range: (0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + T <--> Upper limit of integration of the noncentral t-density. + Input range: ( -infinity, +infinity). + Search range: [ -1E100, 1E100 ] + + DF <--> Degrees of freedom of the noncentral t-distribution. + Input range: (0 , +infinity). + Search range: [1e-100, 1E10] + + PNONC <--> Noncentrality parameter of the noncentral t-distributio + Input range: [-infinity , +infinity). + Search range: [-1e4, 1E4] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Upper tail of the cumulative noncentral t is calculated usin + formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuou + Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995) + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + +************************************************************************/ + +/* {{{ proto float stats_stat_noncentral_t(float par1, float par2, float par3, int which) + Calculates any one parameter of the noncentral t distribution give values for the others. */ +PHP_FUNCTION(stats_cdf_noncentral_t) +{ + double arg1; + double arg2; + double arg3; + double pnonc; + double bound; + double p; + double q; + double t; + double df; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + pnonc = arg3; + } else { + df = arg3; + } + if (which < 3) { + df = arg2; + } else { + t = arg2; + } + + if (which == 1) { + t = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdftnc((int *)&which, &p, &q, &t, &df, &pnonc, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(t); + case 3: RETURN_DOUBLE(df); + case 4: RETURN_DOUBLE(pnonc); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/************************************************************************ + Cumulative Distribution Function Negative BiNomial distribution + + Function + + Calculates any one parameter of the negative binomial + distribution given values for the others. + + The cumulative negative binomial distribution returns the + probability that there will be F or fewer failures before the + XNth success in binomial trials each of which has probability of + success PR. + + The individual term of the negative binomial is the probability of + S failures before XN successes and is + Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S + + Arguments + + WHICH --> Integer indicating which of the next four argument + values is to be calculated from the others. + Legal range: 1..4 + iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR + iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR + iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR + iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN + + P <--> The cumulation from 0 to S of the negative + binomial distribution. + Input range: [0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + S <--> The upper limit of cumulation of the binomial distribution. + There are F or fewer failures before the XNth success. + Input range: [0, +infinity). + Search range: [0, 1E100] + + XN <--> The number of successes. + Input range: [0, +infinity). + Search range: [0, 1E100] + + PR <--> The probability of success in each binomial trial. + Input range: [0,1]. + Search range: [0,1]. + + OMPR <--> 1-PR + Input range: [0,1]. + Search range: [0,1] + PR + OMPR = 1.0 + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + 4 if PR + OMPR .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.5.26 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to reduce calculation of + the cumulative distribution function to that of an incomplete + beta. + + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + +************************************************************************/ + +/* {{{ proto float stats_cdf_negative_binomial(float par1, float par2, float par3, int which) + Calculates any one parameter of the negative binomial distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_negative_binomial) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double bound; + double sn; + double xn; + double pr; + double ompr; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + pr = arg3; + ompr = 1.0 - pr; + } else { + xn = arg3; + } + + if (which < 3) { + xn = arg2; + } else { + sn = arg2; + } + + if (which == 1) { + sn = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfnbn((int *)&which, &p, &q, &sn, &xn, &pr, &ompr, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdfnbn"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(sn); + case 3: RETURN_DOUBLE(xn); + case 4: RETURN_DOUBLE(pr); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/********************************************************************** + Cumulative Distribution Function POIsson distribution + + Function + + Calculates any one parameter of the Poisson + distribution given values for the others. + + Arguments + + WHICH --> Integer indicating which argument + value is to be calculated from the others. + Legal range: 1..3 + iwhich = 1 : Calculate P and Q from S and XLAM + iwhich = 2 : Calculate A from P,Q and XLAM + iwhich = 3 : Calculate XLAM from P,Q and S + + P <--> The cumulation from 0 to S of the poisson density. + Input range: [0,1]. + + Q <--> 1-P. + Input range: (0, 1]. + P + Q = 1.0. + + S <--> Upper limit of cumulation of the Poisson. + Input range: [0, +infinity). + Search range: [0,1E100] + + XLAM <--> Mean of the Poisson distribution. + Input range: [0, +infinity). + Search range: [0,1E100] + + STATUS <-- 0 if calculation completed correctly + -I if input parameter number I is out of range + 1 if answer appears to be lower than lowest + search bound + 2 if answer appears to be higher than greatest + search bound + 3 if P + Q .ne. 1 + + BOUND <-- Undefined if STATUS is 0 + + Bound exceeded by parameter number I if STATUS + is negative. + + Lower search bound if STATUS is 1. + + Upper search bound if STATUS is 2. + + Method + + Formula 26.4.21 of Abramowitz and Stegun, Handbook of + Mathematical Functions (1966) is used to reduce the computation + of the cumulative distribution function to that of computing a + chi-square, hence an incomplete gamma function. + + Cumulative distribution function (P) is calculated directly. + Computation of other parameters involve a seach for a value that + produces the desired value of P. The search relies on the + monotinicity of P with the other parameter. + +**********************************************************************/ + +/* {{{ proto float stats_cdf_poisson(float par1, float par2, float par3, int which) + Calculates any one parameter of the Poisson distribution given values for the others. */ +PHP_FUNCTION(stats_cdf_poisson) +{ + double arg1; + double arg2; + double p; + double q; + double x; + double xlam; + double bound; + long which; + int status = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddl", &arg1, &arg2, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 3) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); + RETURN_FALSE; + } + + if (which < 3 ) { + xlam = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + cdfpoi((int *)&which, &p, &q, &x, &xlam, &status, &bound); + if (status != 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); + RETURN_FALSE; + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(xlam); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + + +static double laplace_quantile(double p) +{ + if (p <= 0.5) { + return log(2.0*p); + } else { + return (-log(2.0*(1.0-p))); + } +} + +static double laplace_cdf(double x) +{ + if (x <= 0) { + return (0.5*exp(x)); + } else { + return (1.0 - 0.5*exp(-x)); + } +} + + +/* {{{ proto float stats_cdf_laplace(float par1, float par2, float par3, int which) + Not documented */ +PHP_FUNCTION(stats_cdf_laplace) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double x; + double t; + double mean; + double sd; + long which; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + sd = arg3; + } else { + mean = arg3; + } + + if (which < 3) { + mean = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + if (which == 1) { + t = (x - mean) / sd; + p = laplace_cdf(t); + } else { + t = laplace_quantile(p); + } + + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(mean + (sd * t)); + case 3: RETURN_DOUBLE(x - (sd * t)); + case 4: RETURN_DOUBLE((x - mean) / t); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +static double cauchy_quantile(double p) +{ + return (tan(STATS_PI*(p-0.5))); +} + +static double cauchy_cdf (double x) +{ + return (0.5+(atan(x)/STATS_PI)); +} + +/* {{{ proto float stats_cdf_cauchy(float par1, float par2, float par3, int which) + Not documented */ +PHP_FUNCTION(stats_cdf_cauchy) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double x; + double t; + double mean; + double sd; + long which; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + sd = arg3; + } else { + mean = arg3; + } + + if (which < 3) { + mean = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + if (which == 1) { + t = (x - mean) / sd; + p = cauchy_cdf(t); + } else { + t = cauchy_quantile(p); + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(mean + (sd * t)); + case 3: RETURN_DOUBLE(x - (sd * t)); + case 4: RETURN_DOUBLE((x - mean) / t); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +static double logistic_cdf(double x) +{ + return (1.0/(1.0+exp(-x))); +} + +static double logistic_quantile (double p) +{ + return (log(p/(1.0-p))); +} + +/* {{{ proto float stats_cdf_logistic(float par1, float par2, float par3, int which) + Not documented */ +PHP_FUNCTION(stats_cdf_logistic) +{ + double arg1; + double arg2; + double arg3; + double sd; + double p; + double q; + double x; + double t; + double mean; + long which; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + sd = arg3; + } else { + mean = arg3; + } + + if (which < 3) { + mean = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + if (which == 1) { + t = (x - mean) / sd; + p = logistic_cdf(t); + } else { + t = logistic_quantile(p); + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(mean + (sd * t)); + case 3: RETURN_DOUBLE(x - (sd * t)); + case 4: RETURN_DOUBLE((x - mean) / t); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/* {{{ proto float stats_cdf_weibull(float par1, float par2, float par3, int which) + Not documented */ +PHP_FUNCTION(stats_cdf_weibull) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double x; + double a; + double b; + long which; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + b = arg3; + } else { + a = arg3; + } + + if (which < 3) { + a = arg2; + } else { + x = arg2; + } + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + if (which == 1) { + p = 1 - exp(-pow(x / b, a)); + } else { + x = b * pow(-log(1.0 - p), 1.0 / a); + } + + switch (which) { + case 1: RETURN_DOUBLE(p); + case 2: RETURN_DOUBLE(x); + case 3: RETURN_DOUBLE(log(-log(1.0 - p)) / log(x / b)); + case 4: RETURN_DOUBLE(x / pow(-log(1.0 - p), 1.0 / a)); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +/* {{{ proto float stats_cdf_uniform(float par1, float par2, float par3, int which) + Not documented */ +PHP_FUNCTION(stats_cdf_uniform) +{ + double arg1; + double arg2; + double arg3; + double p; + double q; + double x; + double a; + double b; + long which; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 4) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); + RETURN_FALSE; + } + + if (which < 4) { + b = arg3; + } else { + a = arg3; + } + + if (which < 3) { + a = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + if (which == 1) { + p = 1 - exp(-pow(x / b, a)); + } else { + x = b * pow(-log(1.0 - p), 1.0 / a); + } + + switch (which) { + case 4: RETURN_DOUBLE((x - (1.0 - p) * a) / p); + case 3: RETURN_DOUBLE((x - p * b) / (1.0 - p)); + case 2: RETURN_DOUBLE(a + p * (b - a)); + case 1: + if (x < a) { + p = 0; + } else { + if (x > b) { + p = 1; + } else { + p = (x - a) / ( b - a); + } + } + RETURN_DOUBLE(p); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + +static double exponential_quantile(double p) +{ + return -log(1.0-p); +} + +static double exponential_cdf(double x) +{ + return (1.0 - exp(-x)); +} + +/* {{{ proto float stats_cdf_exponential(float par1, float par2, int which) + Not documented */ +PHP_FUNCTION(stats_cdf_exponential) +{ + double arg1; + double arg2; + double p; + double q; + double x; + double scale; + long which; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddl", &arg1, &arg2, &which) == FAILURE) { + RETURN_FALSE; + } + + if (which < 1 || which > 3) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); + RETURN_FALSE; + } + + if (which < 3) { + scale = arg2; + } else { + x = arg2; + } + + if (which == 1) { + x = arg1; + } else { + p = arg1; + q = 1.0 - p; + } + + switch (which) { + case 1: RETURN_DOUBLE(exponential_cdf(x / scale)); + case 2: RETURN_DOUBLE(scale * exponential_quantile(p)); + case 3: RETURN_DOUBLE(x / exponential_quantile(p)); + } + RETURN_FALSE; /* never here */ +} +/* }}} */ + + +/*********************/ +/* RANDLIB functions */ +/*********************/ + +/* {{{ proto void stats_rand_setall(int iseed1, int iseed2) + Not documented */ +PHP_FUNCTION(stats_rand_setall) +{ + long iseed_1; + long iseed_2; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ll", &iseed_1, &iseed_2) == FAILURE) { + RETURN_FALSE; + } + + setall(iseed_1, iseed_2); +} +/* }}} */ + +/* {{{ proto array stats_rand_get_seeds(void) + Not documented */ +PHP_FUNCTION(stats_rand_getsd) +{ + long iseed_1; + long iseed_2; + + if (ZEND_NUM_ARGS() != 0) { + WRONG_PARAM_COUNT; + } + getsd(&iseed_1, &iseed_2); + + array_init(return_value); + add_next_index_long(return_value, iseed_1); + add_next_index_long(return_value, iseed_2); +} +/* }}} */ + +/* {{{ proto int stats_rand_gen_iuniform(int low, int high) + Generates integer uniformly distributed between LOW (inclusive) and HIGH (inclusive) */ +PHP_FUNCTION(stats_rand_gen_iuniform) +{ + long low; + long high; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ll", &low, &high) == FAILURE) { + RETURN_FALSE; + } + + if (high - low > 2147483561L) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "high - low too large. low : %16ld high %16ld", low, high); + RETURN_FALSE; + } + if (low > high) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "low greater than high. low : %16ld high %16ld", low, high); + RETURN_FALSE; + } + + RETURN_LONG(ignuin(low, high)); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_funiform(float low, float high) + Generates uniform float between low (exclusive) and high (exclusive) */ +PHP_FUNCTION(stats_rand_gen_funiform) +{ + double low; + double high; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &low, &high) == FAILURE) { + RETURN_FALSE; + } + + if (low > high) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "low greater than high. low : %16.6E high : %16.6E", low, high); + RETURN_FALSE; + } + + RETURN_DOUBLE(genunf(low, high)); +} +/* }}} */ + +/* {{{ proto int stats_rand_gen_int(void) + Generates random integer between 1 and 2147483562 */ +PHP_FUNCTION(stats_rand_ignlgi) +{ + if (ZEND_NUM_ARGS() != 0) { + WRONG_PARAM_COUNT; + } + + RETURN_LONG(ignlgi()); +} +/* }}} */ + +/* {{{ proto float stats_rand_ranf(void) + Returns a random floating point number from a uniform distribution over 0 - 1 (endpoints of this interval are not returned) using the current generator */ +PHP_FUNCTION(stats_rand_ranf) +{ + if (ZEND_NUM_ARGS() != 0) { + WRONG_PARAM_COUNT; + } + + RETURN_DOUBLE(ranf()); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_beta(float a, float b) + Generates beta random deviate. Returns a random deviate from the beta distribution with parameters A and B. The density of the beta is x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 0), r - shape parameter of Gamma distribution (r > 0) */ +PHP_FUNCTION(stats_rand_gen_gamma) +{ + double a; + double r; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &a, &r) == FAILURE) { + RETURN_FALSE; + } + + if (!(a > 0.0 && r > 0.0)) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "A or R nonpositive. A value : %16.6E , R value : %16.6E", a, r); + RETURN_FALSE; + } + + RETURN_DOUBLE(gengam(a, r)); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_noncenral_chisquare(float df, float xnonc) + Generates random deviate from the distribution of a noncentral chisquare with "df" degrees of freedom and noncentrality parameter "xnonc". d must be >= 1.0, xnonc must >= 0.0 */ +PHP_FUNCTION(stats_rand_gen_noncentral_chisquare) +{ + double df; + double xnonc; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &df, &xnonc) == FAILURE) { + RETURN_FALSE; + } + + if (df < 1.0 || xnonc < 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "df < 1 or xnonc < 0. df value : %16.6E xnonc value : %16.6E", df, xnonc); + RETURN_FALSE; + } + + RETURN_DOUBLE(gennch(df, xnonc)); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_noncentral_f(float dfn, float dfd, float xnonc) + Generates a random deviate from the noncentral F (variance ratio) distribution with "dfn" degrees of freedom in the numerator, and "dfd" degrees of freedom in the denominator, and noncentrality parameter "xnonc". Method : directly generates ratio of noncentral numerator chisquare variate to central denominator chisquare variate. */ +PHP_FUNCTION(stats_rand_gen_noncenral_f) +{ + double dfn; + double dfd; + double xnonc; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &dfn, &dfd, &xnonc) == FAILURE) { + RETURN_FALSE; + } + + if (dfn < 1.0 || dfd <= 0.0 || xnonc < 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Either (1) Numerator df < 1.0 or (2) Denominator df <= 0.0 or (3) Noncentrality parameter < 0.0. dfn: %16.6E dfd: %16.6E xnonc: %16.6E", dfn, dfd, xnonc); + RETURN_FALSE; + } + + RETURN_DOUBLE(gennf(dfn, dfd, xnonc)); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_normal(float av, float sd) + Generates a single random deviate from a normal distribution with mean, av, and standard deviation, sd (sd >= 0). Method : Renames SNORM from TOMS as slightly modified by BWB to use RANF instead of SUNIF. */ +PHP_FUNCTION(stats_rand_gen_normal) +{ + double av; + double sd; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &av, &sd) == FAILURE) { + RETURN_FALSE; + } + + if (sd < 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "sd < 0.0 . sd : %16.6E", sd); + RETURN_FALSE; + } + + RETURN_DOUBLE(gennor(av, sd)); +} +/* }}} */ + +/* {{{ proto array stats_rand_phrase_to_seeds(string phrase) + Uses a phrase (characted string) to generate two seeds for the RGN random number generator. Trailing blanks are eliminated before the seeds are generated. Generated seed values will fall in the range 1..2^30. */ +PHP_FUNCTION(stats_rand_phrase_to_seeds) +{ + zval **par1; + char *arg1 = NULL; + long seed_1; + long seed_2; + + if (ZEND_NUM_ARGS() != 1 || zend_get_parameters_ex(1, &par1) == FAILURE) { + WRONG_PARAM_COUNT; + } + convert_to_string_ex(par1); + + arg1 = estrndup(Z_STRVAL_PP(par1), Z_STRLEN_PP(par1)); + phrtsd(arg1, &seed_1, &seed_2); + efree(arg1); + + array_init(return_value); + add_next_index_long(return_value, seed_1); + add_next_index_long(return_value, seed_2); +} +/* }}} */ + +/* {{{ proto int stats_rand_gen_ibinomial(int n, float pp) + Generates a single random deviate from a binomial distribution whose number of trials is "n" (n >= 0) and whose probability of an event in each trial is "pp" ([0;1]). Method : algorithm BTPE */ +PHP_FUNCTION(stats_rand_ibinomial) +{ + long n; + double pp; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ld", &n, &pp) == FAILURE) { + RETURN_FALSE; + } + + if ((n < 0) || (pp < 0.0) || (pp > 1.0)) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Bad values for the arguments. n : %ld pp : %16.6E", n, pp); + RETURN_FALSE; + } + + RETURN_LONG(ignbin(n, pp)); +} +/* }}} */ + +/* {{{ proto int stats_rand_gen_ibinomial_negative(int n, float p) + Generates a single random deviate from a negative binomial distribution. Arguments : n - the number of trials in the negative binomial distribution from which a random deviate is to be generated (n > 0), p - the probability of an event (0 < p < 1)). */ +PHP_FUNCTION(stats_rand_ibinomial_negative) +{ + long n; + double p; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ld", &n, &p) == FAILURE) { + RETURN_FALSE; + } + + if (n <= 0L) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "n < 0. n : %ld", n); + RETURN_FALSE; + } + if (p < 0.0F || p > 1.0F) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "p is out of range. p : %16.E", p); + RETURN_FALSE; + } + + RETURN_LONG(ignnbn(n, p)); +} +/* }}} */ + +/* {{{ proto int stats_rand_gen_ipoisson(float mu) + Generates a single random deviate from a Poisson distribution with mean "mu" (mu >= 0.0). */ +PHP_FUNCTION(stats_rand_gen_ipoisson) +{ + double mu; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "d", &mu) == FAILURE) { + RETURN_FALSE; + } + + if (mu < 0.0F) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "mu < 0.0 . mu : %16.6E", mu); + RETURN_FALSE; + } + + RETURN_LONG(ignpoi(mu)); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_noncentral_t(float df, float xnonc) + Generates a single random deviate from a noncentral T distribution. xnonc - noncentrality parameter. df must be >= 0.0*/ +PHP_FUNCTION(stats_rand_gen_noncentral_t) +{ + double df; + double xnonc; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &df, &xnonc) == FAILURE) { + RETURN_FALSE; + } + + if (df < 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "df <= 0 . df : %16.6E", df); + RETURN_FALSE; + } + + RETURN_DOUBLE(gennor(xnonc, 1) / sqrt(genchi(df) / df) ); +} +/* }}} */ + +/* {{{ proto float stats_rand_gen_t(float df) + Generates a single random deviate from a T distribution. df must be >= 0.0 */ +PHP_FUNCTION(stats_rand_gen_t) +{ + zval **arg1; + double df; + + if (ZEND_NUM_ARGS() != 1 || zend_get_parameters_ex(1, &arg1) == FAILURE) { + WRONG_PARAM_COUNT; + } + + convert_to_double_ex(arg1); + df = Z_DVAL_PP(arg1); + + if (df < 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "df <= 0 . df : %16.6E", df); + RETURN_FALSE; + } + + RETURN_DOUBLE(gennor(0, 1) / sqrt(genchi(df) / df)); +} +/* }}} */ + +/***************************/ +/* Start density functions */ +/***************************/ + +/* {{{ proto float stats_dens_normal(float x, float ave, float stdev) + Not documented */ +PHP_FUNCTION(stats_dens_normal) +{ + double stdev; + double ave; + double x; + double y; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddd", &x, &ave, &stdev) == FAILURE) { + RETURN_FALSE; + } + + if (stdev == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); + RETURN_FALSE; + } + + z = (x - ave) / stdev; + y = (1.0 / (stdev * sqrt(2.0 * STATS_PI))) * exp (-0.5 * z * z); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_cauchy(float x, float ave, float stdev) + Not documented */ +PHP_FUNCTION(stats_dens_cauchy) +{ + double stdev; + double ave; + double x; + double y; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddd", &x, &ave, &stdev) == FAILURE) { + RETURN_FALSE; + } + + if (stdev == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); + RETURN_FALSE; + } + + z = (x - ave) / stdev; + y = 1.0 / (stdev*STATS_PI * (1.0 + (z * z))); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_laplace(float x, float ave, float stdev) + Not documented */ +PHP_FUNCTION(stats_dens_laplace) +{ + double stdev; + double ave; + double x; + double y; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddd", &x, &ave, &stdev) == FAILURE) { + RETURN_FALSE; + } + + if (stdev == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); + RETURN_FALSE; + } + + z = fabs((x - ave) / stdev); + y = (1.0 / (2.0 * stdev)) * exp(- z); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_logistic(float x, float ave, float stdev) + Not documented */ +PHP_FUNCTION(stats_dens_logistic) +{ + double stdev; + double ave; + double x; + double y; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddd", &x, &ave, &stdev) == FAILURE) { + RETURN_FALSE; + } + + if (stdev == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); + RETURN_FALSE; + } + + z = exp((x - ave) / stdev); + y = z / (stdev * pow(1 + z, 2.0)); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_beta(float x, float a, float b) + Not documented */ +PHP_FUNCTION(stats_dens_beta) +{ + double a; + double b; + double beta; + double x; + double y; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &a, &b) == FAILURE) { + RETURN_FALSE; + } + + beta = 1.0 / exp(lgamma(a) + lgamma(b) - lgamma(a + b)); + y = beta * pow(x, a - 1.0) * pow(1.0 - x, b - 1.0); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_weibull(float x, float a, float b) + Not documented */ +PHP_FUNCTION(stats_dens_weibull) +{ + double a; + double b; + double x; + double y; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &a, &b) == FAILURE) { + RETURN_FALSE; + } + + if (b == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "b is 0.0"); + RETURN_FALSE; + } + + y = (a / b) * pow(x / b, a - 1.0) * exp(pow(- x / b, a)); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_uniform(float x, float a, float b) + Not documented */ +PHP_FUNCTION(stats_dens_uniform) +{ + double a; + double b; + double x; + double y; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddd", &x, &a, &b) == FAILURE) { + RETURN_FALSE; + } + + if (a == b) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "b == a == %16.6E", a); + RETURN_FALSE; + } + + if ((x <= b) && (x >= a)) { + y = 1.0 / (b - a); + } else { + y = 0.0; + } + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_chisquare(float x, float dfr) + Not documented */ +PHP_FUNCTION(stats_dens_chisquare) +{ + double dfr; + double e; + double x; + double y; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "dd", &x, &dfr) == FAILURE) { + RETURN_FALSE; + } + + e = dfr / 2.0; + z = ((e - 1.0) * log(x)) - ((x / 2.0) +(e * log(2.0)) + lgamma(e)); + y = exp (z); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_t(float x, float dfr) + Not documented */ +PHP_FUNCTION(stats_dens_t) +{ + double dfr; + double e; + double f; + double fac1; + double fac2; + double fac3; + double x; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &x, &dfr) == FAILURE) { + RETURN_FALSE; + } + + if (dfr == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "dfr == 0.0"); + RETURN_FALSE; + } + + e = dfr / 2.0; + f = e + 0.5; + fac1 = lgamma(f); + fac2 = f * log(1.0 + (x * x) / dfr); + fac3 = lgamma(e) + 0.5 * log(dfr * STATS_PI); + + RETURN_DOUBLE(exp( fac1 - (fac2 + fac3) )); +} +/* }}} */ + +/* {{{ proto float stats_dens_gamma(float x, float shape, float scale) + Not documented */ +PHP_FUNCTION(stats_dens_gamma) +{ + double shape; + double scale; + double x; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &shape, &scale) == FAILURE) { + RETURN_FALSE; + } + + if (scale == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "scale == 0.0"); + RETURN_FALSE; + } + + z = ((shape - 1.0) * log(x)) - + ( + (x / scale) + lgamma(shape) + (shape * log(scale)) + ) + ; + + RETURN_DOUBLE(exp(z)); +} +/* }}} */ + +/* {{{ proto float stats_dens_exponential(float x, float scale) + Not documented */ +PHP_FUNCTION(stats_dens_exponential) +{ + double scale; + double x; + double y; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &x, &scale) == FAILURE) { + RETURN_FALSE; + } + + if (scale == 0.0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "scale == 0.0"); + RETURN_FALSE; + } + + if (x < 0) { + y = 0; + } else { + y = exp(-x / scale) / scale; + } + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_f(float x, float dfr1, float dfr2) + */ +PHP_FUNCTION(stats_dens_f) +{ + double dfr1; + double dfr2; + double efr1; + double efr2; + double fac1; + double fac2; + double fac3; + double fac4; + double x; + double z; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &dfr1, &dfr2) == FAILURE) { + RETURN_FALSE; + } + + efr1 = dfr1 / 2.0; + efr2 = dfr2 / 2.0; + fac1 = (efr1 - 1.0) * log (x); + fac2 = (efr1 + efr2) * log (dfr2 + (dfr1 * x)); + fac3 = (efr1 * log (dfr1)) + (efr2 * log (dfr2)); + fac4 = lgamma (efr1) + lgamma (efr2) - lgamma (efr1 + efr2); + + z = (fac1 + fac3) - (fac2 + fac4); + + RETURN_DOUBLE(exp(z)); +} +/* }}} */ + +static double binom(double x, double n) +{ + int i; + double di; + double s = 1.0; + + for (i = 0; i < x; ++i) { + di = (double) i; + s = (s * (n - di)) / (di + 1.0); + } + + return s; +} + +/* {{{ proto float stats_dens_pmf_binomial(float x, float n, float pi) + Not documented */ +PHP_FUNCTION(stats_dens_pmf_binomial) +{ + double pi; + double y; + double n; + double x; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, + "ddd", &x, &n, &pi) == FAILURE) { + RETURN_FALSE; + } + + if ((x == 0.0 && n == 0.0) || (pi == 0.0 && x == 0.0) + || ( (1.0 - pi) == 0.0 && (n - x) == 0) ) { + + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Params leading to pow(0, 0). x:%16.6E n:%16.6E pi:%16.6E", x, n, pi); + RETURN_FALSE; + } + + y = binom(x,n) * pow(pi,x) * pow((1.0 - pi), (n - x)); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_pmf_poisson(float x, float lb) + Not documented */ +PHP_FUNCTION(stats_dens_pmf_poisson) +{ + double lb; + double z; + double x; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &x, &lb) == FAILURE) { + RETURN_FALSE; + } + + z = (x * log(lb)) - (lb + lgamma(x + 1.0)); + + RETURN_DOUBLE(exp(z)); +} +/* }}} */ + +/* {{{ proto float stats_dens_negative_binomial(float x, float n, float pi) + Not documented */ +PHP_FUNCTION(stats_dens_pmf_negative_binomial) +{ + double pi; + double y; + double n; + double x; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &n, &pi) == FAILURE) { + RETURN_FALSE; + } + + if ((pi == 0.0 && n == 0.0) || ((1.0 - pi) == 0.0 && x == 0.0)) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Params leading to pow(0, 0). x:%16.6E n:%16.6E pi:%16.6E", x, n, pi); + RETURN_FALSE; + } + + y = binom(x, n + x - 1.0) * pow(pi,n) * pow((1.0 - pi), x); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/* {{{ proto float stats_dens_pmf_hypergeometric(float n1, float n2, float N1, float N2) + */ +PHP_FUNCTION(stats_dens_pmf_hypergeometric) +{ + double y; + double N1; + double N2; + double n1; + double n2; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dddd", &n1, &n2, &N1, &N2) == FAILURE) { + RETURN_FALSE; + } + + if ((int)(n1 + n2) >= (int)(N1 + N2)) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "possible division by zero - n1+n2 >= N1+N2"); + /* RETURN_FALSE; */ + } + + y = binom(n1, N1) * binom (n2, N2)/binom(n1 + n2, N1 + N2); + + RETURN_DOUBLE(y); +} +/* }}} */ + +/************************/ +/* Statistics functions */ +/************************/ + +/* {{{ proto float stats_stat_powersum(array arr, float power) + Not documented */ +PHP_FUNCTION(stats_stat_powersum) +{ + zval **arg1, **arg2, **data; /* pointer to array entry */ + HashPosition pos; /* hash iterator */ + double power; + double sum = 0.0; + + if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { + WRONG_PARAM_COUNT; + } + + convert_to_array_ex(arg1); + convert_to_double_ex(arg2); + power = Z_DVAL_PP(arg2); + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data, &pos) == SUCCESS) { + convert_to_double_ex(data); + if (Z_DVAL_PP(data) != 0 && power != 0) { + sum += pow (Z_DVAL_PP(data), power); + } else { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Both value and power are zero"); + } + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos); + } + + RETURN_DOUBLE(sum); +} +/* }}} */ + +/* {{{ proto float stats_stat_innerproduct(array arr1, array arr2) + */ +PHP_FUNCTION(stats_stat_innerproduct) +{ + zval **arg1, **arg2; + zval **data1, **data2; /* pointers to array entries */ + HashPosition pos1; /* hash iterator */ + HashPosition pos2; /* hash iterator */ + double sum = 0.0; + + + if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { + WRONG_PARAM_COUNT; + } + convert_to_array_ex(arg1); + convert_to_array_ex(arg2); + + if (zend_hash_num_elements(Z_ARRVAL_PP(arg1)) != zend_hash_num_elements(Z_ARRVAL_PP(arg2))) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Unequal number of X and Y coordinates"); + RETURN_FALSE; + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS + && zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { + convert_to_double_ex(data1); + convert_to_double_ex(data2); + sum = Z_DVAL_PP(data1) * Z_DVAL_PP(data2); + + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); + } + + RETURN_DOUBLE(sum); +} +/* }}} */ + +/* {{{ proto float stats_stat_independent_t(array arr1, array arr2) + Not documented */ +PHP_FUNCTION(stats_stat_independent_t) +{ + zval **arg1, **arg2; + zval **data1, **data2; /* pointers to array entries */ + HashPosition pos1; /* hash iterator */ + HashPosition pos2; /* hash iterator */ + int xnum = 0, ynum = 0; + double cur; + double sx = 0.0; + double sxx = 0.0; + double sy = 0.0; + double syy = 0.0; + double mx; + double vx; + double my; + double vy; + double sp; + double fc; + double ts; + + if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { + WRONG_PARAM_COUNT; + } + convert_to_array_ex(arg1); + convert_to_array_ex(arg2); + + xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); + ynum = zend_hash_num_elements(Z_ARRVAL_PP(arg2)); + if ( xnum < 2 || ynum < 2) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Each argument should have more than 1 element"); + RETURN_FALSE; + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS) { + convert_to_double_ex(data1); + cur = Z_DVAL_PP(data1); + sx += cur; + sxx += cur * cur; + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { + convert_to_double_ex(data2); + cur = Z_DVAL_PP(data2); + sy += cur; + syy += cur * cur; + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); + } + + mx = sx / xnum; + my = sy / ynum; + vx = (sxx - (xnum * mx * mx)) / (xnum - 1.0); + vy = (syy - (ynum * my * my)) / (ynum - 1.0); + sp = (((xnum - 1.0) * vx) + ((ynum - 1.0) * vy)) / (xnum + ynum - 2.0); + fc = (1.0 / xnum) + (1.0 / ynum); + ts = (mx - my) / sqrt(sp * fc); + + RETURN_DOUBLE(ts); +} +/* }}} */ + +/* {{{ proto float stats_stat_paired_t(array arr1, array arr2) + Not documented */ +PHP_FUNCTION(stats_stat_paired_t) +{ + zval **arg1, **arg2, **data1, **data2; /* pointers to array entries */ + HashPosition pos1; /* hash iterator */ + HashPosition pos2; /* hash iterator */ + int xnum = 0; + int ynum = 0; + double sd = 0.0; + double sdd = 0.0; + double md; + double td; + double ts; + double cur; + + if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { + WRONG_PARAM_COUNT; + } + convert_to_array_ex(arg1); + convert_to_array_ex(arg2); + + xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); + ynum = zend_hash_num_elements(Z_ARRVAL_PP(arg2)); + + if (xnum != ynum) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Unequal number of X and Y coordinates"); + RETURN_FALSE; + } + if (xnum < 2) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "arr1 should have atleast 2 elements"); + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS + && + zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { + convert_to_double_ex(data1); + convert_to_double_ex(data2); + + cur = Z_DVAL_PP(data1) - Z_DVAL_PP(data2); + sd += cur; + sdd += cur * cur; + + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); + } + + md = sd / xnum; + td = sqrt((sdd - (xnum * md * md)) / (xnum - 1.0)); + ts = sqrt((double) xnum) * (md / td); + + RETURN_DOUBLE(ts); +} +/* }}} */ + +/* {{{ proto float stats_stat_percentile(float df, float xnonc) + Not documented */ +PHP_FUNCTION(stats_stat_percentile) +{ + zval **arg1, **arg2; + zval **data1; /* pointers to array entries */ + HashPosition pos1; /* hash iterator */ + long ilow; + long iupp; + int xnum = 0; + int cnt = -1; + double perc; + double low; + double upp; + double val = 0.0; + + if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { + WRONG_PARAM_COUNT; + } + + convert_to_array_ex(arg1); + convert_to_double_ex(arg2); + perc = Z_DVAL_PP(arg2); + + xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); + + if (zend_hash_sort(Z_ARRVAL_PP(arg1), zend_qsort, stats_array_data_compare, 1 TSRMLS_CC) == FAILURE) { + RETURN_FALSE; + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); + + low = .01 * perc * (double)xnum; + upp = .01 * (100.0 - perc) * (double)xnum; + ilow = floor(low); + iupp = floor(upp); + if ((ilow + iupp) == xnum) { + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS) { + if (++cnt == ilow - 1 ) { + convert_to_double_ex(data1); + val = Z_DVAL_PP(data1); + + zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1); + convert_to_double_ex(data1); + val += Z_DVAL_PP(data1); + val = val / 2.0; + break; + } + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); + } + } else { + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS) { + if (++cnt == ilow) { + convert_to_double_ex(data1); + val += Z_DVAL_PP(data1); + break; + } + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); + } + } + + RETURN_DOUBLE(val); +} +/* }}} */ + +/* {{{ proto float stats_stat_correlation(array arr1, array arr2) + Not documented */ +PHP_FUNCTION(stats_stat_correlation) +{ + zval **arg1, **arg2; + zval **data1, **data2; /* pointers to array entries */ + HashPosition pos1; /* hash iterator */ + HashPosition pos2; /* hash iterator */ + int xnum = 0; + int ynum = 0; + double sx = 0.0; + double sy = 0.0; + double sxx = 0.0; + double syy = 0.0; + double sxy = 0.0; + double mx; + double my; + double vx; + double vy; + double cc; + double rr; + + if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { + WRONG_PARAM_COUNT; + } + + convert_to_array_ex(arg1); + convert_to_array_ex(arg2); + + xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); + ynum = zend_hash_num_elements(Z_ARRVAL_PP(arg2)); + + if (xnum != ynum) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "Unequal number of X and Y coordinates"); + RETURN_FALSE; + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); + + while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS + && zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { + + convert_to_double_ex(data1); + convert_to_double_ex(data2); + + sx += Z_DVAL_PP(data1); + sxx += Z_DVAL_PP(data1) * Z_DVAL_PP(data1); + sy += Z_DVAL_PP(data2); + syy += Z_DVAL_PP(data2) * Z_DVAL_PP(data2); + sxy += Z_DVAL_PP(data1) * Z_DVAL_PP(data2); + + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); + zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); + } + + mx = sx / xnum; + my = sy / ynum; + vx = sxx - (xnum * mx * mx); + vy = syy - (ynum * my * my); + cc = sxy - (xnum * mx * my); + rr = cc / sqrt(vx * vy); + + RETURN_DOUBLE(rr); +} +/* }}} */ + +/* {{{ proto float stats_stat_binomial_coef(int x, int n) + Not documented */ +PHP_FUNCTION(stats_stat_binomial_coef) +{ + int i; + int n; + int x; + double bc = 1.0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ll", &x, &n) == FAILURE) { + RETURN_FALSE; + } + + for (i = 0; i < x; ++i) { + bc = (bc * (n - i)) / (i + 1); + } + + RETURN_DOUBLE(bc); +} +/* }}} */ + +/* {{{ proto float stats_stat_gennch(int n) + Not documented */ +PHP_FUNCTION(stats_stat_factorial) +{ + int n; + int i; + double f = 1; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "l", &n) == FAILURE) { + RETURN_FALSE; + } + + for (i = 1; i <= n; ++i) { + f *= i; + } + + RETURN_DOUBLE(f); +} +/* }}} */ + + +/* {{{ php_population_variance +*/ +static long double php_math_mean(zval *arr) +{ + double sum = 0.0; + zval **entry; + HashPosition pos; + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { + convert_to_double_ex(entry); + sum += Z_DVAL_PP(entry); + zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); + } + /* + we don't check whether the array has 0 elements. this is left to the caller - no need + to kill performance by checking on every level. + */ + return sum / zend_hash_num_elements(Z_ARRVAL_P(arr)); +} +/* }}} */ + + +/* {{{ php_population_variance +*/ +static long double php_population_variance(zval *arr, zend_bool sample) +{ + double mean, vr = 0.0; + zval **entry; + HashPosition pos; + int elements_num; + + elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr)); + + mean = php_math_mean(arr); + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { + double d; + convert_to_double_ex(entry); + d = Z_DVAL_PP(entry) - mean; + vr += d*d; + zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); + } + if (sample) { + --elements_num; + } + return (vr / elements_num); +} +/* }}} */ + + +/* {{{ proto float stats_variance(array a [, bool sample]) + Returns the population variance */ +PHP_FUNCTION(stats_variance) +{ + zval *arr; + zend_bool sample = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a|b", &arr, &sample) == FAILURE) { + return; + } + if (zend_hash_num_elements(Z_ARRVAL_P(arr)) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); + RETURN_FALSE; + } + if (sample && zend_hash_num_elements(Z_ARRVAL_P(arr)) == 1) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has only 1 element"); + RETURN_FALSE; + } + RETURN_DOUBLE(php_population_variance(arr, sample)); +} +/* }}} */ + +/* {{{ proto float stats_standard_deviation(array a[, bool sample = false]) + Returns the standard deviation */ +PHP_FUNCTION(stats_standard_deviation) +{ + zval *arr; + zend_bool sample = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a|b", &arr, &sample) == FAILURE) { + return; + } + if (zend_hash_num_elements(Z_ARRVAL_P(arr)) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); + RETURN_FALSE; + } + if (sample && zend_hash_num_elements(Z_ARRVAL_P(arr)) == 1) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has only 1 element"); + RETURN_FALSE; + } + RETURN_DOUBLE(sqrt(php_population_variance(arr, sample))); +} +/* }}} */ + + +/* {{{ proto float stats_absolute_deviation(array a) + Returns the absolute deviation of an array of values*/ +PHP_FUNCTION(stats_absolute_deviation) +{ + zval *arr; + double mean = 0.0, abs_dev = 0.0; + zval **entry; + HashPosition pos; + int elements_num; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { + return; + } + if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); + RETURN_FALSE; + } + + mean = php_math_mean(arr); + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { + convert_to_double_ex(entry); + abs_dev += fabs(Z_DVAL_PP(entry) - mean); + zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); + } + + RETURN_DOUBLE(abs_dev / elements_num); +} +/* }}} */ + +/* {{{ proto float stats_harmonic_mean(array a) + Returns the harmonic mean of an array of values */ +PHP_FUNCTION(stats_harmonic_mean) +{ + zval *arr; + double sum = 0.0; + zval **entry; + HashPosition pos; + int elements_num; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { + return; + } + if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); + RETURN_FALSE; + } + + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { + convert_to_double_ex(entry); + if (Z_DVAL_PP(entry) == 0) { + RETURN_LONG(0); + } + sum += 1 / Z_DVAL_PP(entry); + zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); + } + + RETURN_DOUBLE(elements_num / sum); +} +/* }}} */ + +/* {{{ proto float stats_skew(array a) + Computes the skewness of the data in the array */ +PHP_FUNCTION(stats_skew) +{ + zval *arr; + double mean, std_dev, skew = 0.0; + zval **entry; + HashPosition pos; + int elements_num, i = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { + return; + } + if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); + RETURN_FALSE; + } + + mean = php_math_mean(arr); + std_dev = sqrt(php_population_variance(arr, 0)); + + /* the calculation of the skewness is protected of value "explosion". a bit more + FP operations performed but more accurateness. + */ + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { + double tmp; + convert_to_double_ex(entry); + tmp = ((Z_DVAL_PP(entry) - mean) / std_dev); + skew += (tmp*tmp*tmp - skew) / (i + 1); + zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); + ++i; + } + + RETURN_DOUBLE(skew); +} +/* }}} */ + +/* {{{ proto float stats_kurtosis(array a) + Computes the kurtosis of the data in the array */ +PHP_FUNCTION(stats_kurtosis) +{ + zval *arr; + double mean, std_dev, avg = 0.0; + zval **entry; + HashPosition pos; + int elements_num, i = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { + return; + } + if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); + RETURN_FALSE; + } + + mean = php_math_mean(arr); + std_dev = sqrt(php_population_variance(arr, 0)); + + /* the calculation of the kurtosis is protected of value "explosion". a bit more + FP operations performed but more accurateness. + */ + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { + double tmp; + convert_to_double_ex(entry); + tmp = ((Z_DVAL_PP(entry) - mean) / std_dev); + avg += (tmp*tmp*tmp*tmp - avg) / (i + 1); + + zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); + ++i; + } + + RETURN_DOUBLE(avg - 3); +} +/* }}} */ + + +/* {{{ proto float stats_covariance(array a, array b) + Computes the covariance of two data sets */ +PHP_FUNCTION(stats_covariance) +{ + zval *arr_1, *arr_2; + double mean_1, mean_2, covar = 0.0; + zval **entry; + HashPosition pos_1, pos_2; + int elements_num, i = 0; + + if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "aa", &arr_1, &arr_2) == FAILURE) { + return; + } + if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr_1))) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The first array has zero elements"); + RETURN_FALSE; + } + if (zend_hash_num_elements(Z_ARRVAL_P(arr_2)) == 0) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The second array has zero elements"); + RETURN_FALSE; + } + if (elements_num != zend_hash_num_elements(Z_ARRVAL_P(arr_2))) { + php_error_docref(NULL TSRMLS_CC, E_WARNING, "The datasets are not of the same size"); + RETURN_FALSE; + } + + mean_1 = php_math_mean(arr_1); + mean_2 = php_math_mean(arr_2); + /* the calculation of the covariance is protected of value "explosion". a bit more + FP operations performed but more accurateness. + */ + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr_1), &pos_1); + zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr_2), &pos_2); + while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr_1), (void **)&entry, &pos_1) == SUCCESS) { + double tmp_1, tmp_2; + convert_to_double_ex(entry); + tmp_1 = Z_DVAL_PP(entry) - mean_1; + + if (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr_2), (void **)&entry, &pos_2) != SUCCESS) { + break; + } + convert_to_double_ex(entry); + tmp_2 = Z_DVAL_PP(entry) - mean_2; + + covar += (tmp_1 * tmp_2 - covar) / (i + 1); + + zend_hash_move_forward_ex(Z_ARRVAL_P(arr_1), &pos_1); + zend_hash_move_forward_ex(Z_ARRVAL_P(arr_2), &pos_2); + ++i; + } + + RETURN_DOUBLE(covar); +} +/* }}} */ + + +/* + * Local variables: + * tab-width: 4 + * c-basic-offset: 4 + * indent-tabs-mode: t + * End: + */ diff -dPNur stats-1.0.2/php_stats.h trunk/php_stats.h --- stats-1.0.2/php_stats.h 1970-01-01 01:00:00.000000000 +0100 +++ trunk/php_stats.h 2012-10-29 17:22:36.000000000 +0100 @@ -0,0 +1,129 @@ +/* + +----------------------------------------------------------------------+ + | PHP Version 5 | + +----------------------------------------------------------------------+ + | Copyright (c) 1997-2004 The PHP Group | + +----------------------------------------------------------------------+ + | This source file is subject to version 3.0 of the PHP license, | + | that is bundled with this package in the file LICENSE, and is | + | available through the world-wide-web at the following url: | + | http://www.php.net/license/3_0.txt. | + | If you did not receive a copy of the PHP license and are unable to | + | obtain it through the world-wide-web, please send a note to | + | license@php.net so we can mail you a copy immediately. | + +----------------------------------------------------------------------+ + | Author: Andrey Hristov | + +----------------------------------------------------------------------+ +*/ + +/* $Id: php_stats.h 256977 2008-04-08 16:13:17Z sfox $ */ + +#ifndef PHP_STATS_H +#define PHP_STATS_H + +extern zend_module_entry stats_module_entry; +#define phpext_stats_ptr &stats_module_entry + +#define PHP_STATS_VERSION "1.0.3-dev" + +#ifdef PHP_WIN32 +#define PHP_STATS_API __declspec(dllexport) +#else +#define PHP_STATS_API +#endif + + +PHP_MINFO_FUNCTION(stats); + +PHP_FUNCTION(stats_bin_counts); +PHP_FUNCTION(stats_cdf_t); +PHP_FUNCTION(stats_cdf_normal); +PHP_FUNCTION(stats_cdf_gamma); +PHP_FUNCTION(stats_cdf_chisquare); +PHP_FUNCTION(stats_cdf_beta); +PHP_FUNCTION(stats_cdf_binomial); +PHP_FUNCTION(stats_cdf_noncentral_chisquare); +PHP_FUNCTION(stats_cdf_f); +PHP_FUNCTION(stats_cdf_noncentral_f); +PHP_FUNCTION(stats_cdf_noncentral_t); +PHP_FUNCTION(stats_cdf_negative_binomial); +PHP_FUNCTION(stats_cdf_poisson); +PHP_FUNCTION(stats_cdf_laplace); +PHP_FUNCTION(stats_cdf_cauchy); +PHP_FUNCTION(stats_cdf_logistic); +PHP_FUNCTION(stats_cdf_weibull); +PHP_FUNCTION(stats_cdf_uniform); +PHP_FUNCTION(stats_cdf_exponential); +PHP_FUNCTION(stats_rand_setall); +PHP_FUNCTION(stats_rand_getsd); +PHP_FUNCTION(stats_rand_gen_iuniform); +PHP_FUNCTION(stats_rand_gen_funiform); +PHP_FUNCTION(stats_rand_ignlgi); +PHP_FUNCTION(stats_rand_ranf); +PHP_FUNCTION(stats_rand_gen_beta); +PHP_FUNCTION(stats_rand_gen_chisquare); +PHP_FUNCTION(stats_rand_gen_exponential); +PHP_FUNCTION(stats_rand_gen_f); +PHP_FUNCTION(stats_rand_gen_gamma); +PHP_FUNCTION(stats_rand_gen_noncentral_chisquare); +PHP_FUNCTION(stats_rand_gen_noncenral_f); +PHP_FUNCTION(stats_rand_gen_normal); +PHP_FUNCTION(stats_rand_phrase_to_seeds); +PHP_FUNCTION(stats_rand_ibinomial); +PHP_FUNCTION(stats_rand_ibinomial_negative); +PHP_FUNCTION(stats_rand_gen_ipoisson); +PHP_FUNCTION(stats_rand_gen_noncentral_t); +PHP_FUNCTION(stats_rand_gen_t); +PHP_FUNCTION(stats_dens_normal); +PHP_FUNCTION(stats_dens_cauchy); +PHP_FUNCTION(stats_dens_laplace); +PHP_FUNCTION(stats_dens_logistic); +PHP_FUNCTION(stats_dens_beta); +PHP_FUNCTION(stats_dens_weibull); +PHP_FUNCTION(stats_dens_uniform); +PHP_FUNCTION(stats_dens_chisquare); +PHP_FUNCTION(stats_dens_t); +PHP_FUNCTION(stats_dens_gamma); +PHP_FUNCTION(stats_dens_exponential); +PHP_FUNCTION(stats_dens_f); +PHP_FUNCTION(stats_dens_pmf_binomial); +PHP_FUNCTION(stats_dens_pmf_poisson); +PHP_FUNCTION(stats_dens_pmf_negative_binomial); +PHP_FUNCTION(stats_dens_pmf_hypergeometric); +PHP_FUNCTION(stats_stat_powersum); +PHP_FUNCTION(stats_stat_innerproduct); +PHP_FUNCTION(stats_stat_independent_t); +PHP_FUNCTION(stats_stat_paired_t); +PHP_FUNCTION(stats_stat_percentile); +PHP_FUNCTION(stats_stat_correlation); +PHP_FUNCTION(stats_stat_binomial_coef); +PHP_FUNCTION(stats_stat_factorial); +PHP_FUNCTION(stats_absolute_deviation); +PHP_FUNCTION(stats_standard_deviation); +PHP_FUNCTION(stats_variance); +PHP_FUNCTION(stats_harmonic_mean); +PHP_FUNCTION(stats_skew); +PHP_FUNCTION(stats_kurtosis); +PHP_FUNCTION(stats_covariance); + + +#ifdef ZTS +#define STATS_D zend_stats_globals *stats_globals +#define STATS_G(v) (stats_globals->v) +#define STATS_FETCH() zend_stats_globals *stats_globals = ts_resource(stats_globals_id) +#else +#define STATS_D +#define STATS_G(v) (stats_globals.v) +#define STATS_FETCH() +#endif + +#endif /* PHP_STATS_H */ + + +/* + * Local variables: + * tab-width: 4 + * c-basic-offset: 4 + * indent-tabs-mode: t + * End: + */ diff -dPNur stats-1.0.2/statistics.c trunk/statistics.c --- stats-1.0.2/statistics.c 2006-05-31 19:24:26.000000000 +0200 +++ trunk/statistics.c 1970-01-01 01:00:00.000000000 +0100 @@ -1,3728 +0,0 @@ -/* - +----------------------------------------------------------------------+ - | PHP Version 5 | - +----------------------------------------------------------------------+ - | Copyright (c) 1997-2004 The PHP Group | - +----------------------------------------------------------------------+ - | This source file is subject to version 3.0 of the PHP license, | - | that is bundled with this package in the file LICENSE, and is | - | available through the world-wide-web at the following url: | - | http://www.php.net/license/3_0.txt. | - | If you did not receive a copy of the PHP license and are unable to | - | obtain it through the world-wide-web, please send a note to | - | license@php.net so we can mail you a copy immediately. | - +----------------------------------------------------------------------+ - | Author: Andrey Hristov | - +----------------------------------------------------------------------+ -*/ - -/* $Id: statistics.c,v 1.11 2006/05/30 18:02:36 andrey Exp $ */ - - -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -#include "php.h" -#include "php_statistics.h" -#include "ext/standard/info.h" -#include "ext/standard/head.h" -#include -#include -#include -#include -#include "randlib.h" -#include "cdflib.h" - -#define STATS_PI 3.14159265358979323846 - - -#ifdef PHP_WIN32 -extern double fd_lgamma(double x); -#define lgamma fd_lgamma -#endif - -static double logistic_quantile(double p); -static double logistic_cdf(double x); -static double cauchy_quantile(double p); -static double cauchy_cdf(double x); -static double laplace_quantile(double p); -static double laplace_cdf(double x); -static double exponential_quantile(double p); -static double exponential_cdf(double x); -static double binom(double x, double n); - -zend_function_entry statistics_functions[] = { - PHP_FE(stats_cdf_t, NULL) - PHP_FE(stats_cdf_normal, NULL) - PHP_FE(stats_cdf_gamma, NULL) - PHP_FE(stats_cdf_chisquare, NULL) - PHP_FE(stats_cdf_beta, NULL) - PHP_FE(stats_cdf_binomial, NULL) - PHP_FE(stats_cdf_noncentral_chisquare,NULL) - PHP_FE(stats_cdf_f, NULL) - PHP_FE(stats_cdf_noncentral_f, NULL) - PHP_FE(stats_cdf_noncentral_t, NULL) - PHP_FE(stats_cdf_negative_binomial, NULL) - PHP_FE(stats_cdf_poisson, NULL) - PHP_FE(stats_cdf_laplace, NULL) - PHP_FE(stats_cdf_cauchy, NULL) - PHP_FE(stats_cdf_logistic, NULL) - PHP_FE(stats_cdf_weibull, NULL) - PHP_FE(stats_cdf_uniform, NULL) - PHP_FE(stats_cdf_exponential, NULL) - PHP_FE(stats_rand_setall, NULL) - PHP_FE(stats_rand_getsd, NULL) - PHP_FE(stats_rand_gen_iuniform, NULL) - PHP_FE(stats_rand_gen_funiform, NULL) - PHP_FE(stats_rand_ignlgi, NULL) - PHP_FE(stats_rand_ranf, NULL) - PHP_FE(stats_rand_gen_beta, NULL) - PHP_FE(stats_rand_gen_chisquare, NULL) - PHP_FE(stats_rand_gen_exponential, NULL) - PHP_FE(stats_rand_gen_f, NULL) - PHP_FE(stats_rand_gen_gamma, NULL) - PHP_FE(stats_rand_gen_noncentral_chisquare,NULL) - PHP_FE(stats_rand_gen_noncenral_f, NULL) - PHP_FE(stats_rand_gen_normal, NULL) - PHP_FE(stats_rand_phrase_to_seeds, NULL) - PHP_FE(stats_rand_ibinomial, NULL) - PHP_FE(stats_rand_ibinomial_negative,NULL) - PHP_FE(stats_rand_gen_ipoisson, NULL) - PHP_FE(stats_rand_gen_noncentral_t, NULL) - PHP_FE(stats_rand_gen_t, NULL) - PHP_FE(stats_dens_normal, NULL) - PHP_FE(stats_dens_cauchy, NULL) - PHP_FE(stats_dens_laplace, NULL) - PHP_FE(stats_dens_logistic, NULL) - PHP_FE(stats_dens_beta, NULL) - PHP_FE(stats_dens_weibull, NULL) - PHP_FE(stats_dens_uniform, NULL) - PHP_FE(stats_dens_chisquare, NULL) - PHP_FE(stats_dens_t, NULL) - PHP_FE(stats_dens_gamma, NULL) - PHP_FE(stats_dens_exponential, NULL) - PHP_FE(stats_dens_f, NULL) - PHP_FE(stats_dens_pmf_binomial, NULL) - PHP_FE(stats_dens_pmf_poisson, NULL) - PHP_FE(stats_dens_pmf_negative_binomial,NULL) - PHP_FE(stats_dens_pmf_hypergeometric, NULL) - PHP_FE(stats_stat_powersum, NULL) - PHP_FE(stats_stat_innerproduct, NULL) - PHP_FE(stats_stat_independent_t, NULL) - PHP_FE(stats_stat_paired_t, NULL) - PHP_FE(stats_stat_percentile, NULL) - PHP_FE(stats_stat_correlation, NULL) - PHP_FE(stats_stat_binomial_coef, NULL) - PHP_FE(stats_stat_factorial, NULL) - PHP_FE(stats_standard_deviation, NULL) - PHP_FE(stats_absolute_deviation, NULL) - PHP_FE(stats_variance, NULL) - PHP_FE(stats_harmonic_mean, NULL) - PHP_FE(stats_skew, NULL) - PHP_FE(stats_kurtosis, NULL) - PHP_FE(stats_covariance, NULL) - {NULL, NULL, NULL} -}; - -zend_module_entry stats_module_entry = { - STANDARD_MODULE_HEADER, - "stats", - statistics_functions, - NULL, - NULL, - NULL, - NULL, - PHP_MINFO(stats), - "1.1", - STANDARD_MODULE_PROPERTIES, -}; - -#ifdef COMPILE_DL_STATS -ZEND_GET_MODULE(stats) -#endif - - -PHP_MINFO_FUNCTION(stats) -{ - php_info_print_table_start(); - php_info_print_table_row(2, "Statistics Support", "enabled"); - php_info_print_table_end(); -} - - - -/* Numbers are always smaller than strings int this function as it - * anyway doesn't make much sense to compare two different data types. - * This keeps it consistant and simple. - * - * This is not correct any more, depends on what compare_func is set to. - */ -static int stats_array_data_compare(const void *a, const void *b TSRMLS_DC) -{ - Bucket *f; - Bucket *s; - pval result; - pval *first; - pval *second; - - f = *((Bucket **) a); - s = *((Bucket **) b); - - first = *((pval **) f->pData); - second = *((pval **) s->pData); - - if (numeric_compare_function(&result, first, second TSRMLS_CC) == FAILURE) { - return 0; - } - - if (Z_TYPE(result) == IS_DOUBLE) { - if (Z_DVAL(result) < 0) { - return -1; - } else if (Z_DVAL(result) > 0) { - return 1; - } else { - return 0; - } - } - - convert_to_long(&result); - - if (Z_LVAL(result) < 0) { - return -1; - } else if (Z_LVAL(result) > 0) { - return 1; - } - - return 0; -} - - - -/**************************************/ -/* Cumulative Distributions Functions */ -/**************************************/ - -/****************************************************** - Cumulative Distribution Function - T distribution - - Function - - - Calculates any one parameter of the t distribution given - values for the others. - - Arguments - - WHICH --> Integer indicating which argument - values is to be calculated from the others. - Legal range: 1..3 - iwhich = 1 : Calculate P and Q from T and DF - iwhich = 2 : Calculate T from P,Q and DF - iwhich = 3 : Calculate DF from P,Q and T - - P <--> The integral from -infinity to t of the t-density. - Input range: (0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - T <--> Upper limit of integration of the t-density. - Input range: ( -infinity, +infinity). - Search range: [ -1E100, 1E100 ] - - DF <--> Degrees of freedom of the t-distribution. - Input range: (0 , +infinity). - Search range: [1e-100, 1E10] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.5.27 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to reduce the computation - of the cumulative distribution function to that of an incomplete - beta. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. -******************************************************/ - -/* {{{ proto float stats_cdf_t(float par1, float par2, int which) - Calculates any one parameter of the T distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_t) -{ - double arg1; - double arg2; - double df; - double bound; - double p; - double q; - double t; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddl", &arg1, &arg2, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 3) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); - RETURN_FALSE; - } - - if (which < 3 ) { - df = arg2; - } else { - t = arg2; - } - if (which == 1) { - t = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdft((int *)&which, &p, &q, &t, &df, &status, &bound); - - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(t); - case 3: RETURN_DOUBLE(df); - } - RETURN_FALSE; /* should never be reached */ -} -/* }}} */ - -/********************************************************************* - Cumulative Distribution Function NORmal distribution - - Calculates any one parameter of the normal - distribution given values for the others. - - - Arguments - - - WHICH --> Integer indicating which of the next parameter - values is to be calculated using values of the others. - Legal range: 1..4 - iwhich = 1 : Calculate P and Q from X,MEAN and SD - iwhich = 2 : Calculate X from P,Q,MEAN and SD - iwhich = 3 : Calculate MEAN from P,Q,X and SD - iwhich = 4 : Calculate SD from P,Q,X and MEAN - - P <--> The integral from -infinity to X of the normal density. - Input range: (0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - X < --> Upper limit of integration of the normal-density. - Input range: ( -infinity, +infinity) - - MEAN <--> The mean of the normal density. - Input range: (-infinity, +infinity) - - SD <--> Standard Deviation of the normal density. - Input range: (0, +infinity). - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - A slightly modified version of ANORM from - - Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN - Package of Special Function Routines and Test Drivers" - acm Transactions on Mathematical Software. 19, 22-32. - - is used to calulate the cumulative standard normal distribution. - - The rational functions from pages 90-95 of Kennedy and Gentle, - Statistical Computing, Marcel Dekker, NY, 1980 are used as - starting values to Newton's Iterations which compute the inverse - standard normal. Therefore no searches are necessary for any - parameter. - - For X < -15, the asymptotic expansion for the normal is used as - the starting value in finding the inverse standard normal. - This is formula 26.2.12 of Abramowitz and Stegun. - - Note - - The normal density is proportional to - exp( - 0.5 * (( X - MEAN)/SD)**2) -***********************************************************************/ - -/* {{{ proto float stats_stat_gennch(float par1, float par2, float par3, int which) - Calculates any one parameter of the normal distribution given values for thee others. */ -PHP_FUNCTION(stats_cdf_normal) -{ - double arg1; - double arg2; - double arg3; - double sd; - double bound; - double p; - double q; - double x; - double mean; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - sd = arg3; - } else { - mean = arg3; - } - - if (which < 3) { - mean = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfnor((int *)&which, &p, &q, &x, &mean, &sd, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(mean); - case 4: RETURN_DOUBLE(sd); - } - RETURN_FALSE; /* should never be reached */ -} -/* }}} */ - - -/********************************************************************* - Cumulative Distribution Function - GAMma Distribution - - - Function - - - Calculates any one parameter of the gamma - distribution given values for the others. - - - Arguments - - - WHICH --> Integer indicating which of the next four argument - values is to be calculated from the others. - Legal range: 1..4 - iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE - iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE - iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE - iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE - - P <--> The integral from 0 to X of the gamma density. - Input range: [0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - X <--> The upper limit of integration of the gamma density. - Input range: [0, +infinity). - Search range: [0,1E100] - - SHAPE <--> The shape parameter of the gamma density. - Input range: (0, +infinity). - Search range: [1E-100,1E100] - - SCALE <--> The scale parameter of the gamma density. - Input range: (0, +infinity). - Search range: (1E-100,1E100] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - 10 if the gamma or inverse gamma routine cannot - compute the answer. Usually happens only for - X and SHAPE very large (gt 1E10 or more) - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - - Method - - - Cumulative distribution function (P) is calculated directly by - the code associated with: - - DiDinato, A. R. and Morris, A. H. Computation of the incomplete - gamma function ratios and their inverse. ACM Trans. Math. - Softw. 12 (1986), 377-393. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - - - Note - - - - The gamma density is proportional to - T**(SHAPE - 1) * EXP(- SCALE * T) -**************************************************************************/ -/* {{{ proto float stats_cdf_gamma(float par1, float par2, float par3, int which) - Calculates any one parameter of the gamma distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_gamma) -{ - double arg1; - double arg2; - double arg3; - double bound; - double p; - double q; - double x; - double shape; - double scale; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - scale = arg3; - } else { - shape = arg3; - } - - if (which < 3) { - shape = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfgam((int *)&which, &p, &q, &x, &shape, &scale, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(shape); - case 4: RETURN_DOUBLE(scale); - } - RETURN_FALSE; /* should never be reached */ -} -/* }}} */ - -/***************************************************************** - Cumulative Distribution Function - CHI-Square distribution - - Function - - Calculates any one parameter of the chi-square - distribution given values for the others. - - Arguments - - - WHICH --> Integer indicating which of the next three argument - values is to be calculated from the others. - Legal range: 1..3 - iwhich = 1 : Calculate P and Q from X and DF - iwhich = 2 : Calculate X from P,Q and DF - iwhich = 3 : Calculate DF from P,Q and X - - P <--> The integral from 0 to X of the chi-square - distribution. - Input range: [0, 1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - X <--> Upper limit of integration of the non-central - chi-square distribution. - Input range: [0, +infinity). - Search range: [0,1E100] - - DF <--> Degrees of freedom of the - chi-square distribution. - Input range: (0, +infinity). - Search range: [ 1E-100, 1E100] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - 10 indicates error returned from cumgam. See - references in cdfgam - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - - Method - - - Formula 26.4.19 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to reduce the chisqure - distribution to the incomplete distribution. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. -*****************************************************************/ -/* {{{ proto float stats_cdf_chisquare(float par1, float par2, float par3, int which) - Calculates any one parameter of the chi-square distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_chisquare) -{ - double arg1; - double arg2; - double bound; - double p; - double q; - double x; - double df; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddl", &arg1, &arg2, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 3) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); - RETURN_FALSE; - } - - if (which < 3 ) { - df = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfchi((int *)&which, &p, &q, &x, &df, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(df); - } - RETURN_FALSE; /* should never be here */ -} -/* }}} */ - -/******************************************************************* - Cumulative Distribution Function - BETa Distribution - - Function - - Calculates any one parameter of the beta distribution given - values for the others. - - Arguments - - WHICH --> Integer indicating which of the next four argument - values is to be calculated from the others. - Legal range: 1..4 - iwhich = 1 : Calculate P and Q from X,Y,A and B - iwhich = 2 : Calculate X and Y from P,Q,A and B - iwhich = 3 : Calculate A from P,Q,X,Y and B - iwhich = 4 : Calculate B from P,Q,X,Y and A - - P <--> The integral from 0 to X of the chi-square - distribution. - Input range: [0, 1]. - - Q <--> 1-P. - Input range: [0, 1]. - P + Q = 1.0. - - X <--> Upper limit of integration of beta density. - Input range: [0,1]. - Search range: [0,1] - - Y <--> 1-X. - Input range: [0,1]. - Search range: [0,1] - X + Y = 1.0. - - A <--> The first parameter of the beta density. - Input range: (0, +infinity). - Search range: [1D-100,1D100] - - B <--> The second parameter of the beta density. - Input range: (0, +infinity). - Search range: [1D-100,1D100] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - 4 if X + Y .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Cumulative distribution function (P) is calculated directly by - code associated with the following reference. - - DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant - Digit Computation of the Incomplete Beta Function Ratios. ACM - Trans. Math. Softw. 18 (1993), 360-373. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - - Note - - The beta density is proportional to - t^(A-1) * (1-t)^(B-1) - -*******************************************************************/ - -/* {{{ proto float stats_cdf_beta(float par1, float par2, float par3, int which) - Calculates any one parameter of the beta distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_beta) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double x; - double bound; - double y; - double a; - double b; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - - if (which < 4) { - b = arg3; - } else { - a = arg3; - } - - if (which < 3) { - a = arg2; - } else { - x = arg2; - y = 1.0 - x; - } - - if (which == 1) { - x = arg1; - y = 1.0 - x; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfbet((int *)&which, &p, &q, &x, &y, &a, &b, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(a); - case 4: RETURN_DOUBLE(b); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/********************************************************************* - Cumulative Distribution Function - BINomial distribution - - Function - - Calculates any one parameter of the binomial - distribution given values for the others. - - Arguments - - WHICH --> Integer indicating which of the next four argument - values is to be calculated from the others. - Legal range: 1..4 - iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR - iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR - iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR - iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN - - P <--> The cumulation from 0 to S of the binomial distribution. - (Probablility of S or fewer successes in XN trials each - with probability of success PR.) - Input range: [0,1]. - - Q <--> 1-P. - Input range: [0, 1]. - P + Q = 1.0. - - S <--> The number of successes observed. - Input range: [0, XN] - Search range: [0, XN] - - XN <--> The number of binomial trials. - Input range: (0, +infinity). - Search range: [1E-100, 1E100] - - PR <--> The probability of success in each binomial trial. - Input range: [0,1]. - Search range: [0,1] - - OMPR <--> 1-PR - Input range: [0,1]. - Search range: [0,1] - PR + OMPR = 1.0 - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - 4 if PR + OMPR .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.5.24 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to reduce the binomial - distribution to the cumulative incomplete beta distribution. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - -*********************************************************************/ - -/* {{{ proto float stats_cdf_binomial(float par1, float par2, float par3, int which) - Calculates any one parameter of the binomial distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_binomial) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double xn; - double bound; - double sn; - double pr; - double ompr; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - - if (which < 4) { - pr = arg3; - ompr = 1.0 - pr; - } else { - xn = arg3; - } - - if (which < 3) { - xn = arg2; - } else { - sn = arg2; - } - - if (which == 1) { - sn = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfbin((int *)&which, &p, &q, &sn, &xn, &pr, &ompr, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in binomialcdf"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(sn); - case 3: RETURN_DOUBLE(xn); - case 4: RETURN_DOUBLE(pr); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/***************************************************************** - Cumulative Distribution Function - Non-central Chi-Square - - Function - - Calculates any one parameter of the non-central chi-square - distribution given values for the others. - - Arguments - - WHICH --> Integer indicating which of the next three argument - values is to be calculated from the others. - Input range: 1..4 - iwhich = 1 : Calculate P and Q from X and DF - iwhich = 2 : Calculate X from P,DF and PNONC - iwhich = 3 : Calculate DF from P,X and PNONC - iwhich = 3 : Calculate PNONC from P,X and DF - - P <--> The integral from 0 to X of the non-central chi-square - distribution. - Input range: [0, 1-1E-16). - - Q <--> 1-P. - Q is not used by this subroutine and is only included - for similarity with other cdf* routines. - - X <--> Upper limit of integration of the non-central - chi-square distribution. - Input range: [0, +infinity). - Search range: [0,1E100] - - DF <--> Degrees of freedom of the non-central - chi-square distribution. - Input range: (0, +infinity). - Search range: [ 1E-100, 1E100] - - PNONC <--> Non-centrality parameter of the non-central - chi-square distribution. - Input range: [0, +infinity). - Search range: [0,1E4] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.4.25 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to compute the cumulative - distribution function. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - - WARNING - - The computation time required for this routine is proportional - to the noncentrality parameter (PNONC). Very large values of - this parameter can consume immense computer resources. This is - why the search range is bounded by 10,000. - -*****************************************************************/ - -/* {{{ proto float stats_cdf_noncentral_chisquare(float par1, float par2, float par3, int which) - Calculates any one parameter of the non-central chi-square distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_noncentral_chisquare) -{ - double arg1; - double arg2; - double arg3; - double p; - double bound; - double q; - double x; - double df; - double pnonc; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - pnonc = arg3; - } else { - df = arg3; - } - - if (which < 3) { - df = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfchn((int *)&which, &p, &q, &x, &df, &pnonc, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdfchn"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(df); - case 4: RETURN_DOUBLE(pnonc); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/************************************************************** - Cumulative Distribution Function F distribution - - Function - - Calculates any one parameter of the F distribution - given values for the others. - - Arguments - - WHICH --> Integer indicating which of the next four argument - values is to be calculated from the others. - Legal range: 1..4 - iwhich = 1 : Calculate P and Q from F,DFN and DFD - iwhich = 2 : Calculate F from P,Q,DFN and DFD - iwhich = 3 : Calculate DFN from P,Q,F and DFD - iwhich = 4 : Calculate DFD from P,Q,F and DFN - - P <--> The integral from 0 to F of the f-density. - Input range: [0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - F <--> Upper limit of integration of the f-density. - Input range: [0, +infinity). - Search range: [0,1E100] - - DFN < --> Degrees of freedom of the numerator sum of squares. - Input range: (0, +infinity). - Search range: [ 1E-100, 1E100] - - DFD < --> Degrees of freedom of the denominator sum of squares. - Input range: (0, +infinity). - Search range: [ 1E-100, 1E100] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.6.2 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to reduce the computation - of the cumulative distribution function for the F variate to - that of an incomplete beta. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - - WARNING - - The value of the cumulative F distribution is not necessarily - monotone in either degrees of freedom. There thus may be two - values that provide a given CDF value. This routine assumes - monotonicity and will find an arbitrary one of the two values. - -**********************************************************************/ -/* {{{ proto float stats_cdf_f(float par1, float par2, float par3, int which) - Calculates any one parameter of the F distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_f) -{ - double arg1; - double arg2; - double arg3; - double p; - double bound; - double dfn; - double q; - double f; - double dfd; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - if (which < 4) { - dfd = arg3; - } else { - dfn = arg3; - } - if (which < 3) { - dfn = arg2; - } else { - f = arg2; - } - if (which == 1) { - f = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdff((int *)&which, &p, &q, &f, &dfn, &dfd, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdff"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(f); - case 3: RETURN_DOUBLE(dfn); - case 4: RETURN_DOUBLE(dfd); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/*********************************************************************** - Cumulative Distribution Function - Non-central F distribution - - Function - - Calculates any one parameter of the Non-central F - distribution given values for the others. - - Arguments - - WHICH --> Integer indicating which of the next five argument - values is to be calculated from the others. - Legal range: 1..5 - iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC - iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC - iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC - iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC - iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD - - P <--> The integral from 0 to F of the non-central f-density. - Input range: [0,1-1E-16). - - Q <--> 1-P. - Q is not used by this subroutine and is only included - for similarity with other cdf* routines. - - F <--> Upper limit of integration of the non-central f-density. - Input range: [0, +infinity). - Search range: [0,1E100] - - DFN < --> Degrees of freedom of the numerator sum of squares. - Input range: (0, +infinity). - Search range: [ 1E-100, 1E100] - - DFD < --> Degrees of freedom of the denominator sum of squares. - Must be in range: (0, +infinity). - Input range: (0, +infinity). - Search range: [ 1E-100, 1E100] - - PNONC <-> The non-centrality parameter - Input range: [0,infinity) - Search range: [0,1E4] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.6.20 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to compute the cumulative - distribution function. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - - WARNING - - The computation time required for this routine is proportional - to the noncentrality parameter (PNONC). Very large values of - this parameter can consume immense computer resources. This is - why the search range is bounded by 10,000. - - WARNING - - The value of the cumulative noncentral F distribution is not - necessarily monotone in either degrees of freedom. There thus - may be two values that provide a given CDF value. This routine - assumes monotonicity and will find an arbitrary one of the two - values. - -***********************************************************************/ - -/* {{{ proto float stats_cdf_noncentral_f(float par1, float par2, float par3, float par4, int which) - Calculates any one parameter of the Non-central F distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_noncentral_f) -{ - double arg1; - double arg2; - double arg3; - double arg4; - double p; - double q; - double f; - double dfn; - double dfd; - double pnonc; - double bound; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddddl", &arg1, &arg2, &arg3, &arg4, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 5) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fifth parameter should be in the 1..5 range"); - RETURN_FALSE; - } - - if (which < 5) { - pnonc = arg4; - } else { - dfd = arg4; - } - - if (which < 4) { - dfd = arg3; - } else { - dfn = arg3; - } - - if (which < 3) { - dfn = arg2; - } else { - f = arg2; - } - - if (which == 1) { - f = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdffnc((int *)&which, &p, &q, &f, &dfn, &dfd, &pnonc, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdffnc"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(f); - case 3: RETURN_DOUBLE(dfn); - case 4: RETURN_DOUBLE(dfd); - case 5: RETURN_DOUBLE(pnonc); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/************************************************************************ - Cumulative Distribution Function Non-Central T distribution - - Function - - Calculates any one parameter of the noncentral t distribution give - values for the others. - - Arguments - - WHICH --> Integer indicating which argument - values is to be calculated from the others. - Legal range: 1..3 - iwhich = 1 : Calculate P and Q from T,DF,PNONC - iwhich = 2 : Calculate T from P,Q,DF,PNONC - iwhich = 3 : Calculate DF from P,Q,T - iwhich = 4 : Calculate PNONC from P,Q,DF,T - - P <--> The integral from -infinity to t of the noncentral t-den - Input range: (0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - T <--> Upper limit of integration of the noncentral t-density. - Input range: ( -infinity, +infinity). - Search range: [ -1E100, 1E100 ] - - DF <--> Degrees of freedom of the noncentral t-distribution. - Input range: (0 , +infinity). - Search range: [1e-100, 1E10] - - PNONC <--> Noncentrality parameter of the noncentral t-distributio - Input range: [-infinity , +infinity). - Search range: [-1e4, 1E4] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Upper tail of the cumulative noncentral t is calculated usin - formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuou - Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995) - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - -************************************************************************/ - -/* {{{ proto float stats_stat_noncentral_t(float par1, float par2, float par3, int which) - Calculates any one parameter of the noncentral t distribution give values for the others. */ -PHP_FUNCTION(stats_cdf_noncentral_t) -{ - double arg1; - double arg2; - double arg3; - double pnonc; - double bound; - double p; - double q; - double t; - double df; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - pnonc = arg3; - } else { - df = arg3; - } - if (which < 3) { - df = arg2; - } else { - t = arg2; - } - - if (which == 1) { - t = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdftnc((int *)&which, &p, &q, &t, &df, &pnonc, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(t); - case 3: RETURN_DOUBLE(df); - case 4: RETURN_DOUBLE(pnonc); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/************************************************************************ - Cumulative Distribution Function Negative BiNomial distribution - - Function - - Calculates any one parameter of the negative binomial - distribution given values for the others. - - The cumulative negative binomial distribution returns the - probability that there will be F or fewer failures before the - XNth success in binomial trials each of which has probability of - success PR. - - The individual term of the negative binomial is the probability of - S failures before XN successes and is - Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S - - Arguments - - WHICH --> Integer indicating which of the next four argument - values is to be calculated from the others. - Legal range: 1..4 - iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR - iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR - iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR - iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN - - P <--> The cumulation from 0 to S of the negative - binomial distribution. - Input range: [0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - S <--> The upper limit of cumulation of the binomial distribution. - There are F or fewer failures before the XNth success. - Input range: [0, +infinity). - Search range: [0, 1E100] - - XN <--> The number of successes. - Input range: [0, +infinity). - Search range: [0, 1E100] - - PR <--> The probability of success in each binomial trial. - Input range: [0,1]. - Search range: [0,1]. - - OMPR <--> 1-PR - Input range: [0,1]. - Search range: [0,1] - PR + OMPR = 1.0 - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - 4 if PR + OMPR .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.5.26 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to reduce calculation of - the cumulative distribution function to that of an incomplete - beta. - - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - -************************************************************************/ - -/* {{{ proto float stats_cdf_negative_binomial(float par1, float par2, float par3, int which) - Calculates any one parameter of the negative binomial distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_negative_binomial) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double bound; - double sn; - double xn; - double pr; - double ompr; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - pr = arg3; - ompr = 1.0 - pr; - } else { - xn = arg3; - } - - if (which < 3) { - xn = arg2; - } else { - sn = arg2; - } - - if (which == 1) { - sn = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfnbn((int *)&which, &p, &q, &sn, &xn, &pr, &ompr, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error in cdfnbn"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(sn); - case 3: RETURN_DOUBLE(xn); - case 4: RETURN_DOUBLE(pr); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/********************************************************************** - Cumulative Distribution Function POIsson distribution - - Function - - Calculates any one parameter of the Poisson - distribution given values for the others. - - Arguments - - WHICH --> Integer indicating which argument - value is to be calculated from the others. - Legal range: 1..3 - iwhich = 1 : Calculate P and Q from S and XLAM - iwhich = 2 : Calculate A from P,Q and XLAM - iwhich = 3 : Calculate XLAM from P,Q and S - - P <--> The cumulation from 0 to S of the poisson density. - Input range: [0,1]. - - Q <--> 1-P. - Input range: (0, 1]. - P + Q = 1.0. - - S <--> Upper limit of cumulation of the Poisson. - Input range: [0, +infinity). - Search range: [0,1E100] - - XLAM <--> Mean of the Poisson distribution. - Input range: [0, +infinity). - Search range: [0,1E100] - - STATUS <-- 0 if calculation completed correctly - -I if input parameter number I is out of range - 1 if answer appears to be lower than lowest - search bound - 2 if answer appears to be higher than greatest - search bound - 3 if P + Q .ne. 1 - - BOUND <-- Undefined if STATUS is 0 - - Bound exceeded by parameter number I if STATUS - is negative. - - Lower search bound if STATUS is 1. - - Upper search bound if STATUS is 2. - - Method - - Formula 26.4.21 of Abramowitz and Stegun, Handbook of - Mathematical Functions (1966) is used to reduce the computation - of the cumulative distribution function to that of computing a - chi-square, hence an incomplete gamma function. - - Cumulative distribution function (P) is calculated directly. - Computation of other parameters involve a seach for a value that - produces the desired value of P. The search relies on the - monotinicity of P with the other parameter. - -**********************************************************************/ - -/* {{{ proto float stats_cdf_poisson(float par1, float par2, float par3, int which) - Calculates any one parameter of the Poisson distribution given values for the others. */ -PHP_FUNCTION(stats_cdf_poisson) -{ - double arg1; - double arg2; - double p; - double q; - double x; - double xlam; - double bound; - long which; - int status = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddl", &arg1, &arg2, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 3) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); - RETURN_FALSE; - } - - if (which < 3 ) { - xlam = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - cdfpoi((int *)&which, &p, &q, &x, &xlam, &status, &bound); - if (status != 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Computation Error"); - RETURN_FALSE; - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(xlam); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - - -static double laplace_quantile(double p) -{ - if (p <= 0.5) { - return log(2.0*p); - } else { - return (-log(2.0*(1.0-p))); - } -} - -static double laplace_cdf(double x) -{ - if (x <= 0) { - return (0.5*exp(x)); - } else { - return (1.0 - 0.5*exp(-x)); - } -} - - -/* {{{ proto float stats_cdf_laplace(float par1, float par2, float par3, int which) - Not documented */ -PHP_FUNCTION(stats_cdf_laplace) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double x; - double t; - double mean; - double sd; - long which; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - sd = arg3; - } else { - mean = arg3; - } - - if (which < 3) { - mean = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - if (which == 1) { - t = (x - mean) / sd; - p = laplace_cdf(t); - } else { - t = laplace_quantile(p); - } - - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(mean + (sd * t)); - case 3: RETURN_DOUBLE(x - (sd * t)); - case 4: RETURN_DOUBLE((x - mean) / t); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -static double cauchy_quantile(double p) -{ - return (tan(STATS_PI*(p-0.5))); -} - -static double cauchy_cdf (double x) -{ - return (0.5+(atan(x)/STATS_PI)); -} - -/* {{{ proto float stats_cdf_cauchy(float par1, float par2, float par3, int which) - Not documented */ -PHP_FUNCTION(stats_cdf_cauchy) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double x; - double t; - double mean; - double sd; - long which; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - sd = arg3; - } else { - mean = arg3; - } - - if (which < 3) { - mean = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - if (which == 1) { - t = (x - mean) / sd; - p = cauchy_cdf(t); - } else { - t = cauchy_quantile(p); - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(mean + (sd * t)); - case 3: RETURN_DOUBLE(x - (sd * t)); - case 4: RETURN_DOUBLE((x - mean) / t); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -static double logistic_cdf(double x) -{ - return (1.0/(1.0+exp(-x))); -} - -static double logistic_quantile (double p) -{ - return (log(p/(1.0-p))); -} - -/* {{{ proto float stats_cdf_logistic(float par1, float par2, float par3, int which) - Not documented */ -PHP_FUNCTION(stats_cdf_logistic) -{ - double arg1; - double arg2; - double arg3; - double sd; - double p; - double q; - double x; - double t; - double mean; - long which; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - sd = arg3; - } else { - mean = arg3; - } - - if (which < 3) { - mean = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - if (which == 1) { - t = (x - mean) / sd; - p = logistic_cdf(t); - } else { - t = logistic_quantile(p); - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(mean + (sd * t)); - case 3: RETURN_DOUBLE(x - (sd * t)); - case 4: RETURN_DOUBLE((x - mean) / t); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/* {{{ proto float stats_cdf_weibull(float par1, float par2, float par3, int which) - Not documented */ -PHP_FUNCTION(stats_cdf_weibull) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double x; - double a; - double b; - long which; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - b = arg3; - } else { - a = arg3; - } - - if (which < 3) { - a = arg2; - } else { - x = arg2; - } - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - if (which == 1) { - p = 1 - exp(-pow(x / b, a)); - } else { - x = b * pow(-log(1.0 - p), 1.0 / a); - } - - switch (which) { - case 1: RETURN_DOUBLE(p); - case 2: RETURN_DOUBLE(x); - case 3: RETURN_DOUBLE(log(-log(1.0 - p)) / log(x / b)); - case 4: RETURN_DOUBLE(x / pow(-log(1.0 - p), 1.0 / a)); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -/* {{{ proto float stats_cdf_uniform(float par1, float par2, float par3, int which) - Not documented */ -PHP_FUNCTION(stats_cdf_uniform) -{ - double arg1; - double arg2; - double arg3; - double p; - double q; - double x; - double a; - double b; - long which; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dddl", &arg1, &arg2, &arg3, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 4) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Fourth parameter should be in the 1..4 range"); - RETURN_FALSE; - } - - if (which < 4) { - b = arg3; - } else { - a = arg3; - } - - if (which < 3) { - a = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - if (which == 1) { - p = 1 - exp(-pow(x / b, a)); - } else { - x = b * pow(-log(1.0 - p), 1.0 / a); - } - - switch (which) { - case 4: RETURN_DOUBLE((x - (1.0 - p) * a) / p); - case 3: RETURN_DOUBLE((x - p * b) / (1.0 - p)); - case 2: RETURN_DOUBLE(a + p * (b - a)); - case 1: - if (x < a) { - p = 0; - } else { - if (x > b) { - p = 1; - } else { - p = (x - a) / ( b - a); - } - } - RETURN_DOUBLE(p); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - -static double exponential_quantile(double p) -{ - return -log(1.0-p); -} - -static double exponential_cdf(double x) -{ - return (1.0 - exp(-x)); -} - -/* {{{ proto float stats_cdf_exponential(float par1, float par2, int which) - Not documented */ -PHP_FUNCTION(stats_cdf_exponential) -{ - double arg1; - double arg2; - double p; - double q; - double x; - double scale; - long which; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddl", &arg1, &arg2, &which) == FAILURE) { - RETURN_FALSE; - } - - if (which < 1 || which > 3) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Third parameter should be in the 1..3 range"); - RETURN_FALSE; - } - - if (which < 3) { - scale = arg2; - } else { - x = arg2; - } - - if (which == 1) { - x = arg1; - } else { - p = arg1; - q = 1.0 - p; - } - - switch (which) { - case 1: RETURN_DOUBLE(exponential_cdf(x / scale)); - case 2: RETURN_DOUBLE(scale * exponential_quantile(p)); - case 3: RETURN_DOUBLE(x / exponential_quantile(p)); - } - RETURN_FALSE; /* never here */ -} -/* }}} */ - - -/*********************/ -/* RANDLIB functions */ -/*********************/ - -/* {{{ proto void stats_rand_setall(int iseed1, int iseed2) - Not documented */ -PHP_FUNCTION(stats_rand_setall) -{ - long iseed_1; - long iseed_2; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ll", &iseed_1, &iseed_2) == FAILURE) { - RETURN_FALSE; - } - - setall(iseed_1, iseed_2); -} -/* }}} */ - -/* {{{ proto array stats_rand_get_seeds(void) - Not documented */ -PHP_FUNCTION(stats_rand_getsd) -{ - long iseed_1; - long iseed_2; - - if (ZEND_NUM_ARGS() != 0) { - WRONG_PARAM_COUNT; - } - getsd(&iseed_1, &iseed_2); - - array_init(return_value); - add_next_index_long(return_value, iseed_1); - add_next_index_long(return_value, iseed_2); -} -/* }}} */ - -/* {{{ proto int stats_rand_gen_iuniform(int low, int high) - Generates integer uniformly distributed between LOW (inclusive) and HIGH (inclusive) */ -PHP_FUNCTION(stats_rand_gen_iuniform) -{ - long low; - long high; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ll", &low, &high) == FAILURE) { - RETURN_FALSE; - } - - if (high - low > 2147483561L) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "high - low too large. low : %16ld high %16ld", low, high); - RETURN_FALSE; - } - if (low > high) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "low greater than high. low : %16ld high %16ld", low, high); - RETURN_FALSE; - } - - RETURN_LONG(ignuin(low, high)); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_funiform(float low, float high) - Generates uniform float between low (exclusive) and high (exclusive) */ -PHP_FUNCTION(stats_rand_gen_funiform) -{ - double low; - double high; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &low, &high) == FAILURE) { - RETURN_FALSE; - } - - if (low > high) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "low greater than high. low : %16.6E high : %16.6E", low, high); - RETURN_FALSE; - } - - RETURN_DOUBLE(genunf(low, high)); -} -/* }}} */ - -/* {{{ proto int stats_rand_gen_int(void) - Generates random integer between 1 and 2147483562 */ -PHP_FUNCTION(stats_rand_ignlgi) -{ - if (ZEND_NUM_ARGS() != 0) { - WRONG_PARAM_COUNT; - } - - RETURN_LONG(ignlgi()); -} -/* }}} */ - -/* {{{ proto float stats_rand_ranf(void) - Returns a random floating point number from a uniform distribution over 0 - 1 (endpoints of this interval are not returned) using the current generator */ -PHP_FUNCTION(stats_rand_ranf) -{ - if (ZEND_NUM_ARGS() != 0) { - WRONG_PARAM_COUNT; - } - - RETURN_DOUBLE(ranf()); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_beta(float a, float b) - Generates beta random deviate. Returns a random deviate from the beta distribution with parameters A and B. The density of the beta is x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 0), r - shape parameter of Gamma distribution (r > 0) */ -PHP_FUNCTION(stats_rand_gen_gamma) -{ - double a; - double r; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &a, &r) == FAILURE) { - RETURN_FALSE; - } - - if (!(a > 0.0 && r > 0.0)) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "A or R nonpositive. A value : %16.6E , R value : %16.6E", a, r); - RETURN_FALSE; - } - - RETURN_DOUBLE(gengam(a, r)); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_noncenral_chisquare(float df, float xnonc) - Generates random deviate from the distribution of a noncentral chisquare with "df" degrees of freedom and noncentrality parameter "xnonc". d must be >= 1.0, xnonc must >= 0.0 */ -PHP_FUNCTION(stats_rand_gen_noncentral_chisquare) -{ - double df; - double xnonc; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &df, &xnonc) == FAILURE) { - RETURN_FALSE; - } - - if (df < 1.0 || xnonc < 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "df < 1 or xnonc < 0. df value : %16.6E xnonc value : %16.6E", df, xnonc); - RETURN_FALSE; - } - - RETURN_DOUBLE(gennch(df, xnonc)); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_noncentral_f(float dfn, float dfd, float xnonc) - Generates a random deviate from the noncentral F (variance ratio) distribution with "dfn" degrees of freedom in the numerator, and "dfd" degrees of freedom in the denominator, and noncentrality parameter "xnonc". Method : directly generates ratio of noncentral numerator chisquare variate to central denominator chisquare variate. */ -PHP_FUNCTION(stats_rand_gen_noncenral_f) -{ - double dfn; - double dfd; - double xnonc; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &dfn, &dfd, &xnonc) == FAILURE) { - RETURN_FALSE; - } - - if (dfn < 1.0 || dfd <= 0.0 || xnonc < 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Either (1) Numerator df < 1.0 or (2) Denominator df <= 0.0 or (3) Noncentrality parameter < 0.0. dfn: %16.6E dfd: %16.6E xnonc: %16.6E", dfn, dfd, xnonc); - RETURN_FALSE; - } - - RETURN_DOUBLE(gennf(dfn, dfd, xnonc)); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_normal(float av, float sd) - Generates a single random deviate from a normal distribution with mean, av, and standard deviation, sd (sd >= 0). Method : Renames SNORM from TOMS as slightly modified by BWB to use RANF instead of SUNIF. */ -PHP_FUNCTION(stats_rand_gen_normal) -{ - double av; - double sd; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &av, &sd) == FAILURE) { - RETURN_FALSE; - } - - if (sd < 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "sd < 0.0 . sd : %16.6E", sd); - RETURN_FALSE; - } - - RETURN_DOUBLE(gennor(av, sd)); -} -/* }}} */ - -/* {{{ proto array stats_rand_phrase_to_seeds(string phrase) - Uses a phrase (characted string) to generate two seeds for the RGN random number generator. Trailing blanks are eliminated before the seeds are generated. Generated seed values will fall in the range 1..2^30. */ -PHP_FUNCTION(stats_rand_phrase_to_seeds) -{ - zval **par1; - char *arg1 = NULL; - long seed_1; - long seed_2; - - if (ZEND_NUM_ARGS() != 1 || zend_get_parameters_ex(1, &par1) == FAILURE) { - WRONG_PARAM_COUNT; - } - convert_to_string_ex(par1); - - arg1 = estrndup(Z_STRVAL_PP(par1), Z_STRLEN_PP(par1)); - phrtsd(arg1, &seed_1, &seed_2); - efree(arg1); - - array_init(return_value); - add_next_index_long(return_value, seed_1); - add_next_index_long(return_value, seed_2); -} -/* }}} */ - -/* {{{ proto int stats_rand_gen_ibinomial(int n, float pp) - Generates a single random deviate from a binomial distribution whose number of trials is "n" (n >= 0) and whose probability of an event in each trial is "pp" ([0;1]). Method : algorithm BTPE */ -PHP_FUNCTION(stats_rand_ibinomial) -{ - long n; - double pp; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ld", &n, &pp) == FAILURE) { - RETURN_FALSE; - } - - if ((n < 0) || (pp < 0.0) || (pp > 1.0)) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Bad values for the arguments. n : %ld pp : %16.6E", n, pp); - RETURN_FALSE; - } - - RETURN_LONG(ignbin(n, pp)); -} -/* }}} */ - -/* {{{ proto int stats_rand_gen_ibinomial_negative(int n, float p) - Generates a single random deviate from a negative binomial distribution. Arguments : n - the number of trials in the negative binomial distribution from which a random deviate is to be generated (n > 0), p - the probability of an event (0 < p < 1)). */ -PHP_FUNCTION(stats_rand_ibinomial_negative) -{ - long n; - double p; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ld", &n, &p) == FAILURE) { - RETURN_FALSE; - } - - if (n <= 0L) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "n < 0. n : %ld", n); - RETURN_FALSE; - } - if (p < 0.0F || p > 1.0F) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "p is out of range. p : %16.E", p); - RETURN_FALSE; - } - - RETURN_LONG(ignnbn(n, p)); -} -/* }}} */ - -/* {{{ proto int stats_rand_gen_ipoisson(float mu) - Generates a single random deviate from a Poisson distribution with mean "mu" (mu >= 0.0). */ -PHP_FUNCTION(stats_rand_gen_ipoisson) -{ - double mu; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "d", &mu) == FAILURE) { - RETURN_FALSE; - } - - if (mu < 0.0F) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "mu < 0.0 . mu : %16.6E", mu); - RETURN_FALSE; - } - - RETURN_LONG(ignpoi(mu)); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_noncentral_t(float df, float xnonc) - Generates a single random deviate from a noncentral T distribution. xnonc - noncentrality parameter. df must be >= 0.0*/ -PHP_FUNCTION(stats_rand_gen_noncentral_t) -{ - double df; - double xnonc; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &df, &xnonc) == FAILURE) { - RETURN_FALSE; - } - - if (df < 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "df <= 0 . df : %16.6E", df); - RETURN_FALSE; - } - - RETURN_DOUBLE(gennor(xnonc, 1) / sqrt(genchi(df) / df) ); -} -/* }}} */ - -/* {{{ proto float stats_rand_gen_t(float df) - Generates a single random deviate from a T distribution. df must be >= 0.0 */ -PHP_FUNCTION(stats_rand_gen_t) -{ - zval **arg1; - double df; - - if (ZEND_NUM_ARGS() != 1 || zend_get_parameters_ex(1, &arg1) == FAILURE) { - WRONG_PARAM_COUNT; - } - - convert_to_double_ex(arg1); - df = Z_DVAL_PP(arg1); - - if (df < 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "df <= 0 . df : %16.6E", df); - RETURN_FALSE; - } - - RETURN_DOUBLE(gennor(0, 1) / sqrt(genchi(df) / df)); -} -/* }}} */ - -/***************************/ -/* Start density functions */ -/***************************/ - -/* {{{ proto float stats_dens_normal(float x, float ave, float stdev) - Not documented */ -PHP_FUNCTION(stats_dens_normal) -{ - double stdev; - double ave; - double x; - double y; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddd", &x, &ave, &stdev) == FAILURE) { - RETURN_FALSE; - } - - if (stdev == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); - RETURN_FALSE; - } - - z = (x - ave) / stdev; - y = (1.0 / (stdev * sqrt(2.0 * STATS_PI))) * exp (-0.5 * z * z); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_cauchy(float x, float ave, float stdev) - Not documented */ -PHP_FUNCTION(stats_dens_cauchy) -{ - double stdev; - double ave; - double x; - double y; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddd", &x, &ave, &stdev) == FAILURE) { - RETURN_FALSE; - } - - if (stdev == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); - RETURN_FALSE; - } - - z = (x - ave) / stdev; - y = 1.0 / (stdev*STATS_PI * (1.0 + (z * z))); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_laplace(float x, float ave, float stdev) - Not documented */ -PHP_FUNCTION(stats_dens_laplace) -{ - double stdev; - double ave; - double x; - double y; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddd", &x, &ave, &stdev) == FAILURE) { - RETURN_FALSE; - } - - if (stdev == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); - RETURN_FALSE; - } - - z = fabs((x - ave) / stdev); - y = (1.0 / (2.0 * stdev)) * exp(- z); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_logistic(float x, float ave, float stdev) - Not documented */ -PHP_FUNCTION(stats_dens_logistic) -{ - double stdev; - double ave; - double x; - double y; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddd", &x, &ave, &stdev) == FAILURE) { - RETURN_FALSE; - } - - if (stdev == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "stdev is 0.0"); - RETURN_FALSE; - } - - z = exp((x - ave) / stdev); - y = z / (stdev * pow(1 + z, 2.0)); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_beta(float x, float a, float b) - Not documented */ -PHP_FUNCTION(stats_dens_beta) -{ - double a; - double b; - double beta; - double x; - double y; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &a, &b) == FAILURE) { - RETURN_FALSE; - } - - beta = 1.0 / exp(lgamma(a) + lgamma(b) - lgamma(a + b)); - y = beta * pow(x, a - 1.0) * pow(1.0 - x, b - 1.0); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_weibull(float x, float a, float b) - Not documented */ -PHP_FUNCTION(stats_dens_weibull) -{ - double a; - double b; - double x; - double y; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &a, &b) == FAILURE) { - RETURN_FALSE; - } - - if (b == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "b is 0.0"); - RETURN_FALSE; - } - - y = (a / b) * pow(x / b, a - 1.0) * exp(pow(- x / b, a)); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_den_uniform(float x, float a, float b) - Not documented */ -PHP_FUNCTION(stats_dens_uniform) -{ - double a; - double b; - double x; - double y; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddd", &x, &a, &b) == FAILURE) { - RETURN_FALSE; - } - - if (a == b) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "b == a == %16.6E", a); - RETURN_FALSE; - } - - if ((x <= b) && (x >= a)) { - y = 1.0 / (b - a); - } else { - y = 0.0; - } - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_chisquare(float x, float dfr) - Not documented */ -PHP_FUNCTION(stats_dens_chisquare) -{ - double dfr; - double e; - double x; - double y; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "dd", &x, &dfr) == FAILURE) { - RETURN_FALSE; - } - - e = dfr / 2.0; - z = ((e - 1.0) * log(x)) - ((x / 2.0) +(e * log(2.0)) + lgamma(e)); - y = exp (z); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_t(float x, float dfr) - Not documented */ -PHP_FUNCTION(stats_dens_t) -{ - double dfr; - double e; - double f; - double fac1; - double fac2; - double fac3; - double x; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &x, &dfr) == FAILURE) { - RETURN_FALSE; - } - - if (dfr == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "dfr == 0.0"); - RETURN_FALSE; - } - - e = dfr / 2.0; - f = e + 0.5; - fac1 = lgamma(f); - fac2 = f * log(1.0 + (x * x) / dfr); - fac3 = lgamma(e) + 0.5 * log(dfr * STATS_PI); - - RETURN_DOUBLE(exp( fac1 - (fac2 + fac3) )); -} -/* }}} */ - -/* {{{ proto float stats_dens_gamma(float x, float shape, float scale) - Not documented */ -PHP_FUNCTION(stats_dens_gamma) -{ - double shape; - double scale; - double x; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &shape, &scale) == FAILURE) { - RETURN_FALSE; - } - - if (scale == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "scale == 0.0"); - RETURN_FALSE; - } - - z = ((shape - 1.0) * log(x)) - - ( - (x / scale) + lgamma(shape) + (shape * log(scale)) - ) - ; - - RETURN_DOUBLE(exp(z)); -} -/* }}} */ - -/* {{{ proto float stats_dens_exponential(float x, float scale) - Not documented */ -PHP_FUNCTION(stats_dens_exponential) -{ - double scale; - double x; - double y; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &x, &scale) == FAILURE) { - RETURN_FALSE; - } - - if (scale == 0.0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "scale == 0.0"); - RETURN_FALSE; - } - - if (x < 0) { - y = 0; - } else { - y = exp(-x / scale) / scale; - } - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_f(float x, float dfr1, float dfr2) - */ -PHP_FUNCTION(stats_dens_f) -{ - double dfr1; - double dfr2; - double efr1; - double efr2; - double fac1; - double fac2; - double fac3; - double fac4; - double x; - double z; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &dfr1, &dfr2) == FAILURE) { - RETURN_FALSE; - } - - efr1 = dfr1 / 2.0; - efr2 = dfr2 / 2.0; - fac1 = (efr1 - 1.0) * log (x); - fac2 = (efr1 + efr2) * log (dfr2 + (dfr1 * x)); - fac3 = (efr1 * log (dfr1)) + (efr2 * log (dfr2)); - fac4 = lgamma (efr1) + lgamma (efr2) - lgamma (efr1 + efr2); - - z = (fac1 + fac3) - (fac2 + fac4); - - RETURN_DOUBLE(exp(z)); -} -/* }}} */ - -static double binom(double x, double n) -{ - int i; - double di; - double s = 1.0; - - for (i = 0; i < x; ++i) { - di = (double) i; - s = (s * (n - di)) / (di + 1.0); - } - - return s; -} - -/* {{{ proto float stats_dens_pmf_binomial(float x, float n, float pi) - Not documented */ -PHP_FUNCTION(stats_dens_pmf_binomial) -{ - double pi; - double y; - double n; - double x; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, - "ddd", &x, &n, &pi) == FAILURE) { - RETURN_FALSE; - } - - if ((x == 0.0 && n == 0.0) || (pi == 0.0 && x == 0.0) - || ( (1.0 - pi) == 0.0 && (n - x) == 0) ) { - - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Params leading to pow(0, 0). x:%16.6E n:%16.6E pi:%16.6E", x, n, pi); - RETURN_FALSE; - } - - y = binom(x,n) * pow(pi,x) * pow((1.0 - pi), (n - x)); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_pmf_poisson(float x, float lb) - Not documented */ -PHP_FUNCTION(stats_dens_pmf_poisson) -{ - double lb; - double z; - double x; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd", &x, &lb) == FAILURE) { - RETURN_FALSE; - } - - z = (x * log(lb)) - (lb + lgamma(x + 1.0)); - - RETURN_DOUBLE(exp(z)); -} -/* }}} */ - -/* {{{ proto float stats_dens_negative_binomial(float x, float n, float pi) - Not documented */ -PHP_FUNCTION(stats_dens_pmf_negative_binomial) -{ - double pi; - double y; - double n; - double x; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &n, &pi) == FAILURE) { - RETURN_FALSE; - } - - if ((pi == 0.0 && n == 0.0) || ((1.0 - pi) == 0.0 && x == 0.0)) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Params leading to pow(0, 0). x:%16.6E n:%16.6E pi:%16.6E", x, n, pi); - RETURN_FALSE; - } - - y = binom(x, n + x - 1.0) * pow(pi,n) * pow((1.0 - pi), x); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/* {{{ proto float stats_dens_pmf_hypergeometric(float n1, float n2, float N1, float N2) - */ -PHP_FUNCTION(stats_dens_pmf_hypergeometric) -{ - double y; - double N1; - double N2; - double n1; - double n2; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dddd", &n1, &n2, &N1, &N2) == FAILURE) { - RETURN_FALSE; - } - - if ((int)(n1 + n2) >= (int)(N1 + N2)) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "possible division by zero - n1+n2 >= N1+N2"); - /* RETURN_FALSE; */ - } - - y = binom(n1, N1) * binom (n2, N2)/binom(n1 + n2, N1 + N2); - - RETURN_DOUBLE(y); -} -/* }}} */ - -/************************/ -/* Statistics functions */ -/************************/ - -/* {{{ proto float stats_stat_powersum(array arr, float power) - Not documented */ -PHP_FUNCTION(stats_stat_powersum) -{ - zval **arg1, **arg2, **data; /* pointer to array entry */ - HashPosition pos; /* hash iterator */ - double power; - double sum = 0.0; - - if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { - WRONG_PARAM_COUNT; - } - - convert_to_array_ex(arg1); - convert_to_double_ex(arg2); - power = Z_DVAL_PP(arg2); - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data, &pos) == SUCCESS) { - convert_to_double_ex(data); - if (Z_DVAL_PP(data) != 0 && power != 0) { - sum += pow (Z_DVAL_PP(data), power); - } else { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Both value and power are zero"); - } - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos); - } - - RETURN_DOUBLE(sum); -} -/* }}} */ - -/* {{{ proto float stats_stat_innerproduct(array arr1, array arr2) - */ -PHP_FUNCTION(stats_stat_innerproduct) -{ - zval **arg1, **arg2; - zval **data1, **data2; /* pointers to array entries */ - HashPosition pos1; /* hash iterator */ - HashPosition pos2; /* hash iterator */ - double sum = 0.0; - - - if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { - WRONG_PARAM_COUNT; - } - convert_to_array_ex(arg1); - convert_to_array_ex(arg2); - - if (zend_hash_num_elements(Z_ARRVAL_PP(arg1)) != zend_hash_num_elements(Z_ARRVAL_PP(arg2))) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Unequal number of X and Y coordinates"); - RETURN_FALSE; - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS - && zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { - convert_to_double_ex(data1); - convert_to_double_ex(data2); - sum = Z_DVAL_PP(data1) * Z_DVAL_PP(data2); - - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); - } - - RETURN_DOUBLE(sum); -} -/* }}} */ - -/* {{{ proto float stats_stat_independent_t(array arr1, array arr2) - Not documented */ -PHP_FUNCTION(stats_stat_independent_t) -{ - zval **arg1, **arg2; - zval **data1, **data2; /* pointers to array entries */ - HashPosition pos1; /* hash iterator */ - HashPosition pos2; /* hash iterator */ - int xnum = 0, ynum = 0; - double cur; - double sx = 0.0; - double sxx = 0.0; - double sy = 0.0; - double syy = 0.0; - double mx; - double vx; - double my; - double vy; - double sp; - double fc; - double ts; - - if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { - WRONG_PARAM_COUNT; - } - convert_to_array_ex(arg1); - convert_to_array_ex(arg2); - - xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); - ynum = zend_hash_num_elements(Z_ARRVAL_PP(arg2)); - if ( xnum < 2 || ynum < 2) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Each argument should have more than 1 element"); - RETURN_FALSE; - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS) { - convert_to_double_ex(data1); - cur = Z_DVAL_PP(data1); - sx += cur; - sxx += cur * cur; - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { - convert_to_double_ex(data2); - cur = Z_DVAL_PP(data2); - sy += cur; - syy += cur * cur; - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); - } - - mx = sx / xnum; - my = sy / ynum; - vx = (sxx - (xnum * mx * mx)) / (xnum - 1.0); - vy = (syy - (ynum * my * my)) / (ynum - 1.0); - sp = (((xnum - 1.0) * vx) + ((ynum - 1.0) * vy)) / (xnum + ynum - 2.0); - fc = (1.0 / xnum) + (1.0 / ynum); - ts = (mx - my) / sqrt(sp * fc); - - RETURN_DOUBLE(ts); -} -/* }}} */ - -/* {{{ proto float stats_stat_paired_t(array arr1, array arr2) - Not documented */ -PHP_FUNCTION(stats_stat_paired_t) -{ - zval **arg1, **arg2, **data1, **data2; /* pointers to array entries */ - HashPosition pos1; /* hash iterator */ - HashPosition pos2; /* hash iterator */ - int xnum = 0; - int ynum = 0; - double sd = 0.0; - double sdd = 0.0; - double md; - double td; - double ts; - double cur; - - if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { - WRONG_PARAM_COUNT; - } - convert_to_array_ex(arg1); - convert_to_array_ex(arg2); - - xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); - ynum = zend_hash_num_elements(Z_ARRVAL_PP(arg2)); - - if (xnum != ynum) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Unequal number of X and Y coordinates"); - RETURN_FALSE; - } - if (xnum < 2) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "arr1 should have atleast 2 elements"); - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS - && - zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { - convert_to_double_ex(data1); - convert_to_double_ex(data2); - - cur = Z_DVAL_PP(data1) - Z_DVAL_PP(data2); - sd += cur; - sdd += cur * cur; - - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); - } - - md = sd / xnum; - td = sqrt((sdd - (xnum * md * md)) / (xnum - 1.0)); - ts = sqrt((double) xnum) * (md / td); - - RETURN_DOUBLE(ts); -} -/* }}} */ - -/* {{{ proto float stats_stat_percentile(float df, float xnonc) - Not documented */ -PHP_FUNCTION(stats_stat_percentile) -{ - zval **arg1, **arg2; - zval **data1; /* pointers to array entries */ - HashPosition pos1; /* hash iterator */ - long ilow; - long iupp; - int xnum = 0; - int cnt = -1; - double perc; - double low; - double upp; - double val = 0.0; - - if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { - WRONG_PARAM_COUNT; - } - - convert_to_array_ex(arg1); - convert_to_double_ex(arg2); - perc = Z_DVAL_PP(arg2); - - xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); - - if (zend_hash_sort(Z_ARRVAL_PP(arg1), zend_qsort, stats_array_data_compare, 1 TSRMLS_CC) == FAILURE) { - RETURN_FALSE; - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); - - low = .01 * perc * (double)xnum; - upp = .01 * (100.0 - perc) * (double)xnum; - ilow = floor(low); - iupp = floor(upp); - if ((ilow + iupp) == xnum) { - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS) { - if (++cnt == ilow - 1 ) { - convert_to_double_ex(data1); - val = Z_DVAL_PP(data1); - - zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1); - convert_to_double_ex(data1); - val += Z_DVAL_PP(data1); - val = val / 2.0; - break; - } - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); - } - } else { - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS) { - if (++cnt == ilow) { - convert_to_double_ex(data1); - val += Z_DVAL_PP(data1); - break; - } - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); - } - } - - RETURN_DOUBLE(val); -} -/* }}} */ - -/* {{{ proto float stats_stat_correlation(array arr1, array arr2) - Not documented */ -PHP_FUNCTION(stats_stat_correlation) -{ - zval **arg1, **arg2; - zval **data1, **data2; /* pointers to array entries */ - HashPosition pos1; /* hash iterator */ - HashPosition pos2; /* hash iterator */ - int xnum = 0; - int ynum = 0; - double sx = 0.0; - double sy = 0.0; - double sxx = 0.0; - double syy = 0.0; - double sxy = 0.0; - double mx; - double my; - double vx; - double vy; - double cc; - double rr; - - if (ZEND_NUM_ARGS() != 2 || zend_get_parameters_ex(2, &arg1, &arg2) == FAILURE) { - WRONG_PARAM_COUNT; - } - - convert_to_array_ex(arg1); - convert_to_array_ex(arg2); - - xnum = zend_hash_num_elements(Z_ARRVAL_PP(arg1)); - ynum = zend_hash_num_elements(Z_ARRVAL_PP(arg2)); - - if (xnum != ynum) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "Unequal number of X and Y coordinates"); - RETURN_FALSE; - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg1), &pos1); - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_PP(arg2), &pos2); - - while (zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg1), (void **)&data1, &pos1) == SUCCESS - && zend_hash_get_current_data_ex(Z_ARRVAL_PP(arg2), (void **)&data2, &pos2) == SUCCESS) { - - convert_to_double_ex(data1); - convert_to_double_ex(data2); - - sx += Z_DVAL_PP(data1); - sxx += Z_DVAL_PP(data1) * Z_DVAL_PP(data1); - sy += Z_DVAL_PP(data2); - syy += Z_DVAL_PP(data2) * Z_DVAL_PP(data2); - sxy += Z_DVAL_PP(data1) * Z_DVAL_PP(data2); - - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg1), &pos1); - zend_hash_move_forward_ex(Z_ARRVAL_PP(arg2), &pos2); - } - - mx = sx / xnum; - my = sy / ynum; - vx = sxx - (xnum * mx * mx); - vy = syy - (ynum * my * my); - cc = sxy - (xnum * mx * my); - rr = cc / sqrt(vx * vy); - - RETURN_DOUBLE(rr); -} -/* }}} */ - -/* {{{ proto float stats_stat_binomial_coef(int x, int n) - Not documented */ -PHP_FUNCTION(stats_stat_binomial_coef) -{ - int i; - int n; - int x; - double bc = 1.0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ll", &x, &n) == FAILURE) { - RETURN_FALSE; - } - - for (i = 0; i < x; ++i) { - bc = (bc * (n - i)) / (i + 1); - } - - RETURN_DOUBLE(bc); -} -/* }}} */ - -/* {{{ proto float stats_stat_gennch(int n) - Not documented */ -PHP_FUNCTION(stats_stat_factorial) -{ - int n; - int i; - double f = 1; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "l", &n) == FAILURE) { - RETURN_FALSE; - } - - for (i = 1; i <= n; ++i) { - f *= i; - } - - RETURN_DOUBLE(f); -} -/* }}} */ - - -/* {{{ php_population_variance -*/ -static long double php_math_mean(zval *arr) -{ - double sum = 0.0; - zval **entry; - HashPosition pos; - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { - convert_to_double_ex(entry); - sum += Z_DVAL_PP(entry); - zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); - } - /* - we don't check whether the array has 0 elements. this is left to the caller - no need - to kill performance by checking on every level. - */ - return sum / zend_hash_num_elements(Z_ARRVAL_P(arr)); -} -/* }}} */ - - -/* {{{ php_population_variance -*/ -static long double php_population_variance(zval *arr, zend_bool sample) -{ - double mean, vr = 0.0; - zval **entry; - HashPosition pos; - int elements_num; - - elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr)); - - mean = php_math_mean(arr); - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { - double d; - convert_to_double_ex(entry); - d = Z_DVAL_PP(entry) - mean; - vr += d*d; - zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); - } - if (sample) { - --elements_num; - } - return (vr / elements_num); -} -/* }}} */ - - -/* {{{ proto float stats_variance(array a [, bool sample]) - Returns the population variance */ -PHP_FUNCTION(stats_variance) -{ - zval *arr; - zend_bool sample = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a|b", &arr, &sample) == FAILURE) { - return; - } - if (zend_hash_num_elements(Z_ARRVAL_P(arr)) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); - RETURN_FALSE; - } - if (sample && zend_hash_num_elements(Z_ARRVAL_P(arr)) == 1) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has only 1 element"); - RETURN_FALSE; - } - RETURN_DOUBLE(php_population_variance(arr, sample)); -} -/* }}} */ - -/* {{{ proto float stats_standard_deviation(array a[, bool sample = false]) - Returns the standard deviation */ -PHP_FUNCTION(stats_standard_deviation) -{ - zval *arr; - zend_bool sample = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a|b", &arr, &sample) == FAILURE) { - return; - } - if (zend_hash_num_elements(Z_ARRVAL_P(arr)) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); - RETURN_FALSE; - } - if (sample && zend_hash_num_elements(Z_ARRVAL_P(arr)) == 1) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has only 1 element"); - RETURN_FALSE; - } - RETURN_DOUBLE(sqrt(php_population_variance(arr, sample))); -} -/* }}} */ - - -/* {{{ proto float stats_absolute_deviation(array a) - Returns the absolute deviation of an array of values*/ -PHP_FUNCTION(stats_absolute_deviation) -{ - zval *arr; - double mean = 0.0, abs_dev = 0.0; - zval **entry; - HashPosition pos; - int elements_num; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { - return; - } - if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); - RETURN_FALSE; - } - - mean = php_math_mean(arr); - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { - convert_to_double_ex(entry); - abs_dev += fabs(Z_DVAL_PP(entry) - mean); - zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); - } - - RETURN_DOUBLE(abs_dev / elements_num); -} -/* }}} */ - -/* {{{ proto float stats_harmonic_mean(array a) - Returns the harmonic mean of an array of values */ -PHP_FUNCTION(stats_harmonic_mean) -{ - zval *arr; - double sum = 0.0; - zval **entry; - HashPosition pos; - int elements_num; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { - return; - } - if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); - RETURN_FALSE; - } - - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { - convert_to_double_ex(entry); - if (Z_DVAL_PP(entry) == 0) { - RETURN_LONG(0); - } - sum += 1 / Z_DVAL_PP(entry); - zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); - } - - RETURN_DOUBLE(elements_num / sum); -} -/* }}} */ - -/* {{{ proto float stats_skew(array a) - Computes the skewness of the data in the array */ -PHP_FUNCTION(stats_skew) -{ - zval *arr; - double mean, std_dev, skew = 0.0; - zval **entry; - HashPosition pos; - int elements_num, i = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { - return; - } - if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); - RETURN_FALSE; - } - - mean = php_math_mean(arr); - std_dev = sqrt(php_population_variance(arr, 0)); - - /* the calculation of the skewness is protected of value "explosion". a bit more - FP operations performed but more accurateness. - */ - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { - double tmp; - convert_to_double_ex(entry); - tmp = ((Z_DVAL_PP(entry) - mean) / std_dev); - skew += (tmp*tmp*tmp - skew) / (i + 1); - zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); - ++i; - } - - RETURN_DOUBLE(skew); -} -/* }}} */ - -/* {{{ proto float stats_kurtosis(array a) - Computes the kurtosis of the data in the array */ -PHP_FUNCTION(stats_kurtosis) -{ - zval *arr; - double mean, std_dev, avg = 0.0; - zval **entry; - HashPosition pos; - int elements_num, i = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a", &arr) == FAILURE) { - return; - } - if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements"); - RETURN_FALSE; - } - - mean = php_math_mean(arr); - std_dev = sqrt(php_population_variance(arr, 0)); - - /* the calculation of the kurtosis is protected of value "explosion". a bit more - FP operations performed but more accurateness. - */ - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) { - double tmp; - convert_to_double_ex(entry); - tmp = ((Z_DVAL_PP(entry) - mean) / std_dev); - avg += (tmp*tmp*tmp*tmp - avg) / (i + 1); - - zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos); - ++i; - } - - RETURN_DOUBLE(avg - 3); -} -/* }}} */ - - -/* {{{ proto float stats_covariance(array a, array b) - Computes the covariance of two data sets */ -PHP_FUNCTION(stats_covariance) -{ - zval *arr_1, *arr_2; - double mean_1, mean_2, covar = 0.0; - zval **entry; - HashPosition pos_1, pos_2; - int elements_num, i = 0; - - if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "aa", &arr_1, &arr_2) == FAILURE) { - return; - } - if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr_1))) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The first array has zero elements"); - RETURN_FALSE; - } - if (zend_hash_num_elements(Z_ARRVAL_P(arr_2)) == 0) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The second array has zero elements"); - RETURN_FALSE; - } - if (elements_num != zend_hash_num_elements(Z_ARRVAL_P(arr_2))) { - php_error_docref(NULL TSRMLS_CC, E_WARNING, "The datasets are not of the same size"); - RETURN_FALSE; - } - - mean_1 = php_math_mean(arr_1); - mean_2 = php_math_mean(arr_2); - /* the calculation of the covariance is protected of value "explosion". a bit more - FP operations performed but more accurateness. - */ - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr_1), &pos_1); - zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr_2), &pos_2); - while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr_1), (void **)&entry, &pos_1) == SUCCESS) { - double tmp_1, tmp_2; - convert_to_double_ex(entry); - tmp_1 = Z_DVAL_PP(entry) - mean_1; - - if (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr_2), (void **)&entry, &pos_2) != SUCCESS) { - break; - } - convert_to_double_ex(entry); - tmp_2 = Z_DVAL_PP(entry) - mean_2; - - covar += (tmp_1 * tmp_2 - covar) / (i + 1); - - zend_hash_move_forward_ex(Z_ARRVAL_P(arr_1), &pos_1); - zend_hash_move_forward_ex(Z_ARRVAL_P(arr_2), &pos_2); - ++i; - } - - RETURN_DOUBLE(covar); -} -/* }}} */ - - -/* - * Local variables: - * tab-width: 4 - * c-basic-offset: 4 - * indent-tabs-mode: t - * End: - */ diff -dPNur stats-1.0.2/.svn/all-wcprops trunk/.svn/all-wcprops diff -dPNur stats-1.0.2/.svn/entries trunk/.svn/entries diff -dPNur stats-1.0.2/.svn/prop-base/cdflib.h.svn-base trunk/.svn/prop-base/cdflib.h.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/com.c.svn-base trunk/.svn/prop-base/com.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/config.m4.svn-base trunk/.svn/prop-base/config.m4.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/config.w32.svn-base trunk/.svn/prop-base/config.w32.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/CREDITS.svn-base trunk/.svn/prop-base/CREDITS.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/dcdflib.c.svn-base trunk/.svn/prop-base/dcdflib.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/fd_e_lgamma_r.c.svn-base trunk/.svn/prop-base/fd_e_lgamma_r.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/fd_e_log.c.svn-base trunk/.svn/prop-base/fd_e_log.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/fd_k_cos.c.svn-base trunk/.svn/prop-base/fd_k_cos.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/fd_k_sin.c.svn-base trunk/.svn/prop-base/fd_k_sin.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/fdlibm.h.svn-base trunk/.svn/prop-base/fdlibm.h.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/fd_w_lgamma.c.svn-base trunk/.svn/prop-base/fd_w_lgamma.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/ipmpar.c.svn-base trunk/.svn/prop-base/ipmpar.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/linpack.c.svn-base trunk/.svn/prop-base/linpack.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/package.xml.svn-base trunk/.svn/prop-base/package.xml.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/php_stats.c.svn-base trunk/.svn/prop-base/php_stats.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/php_stats.h.svn-base trunk/.svn/prop-base/php_stats.h.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/randlib.c.svn-base trunk/.svn/prop-base/randlib.c.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/randlib.h.svn-base trunk/.svn/prop-base/randlib.h.svn-base diff -dPNur stats-1.0.2/.svn/prop-base/TODO.svn-base trunk/.svn/prop-base/TODO.svn-base diff -dPNur stats-1.0.2/.svn/text-base/cdflib.h.svn-base trunk/.svn/text-base/cdflib.h.svn-base diff -dPNur stats-1.0.2/.svn/text-base/com.c.svn-base trunk/.svn/text-base/com.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/config.m4.svn-base trunk/.svn/text-base/config.m4.svn-base diff -dPNur stats-1.0.2/.svn/text-base/config.w32.svn-base trunk/.svn/text-base/config.w32.svn-base diff -dPNur stats-1.0.2/.svn/text-base/CREDITS.svn-base trunk/.svn/text-base/CREDITS.svn-base \ No newline at end of file diff -dPNur stats-1.0.2/.svn/text-base/dcdflib.c.svn-base trunk/.svn/text-base/dcdflib.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/fd_e_lgamma_r.c.svn-base trunk/.svn/text-base/fd_e_lgamma_r.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/fd_e_log.c.svn-base trunk/.svn/text-base/fd_e_log.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/fd_k_cos.c.svn-base trunk/.svn/text-base/fd_k_cos.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/fd_k_sin.c.svn-base trunk/.svn/text-base/fd_k_sin.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/fdlibm.h.svn-base trunk/.svn/text-base/fdlibm.h.svn-base diff -dPNur stats-1.0.2/.svn/text-base/fd_w_lgamma.c.svn-base trunk/.svn/text-base/fd_w_lgamma.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/ipmpar.c.svn-base trunk/.svn/text-base/ipmpar.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/linpack.c.svn-base trunk/.svn/text-base/linpack.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/package.xml.svn-base trunk/.svn/text-base/package.xml.svn-base diff -dPNur stats-1.0.2/.svn/text-base/php_stats.c.svn-base trunk/.svn/text-base/php_stats.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/php_stats.h.svn-base trunk/.svn/text-base/php_stats.h.svn-base diff -dPNur stats-1.0.2/.svn/text-base/randlib.c.svn-base trunk/.svn/text-base/randlib.c.svn-base diff -dPNur stats-1.0.2/.svn/text-base/randlib.h.svn-base trunk/.svn/text-base/randlib.h.svn-base diff -dPNur stats-1.0.2/.svn/text-base/TODO.svn-base trunk/.svn/text-base/TODO.svn-base \ No newline at end of file diff -dPNur stats-1.0.2/tests/.svn/all-wcprops trunk/tests/.svn/all-wcprops diff -dPNur stats-1.0.2/tests/.svn/entries trunk/tests/.svn/entries diff -dPNur stats-1.0.2/tests/.svn/prop-base/common.php.svn-base trunk/tests/.svn/prop-base/common.php.svn-base diff -dPNur stats-1.0.2/tests/.svn/prop-base/math_abs_dev.phpt.svn-base trunk/tests/.svn/prop-base/math_abs_dev.phpt.svn-base diff -dPNur stats-1.0.2/tests/.svn/prop-base/math_covariance.phpt.svn-base trunk/tests/.svn/prop-base/math_covariance.phpt.svn-base diff -dPNur stats-1.0.2/tests/.svn/prop-base/math_harmonic_mean.phpt.svn-base trunk/tests/.svn/prop-base/math_harmonic_mean.phpt.svn-base diff -dPNur stats-1.0.2/tests/.svn/prop-base/math_skew_kurtosis.phpt.svn-base trunk/tests/.svn/prop-base/math_skew_kurtosis.phpt.svn-base diff -dPNur stats-1.0.2/tests/.svn/prop-base/math_std_dev.phpt.svn-base trunk/tests/.svn/prop-base/math_std_dev.phpt.svn-base diff -dPNur stats-1.0.2/tests/.svn/text-base/common.php.svn-base trunk/tests/.svn/text-base/common.php.svn-base diff -dPNur stats-1.0.2/tests/.svn/text-base/math_abs_dev.phpt.svn-base trunk/tests/.svn/text-base/math_abs_dev.phpt.svn-base \ No newline at end of file diff -dPNur stats-1.0.2/tests/.svn/text-base/math_covariance.phpt.svn-base trunk/tests/.svn/text-base/math_covariance.phpt.svn-base \ No newline at end of file diff -dPNur stats-1.0.2/tests/.svn/text-base/math_harmonic_mean.phpt.svn-base trunk/tests/.svn/text-base/math_harmonic_mean.phpt.svn-base \ No newline at end of file diff -dPNur stats-1.0.2/tests/.svn/text-base/math_skew_kurtosis.phpt.svn-base trunk/tests/.svn/text-base/math_skew_kurtosis.phpt.svn-base \ No newline at end of file diff -dPNur stats-1.0.2/tests/.svn/text-base/math_std_dev.phpt.svn-base trunk/tests/.svn/text-base/math_std_dev.phpt.svn-base \ No newline at end of file