https://github.com/ekg/freebayes
Tip revision: 60850dc518fc453622cbb40ad6dd9f67644ed859 authored by Pjotr Prins on 16 December 2020, 10:18:06 UTC
Disable cmake, but leave the file with a message
Disable cmake, but leave the file with a message
Tip revision: 60850dc
Multinomial.cpp
#include "Multinomial.h"
#include "Sum.h"
#include "Product.h"
long double multinomialSamplingProb(const vector<long double>& probs, const vector<int>& obs) {
vector<long double> factorials;
vector<long double> probsPowObs;
factorials.resize(obs.size());
transform(obs.begin(), obs.end(), factorials.begin(), factorial);
vector<long double>::const_iterator p = probs.begin();
vector<int>::const_iterator o = obs.begin();
for (; p != probs.end() && o != obs.end(); ++p, ++o) {
probsPowObs.push_back(pow(*p, *o));
}
return factorial(sum(obs)) / product(factorials) * product(probsPowObs);
}
// TODO rename to reflect the fact that this is the multinomial sampling
// probability for obs counts given probs probabilities
long double multinomialSamplingProbLn(const vector<long double>& probs, const vector<int>& obs) {
vector<long double> factorials;
vector<long double> probsPowObs;
factorials.resize(obs.size());
transform(obs.begin(), obs.end(), factorials.begin(), factorialln);
vector<long double>::const_iterator p = probs.begin();
vector<int>::const_iterator o = obs.begin();
for (; p != probs.end() && o != obs.end(); ++p, ++o) {
probsPowObs.push_back(powln(log(*p), *o));
}
return factorialln(sum(obs)) - sum(factorials) + sum(probsPowObs);
}
long double multinomialCoefficientLn(int n, const vector<int>& counts) {
vector<long double> count_factorials;
count_factorials.resize(counts.size());
transform(counts.begin(), counts.end(), count_factorials.begin(), factorialln);
return factorialln(n) - sum(count_factorials);
}
long double samplingProbLn(const vector<long double>& probs, const vector<int>& obs) {
vector<long double>::const_iterator p = probs.begin();
vector<int>::const_iterator o = obs.begin();
long double r = 0;
for (; p != probs.end() && o != obs.end(); ++p, ++o) {
r += powln(log(*p), *o);
}
return r;
}