Index: trunk/npstat/stat/likelihoodStatisticCumulants.cc =================================================================== --- trunk/npstat/stat/likelihoodStatisticCumulants.cc (revision 598) +++ trunk/npstat/stat/likelihoodStatisticCumulants.cc (revision 599) @@ -1,805 +1,807 @@ #include #include #include #include "npstat/stat/likelihoodStatisticCumulants.hh" #define DEFAULT_CUM_LEN 6U #define MAX_ORDER 4U #define need_parameter(name) /**/ \ const double name = i_ ## name; \ assert(is_valid_shape_parameter(name)); static bool is_valid_shape_parameter(const double p) { return p < 0.9999999*std::numeric_limits::max(); } static unsigned nShapeParameters(const double gamma, const double rho, const double xi, const double zeta) { unsigned count = 0; if (is_valid_shape_parameter(gamma)) { ++count; if (is_valid_shape_parameter(rho)) { ++count; if (is_valid_shape_parameter(xi)) { ++count; if (is_valid_shape_parameter(zeta)) { ++count; } } } } return count; } static unsigned minNCumulants(const npstat::LikelihoodStatisticType t, const unsigned order) { if (t == npstat::L_STAT_R_LR) return order; else return order ? order + 2U : 0U; } static double mixtureCumulantLS( const unsigned cumN, const npstat::LikelihoodStatisticType t, const unsigned order, const unsigned N, const double i_gamma, const double i_rho, const double i_xi, const double i_zeta) { const double second = 1.0/N; const double first = sqrt(second); const double third = first*second; const double fourth = second*second; long double cum = 0.0L; switch (t) { case npstat::L_STAT_R_W: { switch (cumN) { case 1U: { if (order >= 3U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(2*gamma*rho - gamma3 - xi); } } break; case 2U: { cum += 1.0; if (order >= 2U) { need_parameter(gamma); need_parameter(rho); const double gamma2 = gamma*gamma; cum += second*(rho - 1 - gamma2); if (order >= 4U) { need_parameter(xi); need_parameter(zeta); const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(1 + 3*gamma2 - 10*gamma4 - 3*rho + 30*gamma2*rho - 8*rho2 - 22*gamma*xi + 10*zeta); } } } break; case 3U: { if (order) { need_parameter(gamma); cum += first*gamma; if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(21*gamma*rho - 3*gamma - 12*gamma3 - 9*xi); } } } break; case 4U: { if (order >= 2U) { need_parameter(rho); cum += second*(rho - 3); if (order >= 4U) { need_parameter(gamma); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(12 + 12*gamma2 - 114*gamma4 - 18*rho + 288*gamma2*rho - 60*rho2 - 180*gamma*xi + 66*zeta); } } } break; case 5U: { if (order >= 3U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(50*gamma*rho - 20*gamma - 30*gamma3 - 19*xi); } } break; case 6U: { if (order >= 4U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(60 - 40*gamma2 - 420*gamma4 - 45*rho + 930*gamma2*rho - 150*rho2 - 510*gamma*xi + 151*zeta); } } break; default: assert(0); } } break; case npstat::L_STAT_R_W2: { switch (cumN) { case 1U: { if (order) { need_parameter(gamma); cum += first*(-gamma/2); if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*((3*gamma)/16 + gamma3/4 + (5*gamma*rho)/16 - (9*xi)/8); } } } break; case 2U: { cum += 1.0; if (order >= 2U) { need_parameter(gamma); need_parameter(rho); const double gamma2 = gamma*gamma; cum += second*(3*rho - 5*gamma2/4); if (order >= 4U) { need_parameter(xi); need_parameter(zeta); const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*((19*gamma2)/16 + (41*gamma4)/4 - 3*rho - (443*gamma2*rho)/16 + 10*rho2 + (55*gamma*xi)/8 + 6*zeta); } } } break; case 3U: { if (order) { need_parameter(gamma); cum += first*(-2*gamma); if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*((3*gamma)/4 + (25*gamma3)/8 + (3*gamma*rho)/4 - 15*xi); } } } break; case 4U: { if (order >= 2U) { need_parameter(rho); cum += second*(10*rho); if (order >= 4U) { need_parameter(gamma); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(88*zeta - gamma2 + (243*gamma4)/8 - 10*rho - 132*gamma2*rho + 72*rho2 + 45*gamma*xi); } } } break; case 5U: { if (order >= 3U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(15*gamma3 - 40*gamma*rho - 54*xi); } } break; case 6U: { if (order >= 4U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(376*zeta - 10*gamma2 - 60*gamma4 - 150*gamma2*rho + 360*rho2 + 330*gamma*xi); } } break; default: assert(0); } } break; case npstat::L_STAT_R_W3: { switch (cumN) { case 1U: { if (order) { need_parameter(gamma); cum += first*(-gamma/2); if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*((3*gamma)/16 - (4*gamma3)/9 + (259*gamma*rho)/144 - (15*xi)/8); } } } break; case 2U: { cum += 1.0; if (order >= 2U) { need_parameter(gamma); need_parameter(rho); const double gamma2 = gamma*gamma; cum += second*((7*rho)/2 - 65*gamma2/36); if (order >= 4U) { need_parameter(xi); need_parameter(zeta); const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*((251*gamma2)/144 + (170*gamma4)/27 - (7*rho)/2 - (583*gamma2*rho)/48 + (131*rho2)/36 - (259*gamma*xi)/40 + (124*zeta)/9); } } } break; case 3U: { if (order) { need_parameter(gamma); cum += first*(-2*gamma); if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*((3*gamma)/4 - (5*gamma3)/24 + (107*gamma*rho)/12 - (39*xi)/2); } } } break; case 4U: { if (order >= 2U) { need_parameter(rho); cum += second*(10*rho); if (order >= 4U) { need_parameter(gamma); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(124*zeta - gamma2 + (163*gamma4)/8 - 10*rho - (748*gamma2*rho)/9 + (130*rho2)/3 - 4*gamma*xi); } } } break; case 5U: { if (order >= 3U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(15*gamma3 - 40*gamma*rho - 54*xi); } } break; case 6U: { if (order >= 4U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(376*zeta - 10*gamma2 - 60*gamma4 - 150*gamma2*rho + 360*rho2 + 330*gamma*xi); } } break; default: assert(0); } } break; case npstat::L_STAT_R_S: { switch (cumN) { case 1U: break; case 2U: cum += 1.0; break; case 3U: { if (order) { need_parameter(gamma); cum += first*gamma; } } break; case 4U: { if (order >= 2U) { need_parameter(rho); cum += second*(rho - 3); } } break; case 5U: { if (order >= 3U) { need_parameter(gamma); need_parameter(xi); cum += third*(xi - 10*gamma); } } break; case 6U: { if (order >= 4U) { need_parameter(gamma); need_parameter(rho); need_parameter(zeta); const double gamma2 = gamma*gamma; cum += fourth*(30 - 10*gamma2 - 15*rho + zeta); } } break; default: assert(0); } } break; case npstat::L_STAT_R_LR: { switch (cumN) { case 1U: { if (order) { need_parameter(gamma); cum -= first*gamma/6; if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(gamma/16 - gamma3/12 + (5*gamma*rho)/16 - (11*xi)/40); } } } break; case 2U: { cum += 1.0; if (order >= 2U) { need_parameter(gamma); need_parameter(rho); const double gamma2 = gamma*gamma; cum += second*(0.5*rho - 13.0/36*gamma2); if (order >= 4U) { need_parameter(xi); need_parameter(zeta); cum += fourth*(gamma2*(17.0/48 - gamma2/36 + 53*rho/48) - rho/2*(1 + rho) - 251*gamma*xi/120 + 5*zeta/3); } } } break; case 3U: { if (order >= 3U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*((-251*gamma3)/216 + (11*gamma*rho)/4 - (9*xi)/5); } } break; case 4U: { if (order >= 4U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(-1079*gamma4/216 + 33*gamma2*rho/2 - 35*rho2/6 - 66*gamma*xi/5 + 8*zeta); } } break; case 5U: break; case 6U: break; default: assert(0); } } break; default: assert(0); } return cum; } static double poissonCumulantLS( const unsigned cumN, const npstat::LikelihoodStatisticType t, const unsigned order, const double mu, const double i_gamma, const double i_rho, const double i_xi, const double i_zeta) { const double second = 1.0/mu; const double first = sqrt(second); const double third = first*second; const double fourth = second*second; long double cum = 0.0L; switch (t) { case npstat::L_STAT_R_S: { switch (cumN) { case 1U: break; case 2U: cum += 1.0; break; case 3U: { if (order) { need_parameter(gamma); cum += first*gamma; } } break; case 4U: { if (order >= 2U) { need_parameter(rho); cum += second*rho; } } break; case 5U: { if (order >= 3U) { need_parameter(xi); cum += third*xi; } } break; case 6U: { if (order >= 4U) { need_parameter(zeta); cum += fourth*zeta; } } break; default: assert(0); } } break; case npstat::L_STAT_R_LR: { switch (cumN) { case 1U: { if (order) { need_parameter(gamma); cum -= first*gamma/6; if (order >= 3U) { need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*(5*gamma*rho/16 - gamma3/12 - 11*xi/40); } } } break; case 2U: { cum += 1.0; if (order >= 2U) { need_parameter(gamma); need_parameter(rho); const double gamma2 = gamma*gamma; cum += second*(0.5*rho - 13.0/36*gamma2); if (order >= 4U) { need_parameter(xi); need_parameter(zeta); const double rho2 = rho*rho; cum += fourth*(gamma2*(53*rho/48 - gamma2/36) - rho2/2 - 251*gamma*xi/120 + 5*zeta/3); } } } break; case 3U: { if (order >= 3U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); const double gamma3 = gamma*gamma*gamma; cum += third*((-251*gamma3)/216 + (11*gamma*rho)/4 - (9*xi)/5); } } break; case 4U: { if (order >= 4U) { need_parameter(gamma); need_parameter(rho); need_parameter(xi); need_parameter(zeta); const double gamma2 = gamma*gamma; const double gamma4 = gamma2*gamma2; const double rho2 = rho*rho; cum += fourth*(-1079*gamma4/216 + 33*gamma2*rho/2 - 35*rho2/6 - 66*gamma*xi/5 + 8*zeta); } } break; case 5U: break; case 6U: break; default: assert(0); } } break; default: assert(0); } return cum; } namespace npstat { std::vector mixtureModelCumulants( const LikelihoodStatisticType t, const unsigned order, const unsigned N, const double gamma, const double rho, const double xi, const double zeta) { std::vector cumulants(DEFAULT_CUM_LEN, 0.0); + cumulants[1] = 1.0; if (!N) throw std::invalid_argument("In npstat::mixtureModelCumulants: " "sample size must be positive"); // Check that we know how to build the cumulants up to the order requested if (order > MAX_ORDER) throw std::invalid_argument("In npstat::mixtureModelCumulants: " "requested approximation order is " "too large"); // Check that enough shape parameters were provided if (nShapeParameters(gamma, rho, xi, zeta) < order) throw std::invalid_argument("In npstat::mixtureModelCumulants: " "not enough shape parameters for " "the requested approximation order"); const unsigned cumulantsNeeded = minNCumulants(t, order); assert(cumulants.size() >= cumulantsNeeded); for (unsigned i=0; i poissonProcessCumulants( const LikelihoodStatisticType t, const unsigned order, const double mu, const double gamma, const double rho, const double xi, const double zeta) { std::vector cumulants(DEFAULT_CUM_LEN, 0.0); + cumulants[1] = 1.0; // Check that we know how to build the cumulants for this statistic if (!(t == L_STAT_R_S || t == L_STAT_R_LR)) throw std::invalid_argument("In npstat::poissonProcessCumulants: " "requested statistic is not supported"); if (mu <= 0.0) throw std::invalid_argument("In npstat::poissonProcessCumulants: " "background rate must be positive"); // Check that we know how to build the cumulants up to the order requested if (order > MAX_ORDER) throw std::invalid_argument("In npstat::poissonProcessCumulants: " "requested approximation order is " "too large"); // Check that enough shape parameters were provided if (nShapeParameters(gamma, rho, xi, zeta) < order) throw std::invalid_argument("In npstat::poissonProcessCumulants: " "not enough shape parameters for " "the requested approximation order"); const unsigned cumulantsNeeded = minNCumulants(t, order); assert(cumulants.size() >= cumulantsNeeded); for (unsigned i=0; i