#! /usr/bin/python
# -*- coding: utf-8 -*-
#
# Copyright 2016 Google Inc. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""One-sided beta-binomial test.

Performs a test about the null hypothesis that the observed number of
successes in a given number of trials was generated by a beta-binomial
distribution with fixed and known parameters alpha and beta.

Computes the p-value of seeing at least as many successes as were actually
observed (one-sided, upper-tail alternative).

Prints verbose test results.

"""

from __future__ import division, print_function

import math
import sys


def LogBetaBinomPMF(k, n, alpha, beta):
  lp = math.lgamma(n + 1)
  lp -= math.lgamma(k + 1) + math.lgamma(n - k + 1)
  lp += math.lgamma(k + alpha) + math.lgamma(n - k + beta)
  lp -= math.lgamma(n + alpha + beta)
  lp += math.lgamma(alpha + beta)
  lp -= math.lgamma(alpha) + math.lgamma(beta)
  return lp


def UpperTailPValue(successes, trials, alpha, beta):
  assert 0 <= successes <= trials
  assert alpha >= 0 and beta >= 0
  p_value = 0
  for k in range(successes, trials + 1):
    p_value += math.exp(LogBetaBinomPMF(k, trials, alpha, beta))
  return p_value


def UpperTailTest(x, n, alpha, beta, conf_level=0.95):
  p_value = UpperTailPValue(x, n, alpha, beta)
  print()
  print('\tExact beta-binomial test')
  print()
  print('data:  %d and %d' % (x, n))
  print('number of successes = %d, number of trials = %d, p-value = %.7g'
        % (x, n, p_value))
  print('model of null hypothesis: successes ~ betabinom(%d; %g, %g)'
        % (n, alpha, beta))
  print()
  return p_value


def main(argv):
  successes = int(argv[1])
  trials = int(argv[2])
  alpha = float(argv[3])
  beta = float(argv[4])
  conf_level = 0.95 if len(argv) < 6 else float(argv[5])
  p_value = UpperTailTest(successes, trials, alpha, beta, conf_level)
  significance_level = 1 - conf_level
  sys.exit(0 if p_value > significance_level else 1)
  return


if __name__ == '__main__':
  main(sys.argv)

#! /usr/bin/python
# -*- coding: utf-8 -*-
#
# Copyright 2016 Google Inc. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""One-sided beta-binomial test.

Performs a test about the null hypothesis that the observed number of
successes in a given number of trials was generated by a beta-binomial
distribution with fixed and known parameters alpha and beta.

Computes the p-value of seeing at least as many successes as were actually
observed (one-sided, upper-tail alternative).

Prints verbose test results.

"""

from __future__ import division, print_function

import math
import sys


def LogBetaBinomPMF(k, n, alpha, beta):
  lp = math.lgamma(n + 1)
  lp -= math.lgamma(k + 1) + math.lgamma(n - k + 1)
  lp += math.lgamma(k + alpha) + math.lgamma(n - k + beta)
  lp -= math.lgamma(n + alpha + beta)
  lp += math.lgamma(alpha + beta)
  lp -= math.lgamma(alpha) + math.lgamma(beta)
  return lp


def UpperTailPValue(successes, trials, alpha, beta):
  assert 0 <= successes <= trials
  assert alpha >= 0 and beta >= 0
  # This evaluates the cumulative distribution function in the upper tail by
  # naive summation of the probability mass function. This is linear in the
  # length of the tail. To make it a constant-time computation, we would need a
  # way to evaluate the generalized hypergeometric function 3F2 at unity.
  p_value = 0
  for k in range(successes, trials + 1):
    p_value += math.exp(LogBetaBinomPMF(k, trials, alpha, beta))
  return p_value


def UpperTailTest(x, n, alpha, beta, conf_level=0.95):
  p_value = UpperTailPValue(x, n, alpha, beta)
  print()
  print('\tExact beta-binomial test')
  print()
  print('data:  %d and %d' % (x, n))
  print('number of successes = %d, number of trials = %d, p-value = %.7g'
        % (x, n, p_value))
  print('model of null hypothesis: successes ~ betabinom(%d; %g, %g)'
        % (n, alpha, beta))
  print()
  return p_value


def main(argv):
  successes = int(argv[1])
  trials = int(argv[2])
  alpha = float(argv[3])
  beta = float(argv[4])
  conf_level = 0.95 if len(argv) < 6 else float(argv[5])
  p_value = UpperTailTest(successes, trials, alpha, beta, conf_level)
  significance_level = 1 - conf_level
  sys.exit(0 if p_value > significance_level else 1)
  return


if __name__ == '__main__':
  main(sys.argv)

Clarify how CDF is computed.