binomial.h

/**************************************************************************
 * Copyright (c) 2002-2011 T. M. Murali                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by   *
 * the Free Software Foundation, either version 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 * GNU General Public License for more details.                           *
 *                                                                        *
 * You should have received a copy of the GNU General Public License      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef _BINOMIAL_H
#define _BINOMIAL_H

#include <iostream>

#include <math.h>

#include "constants.h"

inline double log_fac(int n)
{
  // take care of special cases first.
  if (0 == n)
    return(0);
  if (1 == n)
    return(0);
  return(0.5*log(2 * PI * n) + n*log(n/E));
}

inline double fac(int n)
{
  return(exp(log_fac(n)));
}

inline double log_choose(int n, int k)
{
  return(log_fac(n) - log_fac(k) - log_fac(n - k));
}

inline double choose(int n, int k)
{
  return(exp(log_choose(n, k)));
}

inline double log_prob_binomial(int n, int k, double prob)
{
    if (n < k)
      {
#ifdef DEBUG
        cerr << "log_prob_binomial(): n = " << n << " < k = " << k << ", returning 0.\n";
#endif
        return(0);
      }
    if (0 > n )
      {
#ifdef DEBUG
        cerr << "log_prob_binomial(): n = " << n << " < 0, returning 0.\n";
#endif
        return(0);
      }
    if (1 < prob)
      {
#ifdef DEBUG
        cerr << "log_prob_binomial(): p = " << prob << " >= 1, returning 0.\n";
#endif
        return(0);
      }
    if (1 == prob)
      {
        if (n != k)
          {
#ifdef DEBUG
//            cerr << "log_prob_binomial(): p = " << prob << " == 1 and n = "
//            << n << " != k = " << k << ", returning 0.\n";
#endif
            return(0);
          }
        else
          {
#ifdef DEBUG
//            cerr << "log_prob_binomial(): p = " << prob << " == 1 and n = "
//            << n << " == k = " << k << ", returning 1.\n";
#endif
            return(1);
          }
      }
    if (0 > prob)
         {
#ifdef DEBUG
           cerr << "log_prob_binomial(): p = " << prob << " <= 0, returning 0.\n";
#endif
           return(0);
         }
    if (0 == prob)
      {
        if (0 != k)
          {
#ifdef DEBUG
//             cerr << "log_prob_binomial(): p = " << prob << " == 1 and k = " << k
//                  << " != 0, returning 0.\n";
#endif
            return(0);
          }
        else
          {
#ifdef DEBUG
//             cerr << "log_prob_binomial(): p = " << prob << " == 1 and k = " << k
//                  << " == 0, returning 1.\n";
#endif
            return(1);
          }
      }
    return(log_choose(n, k) + k*log(prob) + (n - k)*log(1 - prob));
  }

inline double sum_binomial_probs(int n, int start, int end, double prob)
  {
    double sum = 0;
    for (int i = start; i <= end; i++)
      {
        sum += exp(log_prob_binomial(n, i, prob));
      }
    return(sum);
  }

// there are total items. of these, specialSize items are special.
// pick a sample of size sampleSize. what is the probability that at
// specialSizeInSample in the sample are special?
inline double sum_binomial_product(int total, int specialSize, int sampleSize,
                            int specialSizeInSample)
{
  int i;
  double sum = 0;
  for (i = specialSizeInSample; i < sampleSize; i++)
    // pick i special and sampleSize - i non-special
    sum += choose(specialSize, i)*choose(total - specialSize, sampleSize - i);
  // there are choose(total, sampleSize) possible choices.
  sum = sum/choose(total, sampleSize);
  return(sum);
}

#endif // _BINOMIAL_H