diff options
Diffstat (limited to 'lib/binomial_pdf.c')
| -rw-r--r-- | lib/binomial_pdf.c | 63 |
1 files changed, 63 insertions, 0 deletions
diff --git a/lib/binomial_pdf.c b/lib/binomial_pdf.c new file mode 100644 index 0000000..8b3a442 --- /dev/null +++ b/lib/binomial_pdf.c @@ -0,0 +1,63 @@ +// Binomial PDF +// +// binomial_pdf - binomial probability density +// +// P(X=k) = C(n, k) * p^k * (1-p)^(n-k) +// +// Inputs: 1 (k values) +// Outputs: 1 (PDF value) +// Constants: 1 (n trials, p probability) + +#include <math.h> +#include "ts.h" + +typedef struct { + uint_t n_trials; + fp_t p; +} ts_binomial_pdf_params_t; + +ts_parameter_t const ts_binomial_pdf_p = { .n_inputs = 1, .n_outputs = 1, .n_constants = 1, .params_size = sizeof(ts_binomial_pdf_params_t) }; + +uint_t ts_binomial_pdf_offset( ts_constant_t const *c ) +{ + ( void )c; + return ( 0 ); +} + +int ts_binomial_pdf( + fp_t * const *y, uint_t y_n, + fp_t const * const *x, uint_t x_n, + uint_t n, ts_constant_t const *c +) +{ + TS_PROLOGUE( ts_binomial_pdf_p ); + if ( !c ) return -1; + uint_t n_trials = ((ts_binomial_pdf_params_t *)c)->n_trials; + fp_t p = ((ts_binomial_pdf_params_t *)c)->p; + if (p < 0.0 || p > 1.0 || isnan(p) || isinf(p)) return -1; + for( uint_t i = 0; i < n; ++i ) + { + uint_t k = (uint_t)x[0][i]; + if( k > n_trials ) + { + y[0][i] = 0.0; + } + else if( p == 0.0 ) + { + y[0][i] = k == 0 ? 1.0 : 0.0; + } + else if( p == 1.0 ) + { + y[0][i] = k == n_trials ? 1.0 : 0.0; + } + else + { + fp_t log_comb = lgamma( n_trials + 1.0 ) - lgamma( k + 1.0 ) - lgamma( n_trials - k + 1.0 ); + y[0][i] = exp( log_comb + k * log( p ) + ( n_trials - k ) * log( 1.0 - p ) ); + } + } + return 0; +} + + + |