Wilson Binomial Proportion Confidence Interval

wilson_interval.h

#ifndef WILSON_INTERVAL_H                      /* Prevent Multiple Inclusion */
#define WILSON_INTERVAL_H
/* ========================================================================= **
**              WILSON BINOMIAL PROPORTION CONFIDENCE INTERVAL               **
** ========================================================================= **
** Copyright (C) Jonah H. Harris <jonah.harris@gmail.com>                    **
** All Rights Reserved.                                                      **
**                                                                           **
** Permission is hereby granted, free of charge, to any person obtaining a   **
** copy of this software and associated documentation files (the "Software") **
** to deal in the Software without restriction, including without limitation **
** the rights to use, copy, modify, merge, publish, distribute, sublicense,  **
** and/or sell copies of the Software, and to permit persons to whom the     **
** Software is furnished to do so, subject to the following conditions:      **
**                                                                           **
** The above copyright notice and this permission notice shall be included   **
** in all copies or substantial portions of the Software.                    **
**                                                                           **
** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS   **
** OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF                **
** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN **
** NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,  **
** DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR     **
** OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **
** USE OR OTHER DEALINGS IN THE SOFTWARE.                                    **
** ========================================================================= */

/* ========================================================================= */
/* -- INCLUSIONS ----------------------------------------------------------- */
/* ========================================================================= */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>

/* ========================================================================= */
/* -- INLINE FUNCTIONS ----------------------------------------------------- */
/* ========================================================================= */

static inline double
wilson_pnormdist (
  double qn
) {
  static double b[11] = {
    1.570796288,       0.03706987906,    -0.8364353589e-3,
    -0.2250947176e-3,  0.6841218299e-5,   0.5824238515e-5,
    -0.104527497e-5,   0.8360937017e-7,  -0.3231081277e-8,
    0.3657763036e-10,  0.6936233982e-12
  };
  double w1, w3;
  unsigned int i;

  if ((qn < 0) || (qn > 1)) {
    return (0);
  }
  if (qn == 0.5) {
    return (0);
  }

  w1 = qn;
  if (qn > 0.5) {
    w1 = 1.0 - w1;
  }
  w3 = -log(4.0 * w1 * (1.0 - w1));
  w1 = b[0];
  for (i = 1; i < 11; i++) {
    w1 = w1 + (b[i] * pow(w3, i));
  }

  if (qn > 0.5) {
    return (sqrt(w1 * w3));
  } else {
    return (-sqrt(w1 * w3));
  }
} /* wilson_pnormdist() */

/* ------------------------------------------------------------------------- */

static inline int
wilson_interval (
  double  k,
  double  n,
  double  confidence,
  bool    use_continuity_correction,
  double *lower_bound,
  double *upper_bound
) {
#define WSDBLAEQ(x, y)    (fabs(x - y) < FLT_EPSILON)
#define WSMAX(X, Y)       (((X) > (Y)) ? (X) : (Y))
#define WSMIN(X, Y)       (((X) < (Y)) ? (X) : (Y))
#define WSSET_IF_NOT_NULL(dstp, src)                                          \
  do {                                                                        \
    if (NULL != dstp) {                                                       \
      *dstp = src;                                                            \
    }                                                                         \
  } while (0)

  if (0 == n) {
    return (EXIT_FAILURE);
  }

  double phat = (k / n);
  double z = 1.959963984540;
  if (!(WSDBLAEQ(confidence, 0.95))) {
    z = wilson_pnormdist(1.0 - ((1.0 - confidence) / 2.0));
  }
  double z2 = (z * z);

  if (use_continuity_correction) {
    double a = 2 * (n + z2);
    double b = 2 * n * phat + z2;
    double c = z *
      sqrt(z2 - 1.0 / n + 4 * n * phat * (1 - phat) + (4 * phat - 2)) + 1;
    double d = z *
      sqrt(z2 - 1.0 / n + 4 * n * phat * (1 - phat) - (4 * phat - 2)) + 1;
    double lower = (WSDBLAEQ(0.0, phat)) ? 0.0 : WSMAX(0.0, (b - c) / a);
    double upper = (WSDBLAEQ(1.0, phat)) ? 1.0 : WSMIN(1.0, (b + d) / a);
    WSSET_IF_NOT_NULL(lower_bound, lower);
    WSSET_IF_NOT_NULL(upper_bound, upper);
  } else {
    double a = 1 + z2 / n;
    double b = phat + z2 / (2 * n);
    double c = z * sqrt((phat * (1 - phat) + z2 / (4 * n)) / n);
    WSSET_IF_NOT_NULL(lower_bound, ((b - c) / a));
    WSSET_IF_NOT_NULL(upper_bound, ((b + c) / a));
  }

  return (EXIT_SUCCESS);
#undef WSSET_IF_NOT_NULL
#undef WSMIN
#undef WSMAX
#undef WSDBLAEQ
} /* wilson_interval() */

#endif /* WILSON_INTERVAL_H */

wilson_test.c

#include "wilson_interval.h"
#include <sys/time.h>

int
main (
  int   argc,
  char *argv[]
) {
  ((void) argc);
  ((void) argv);

  srand(time(NULL));
  for (int ii = 0; ii < 100; ++ii) {
    int k = 0;
    int n = 0;

    int x = ((int) rand() % (int) 10e2) + 1;
    int y = ((int) rand() % (int) 10e2) + 1;
    if (x > y) {
      n = x;
      k = y;
    } else {
      n = y;
      k = x;
    }

    double lower = 0.0;
    double upper = 0.0;
    if (wilson_interval(k, n, 0.95, false, &lower, &upper)) {
      printf("yikes!\n");
      exit(EXIT_FAILURE);
    } else {
      printf("K: %d\tN: %d\tAverage: %.2f%% Wilson: %.2f%% [%.2f%%, %.2f%%]\n",
        k, n,
        (1.0 * k/n * 100),
        (((lower + upper) / 2) * 100),
        (lower * 100), (upper * 100));
    }
  }
}