1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
#include "Dirichlet.h"
#include "Sum.h"
#include "Product.h"
#include <iostream>
long double dirichlet(const vector<long double>& probs,
const vector<int>& obs,
long double s) {
vector<long double> alphas;
for (vector<int>::const_iterator o = obs.begin(); o != obs.end(); ++o)
alphas.push_back(*o + 1 * s);
vector<long double> obsProbs;
vector<long double>::const_iterator a = alphas.begin();
vector<long double>::const_iterator p = probs.begin();
for (; p != probs.end() && a != alphas.end(); ++p, ++a) {
obsProbs.push_back(pow(*p, *a - 1));
}
return 1.0 / beta(alphas) * product(obsProbs);
}
long double dirichletMaximumLikelihoodRatio(const vector<long double>& probs,
const vector<int>& obs,
long double s) {
long double maximizingObs = obs.size() / sum(obs);
vector<int> m(obs.size(), maximizingObs);
return dirichlet(probs, obs, s) / dirichlet(probs, m, s);
}
// XXX the logspace versions are broken
long double dirichletln(const vector<long double>& probs,
const vector<int>& obs,
long double s) {
vector<long double> alphas;
for (vector<int>::const_iterator o = obs.begin(); o != obs.end(); ++o)
alphas.push_back(*o + 1 * s);
vector<long double> obsProbs;
vector<long double>::const_iterator a = alphas.begin();
vector<long double>::const_iterator p = probs.begin();
for (; p != probs.end() && a != alphas.end(); ++p, ++a) {
obsProbs.push_back(powln(log(*p), *a - 1));
}
return log(1.0) - (betaln(alphas) + sum(obsProbs));
}
long double dirichletMaximumLikelihoodRatioln(const vector<long double>& probs,
const vector<int>& obs,
long double s) {
long double maximizingObs = (long double) obs.size() / (long double) sum(obs);
vector<int> m(obs.size(), maximizingObs);
return dirichletln(probs, obs, s) - dirichletln(probs, m, s);
}
|