1#!/usr/bin/env python3
2
3from scipy.special import gamma, gammaln
4import operator
5import math
6from logsumexp import logsumexp
7from factorialln import factorialln
8
9def fold(func, iterable, initial=None, reverse=False):
10    x=initial
11    if reverse:
12        iterable=reversed(iterable)
13    for e in iterable:
14        x=func(x,e) if x is not None else e
15    return x
16
17def product(listy):
18    return fold(operator.mul, listy)
19
20def beta(alphas):
21    """Multivariate beta function"""
22    return math.exp(sum(map(gammaln, alphas)) - gammaln(sum(alphas)))
23
24def dirichlet(probs, obs, s=1):
25    """Dirichlet probability density for a given set probabilities and prior
26    observation counts.  s serves as a concentration parameter between 0 and 1,
27    with smaller s yielding progressively more diffuse probability density
28    distributions."""
29    alphas = [(a + 1) * float(s) for a in obs]
30    return 1 / beta(alphas) * product([math.pow(x, a - 1) for x, a in zip(probs, alphas)])
31
32def dirichlet_maximum_likelihood_ratio(probs, obs, s=1):
33    """Ratio between the dirichlet for the specific probs and obs and the
34    maximum likelihood value for the dirichlet over the given probs (this can
35    be determined by partitioning the observations according to the
36    probabilities)"""
37    maximum_likelihood = dirichlet(probs, [float(sum(obs)) / len(obs) for o in obs], s)
38    return dirichlet(probs, obs, s) / float(maximum_likelihood)
39
40def multinomial(probs, obs):
41    return math.factorial(sum(obs)) / product(list(map(math.factorial, obs))) * product([math.pow(p, x) for p,x in zip(probs, obs)])
42
43def multinomialln(probs, obs):
44    return factorialln(sum(obs)) - sum(map(factorialln, obs)) + sum([math.log(math.pow(p, x)) for p,x in zip(probs, obs)])
45
46def multinomial_coefficientln(n, counts):
47    return factorialln(n) - sum(map(factorialln, counts))
48
49def multinomial_coefficient(n, counts):
50    return math.exp(multinomial_coefficientln(n, counts))
51
52def multinomial_dirichlet(probs, obs): return multinomial(probs, obs) * dirichlet(probs, obs)
53
54# NOTE:
55# I started exploring the multinomial_maximum_likelihood_ratio to see if the
56# same maximim likelihood estimation approach could be applied to multinomials.
57# It *can't* for the obvious reason that you can't have non-integral
58# observation counts while the dirichlet distribution is defined across
59# non-integral alphas!  I see no clean way to resolve this using the
60# multinomial distribution and now have a better understanding of the use and
61# abuse of the dirichlet distribution as a conjugate prior for multinomial
62# posteriors.
63def multinomial_maximum_likelihood_ratio(probs, obs):
64    maximum_likelihood = multinomial(probs, [float(sum(obs)) / len(obs) for o in obs])
65    return multinomial(probs, obs) / float(maximum_likelihood)
66
67def binomial(successes, trials, prob):
68    return math.factorial(trials) / (math.factorial(successes) * math.factorial(trials - successes)) * math.pow(prob, successes) * math.pow(1 - prob, trials - successes)
69