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