/* * Adapted from R code published in conjunction with: * * Liao, J.G. And Rosen, O. (2001) Fast and Stable Algorithms for * Computing and Sampling from the Noncentral Hypergeometric * Distribution. The American Statistician 55, 366-369. */ #include #include "DHyper.h" #include #include #include #include #include #include using std::vector; using std::max; using std::min; using std::vector; using std::logic_error; namespace jags { namespace bugs { DHyper::DHyper() : RScalarDist("dhyper", 4, DIST_SPECIAL, true) {} bool DHyper::canBound() const { return false; } static void getParameters(int &n1, int &n2, int &m1, double &psi, vector const ¶meters) { n1 = static_cast(*parameters[0]); n2 = static_cast(*parameters[1]); m1 = static_cast(*parameters[2]); psi = *parameters[3]; } bool DHyper::checkParameterDiscrete (vector const &mask) const { // Check that n1, n2, m1 are discrete-valued for (unsigned int i = 0; i < 3; ++i) { if (mask[i] == false) return false; } return true; } bool DHyper::checkParameterValue (vector const ¶ms) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, params); if (n1 < 0 || n2 < 0) return false; else if (m1 < 0 || m1 > n1 + n2) return false; else if (psi <= 0) return false; else return true; } /* Calculates the mode of the hypergeometric distribution We solve the equation p(x) = p(x-1) for continuous x and then take the floor of x. This reduces to solving a quadratic equation in x, by the recurrence relation embodied in the function "rfunction", i.e. p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x) */ static int modeCompute(int n1, int n2, int m1, double psi) { double a = psi - 1; double b = -((m1 + n1 + 2) * psi + n2 - m1); double c = psi * (n1 + 1) * (m1 + 1); double q = b; if (b > 0) { q += sqrt(b * b - 4 * a * c); } else { q -= sqrt(b * b - 4 * a * c); } q = -q/2; int mode = static_cast(c/q); if (mode >= 0 && mode >= m1 - n2 && mode <= n1 && mode <= m1) { return mode; } else { return static_cast(q/a); } } /* The recurrence relation p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x) allows us to calculate the hypergeometric probabilities recursively, avoiding combinatoric problems. */ double rfunction(int n1, int n2, int m1, double psi, int i) { return psi * (n1 - i + 1) * (m1 - i + 1)/(i * (n2 - m1 + i)); } /** * Returns a vector of normalized probabilities for a hypergeometric * distribution. If the returned vector p is of length N, then p[0] * represents the probability of the lowest possible value ll = max(0, * m1 - n2) and p[N - 1] represents the probability of the largest * possible value uu = min(n1, m1) */ static vector density(int n1, int n2, int m1, double psi) { int ll = max(0U, m1 - n2); int uu = min(n1, m1); int N = uu - ll + 1; vector p(N); // Density at mode has reference value 1 int mode = modeCompute(n1, n2, m1, psi); p[mode - ll] = 1; // Calculate relative density above the mode if (mode < uu) { double r = 1; for (int i = mode + 1; i <= uu; ++i) { r *= rfunction(n1, n2, m1, psi, i); p[i - ll] = r; } } // Calculate relative density below the node if (mode > ll) { double r = 1; for (int i = mode - 1; i >= ll; --i) { r /= rfunction(n1, n2, m1, psi, i + 1); p[i - ll] = r; } } //Normalize double sump = 0; for (int i = 0; i < N; ++i) { sump += p[i]; } for (int i = 0; i < N; ++i) { p[i] /= sump; } return p; } /* * Sample from a unimodal discrete distribution given a vector of normalized * probabilities pi[] and an estimate of the mode. U should be drawn from a * uniform(0,1) distribution. */ static int sampleWithMode(int mode, vector const &pi, double U) { //We try to spend our probability allocation U as quickly as //possible, starting with the mode. U -= pi[mode]; if (U <= 0) return mode; //Then we work our way outwards from the mode, one step at a time. int lower = mode - 1; int upper = mode + 1; while (true) { if (lower < 0) { for ( ; upper < pi.size() - 1; ++upper) { U -= p[upper]; if (U <= 0) break; } return upper; } else if (upper >= pi.size()) { for ( ; lower > 0; --lower) { U -= p[lower]; if (U <= 0) break; } return lower; } else if (p[upper] > p[lower]) { U -= p[upper]; if (U <= 0) return upper; else ++upper; } else { U -= p[lower]; if (U <= 0) return lower; else --lower; } } return mode; //-Wall } double DHyper::d(double z, PDFType type, vector const ¶meters, bool give_log) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, parameters); int x = static_cast(z); int ll = max(0U, m1 - n2); int uu = min(n1, m1); double den = 0; if (x >= ll && x <= uu) { den = density(n1, n2, m1, psi)[x - ll]; } if (give_log) { return den == 0 ? JAGS_NEGINF : log(den); } else { return den; } } double DHyper::p(double x, vector const ¶meters, bool lower, bool give_log) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, parameters); int ll = max((int) 0, m1 - n2); int uu = min(n1, m1); double sumpi = 0; if (x >= ll) { if (x >= uu) { sumpi = 1; } else { vector pi = density(n1, n2, m1, psi); for (int i = ll; i <= x; ++i) { sumpi += pi[i - ll]; } } } if (!lower) sumpi = max(1 - sumpi, 0.0); if (give_log) { return sumpi == 0 ? JAGS_NEGINF : log(sumpi); } else { return sumpi; } } double DHyper::q(double p, vector const ¶meters, bool lower, bool log_p) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, parameters); int ll = max((int) 0, m1 - n2); int uu = min(n1, m1); vector pi = density(n1, n2, m1, psi); if (log_p) p = exp(p); if (!lower) p = 1 - p; double sumpi = 0; for (int i = ll; i < uu; ++i) { sumpi += pi[i - ll]; if (sumpi >= p) { return i; } } return uu; } double DHyper::r(vector const ¶meters, RNG *rng) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, parameters); int ll = max(0U, m1 - n2); int mode = modeCompute(n1, n2, m1, psi); vector pi = density(n1, n2, m1, psi); return sampleWithMode(mode - ll, pi, rng->uniform()) + ll; } double DHyper::l(vector const ¶meters) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, parameters); return max(0U, m1 - n2); } double DHyper::u(vector const ¶meters) const { int n1,n2,m1; double psi; getParameters(n1, n2, m1, psi, parameters); return min(n1, m1); } bool DHyper::isSupportFixed(vector const &fixmask) const { return fixmask[0] && fixmask[1] && fixmask[2]; //Margins fixed } double DHyper::KL(vector const ¶, vector const &parb) const { int n1a,n2a,m1a; double psia; getParameters(n1a, n2a, m1a, psia, para); int lla = max(0U, m1a - n2a); int uua = min(n1a, m1a); int n1b,n2b,m1b; double psib; getParameters(n1b, n2b, m1b, psib, para); int llb = max(0U, m1b - n2b); int uub = min(n1b, m1b); if (lla < llb || uua > uub) return JAGS_POSINF; vector da = density(n1a, n2a, m1a, psia); vector db = density(n1b, n2b, m1b, psib); double y = 0; for (int i = lla; i <= uua; ++i) { double proba = da[i - lla]; double probb = db[i - llb]; y += proba * (log(proba) - log(probb)); } return y; } }}