get_k.py

import numpy as np


def p_max(k, f):

    rnd = int(round(1. / f))
    sum = 0

    for n in range(1, rnd + 1):
        sum += (-1) ** (n - 1) * (float(rnd) / n) * (1 - n * f) ** k

    return sum


def bestk(f, P, guesslow, guesshigh):

    ks = np.array(range(guesslow, guesshigh + 1))
    Plist = []

    for k in ks:
        Plist.append(p_max(k, f))

    Plist = np.array(Plist)

    Plist_neg_inds, = np.where(Plist < P)
    bestk = ks[Plist_neg_inds[0]]

    print bestk


def bestf(k, P):

    fs = np.logspace(-3., -1.)
    Plist = []

    for f in fs:
        Plist.append(p_max(k, f))

    Plist = np.array(Plist)

    Plist_neg_inds, = np.where(Plist < P)
    bestf = fs[Plist_neg_inds[0]]

    print bestf


if __name__ == "__main__":

    import sys

    if not (len(sys.argv) == 4 or len(sys.argv) == 6):

        print "usage: python get_k.py <f | k (which are you searching)> \
                                 \n\t <f | k (which you are given, float or int respectively)> \
                                 \n\t <P> \
                                 \n\t <kguess_low> <kguess_high> (if searching ks)"
        assert(0)

    if sys.argv[1] == "f":
        bestf(float(sys.argv[2]), float(sys.argv[3]))
    elif sys.argv[1] == "k":
        bestk(float(sys.argv[2]), float(sys.argv[3]),
              int(sys.argv[4]), int(sys.argv[5]))

code for Gavin's geology research, finds the minimum number of samples needed for an analysis to achieve a specified confidence level for the age fraction of a population having been found or, alternatively, for finding the lowest population percentages you can find with given confidence, if given a number of samples