Index: trunk/tests/test_likelihoodStatisticCumulants.cc
===================================================================
--- trunk/tests/test_likelihoodStatisticCumulants.cc	(revision 602)
+++ trunk/tests/test_likelihoodStatisticCumulants.cc	(revision 603)
@@ -1,109 +1,107 @@
 #include "test_utils.hh"
 #include "UnitTest++.h"
 
 #include "npstat/stat/likelihoodStatisticCumulants.hh"
 
 using namespace npstat;
 using namespace std;
 
 #define N_STATISTICS 5
 #define N_CUMULANTS 6
 
 namespace {
     const double gamma = 3.0;
     const double rho = 5.0;
     const double xi = 7.0;
     const double zeta = 11.0;
     const unsigned N = 13;
 
     static void slr_cumulants_of_poisson_counts(const double mu, double cumulants[6])
     {
         cumulants[0] = -1.0/6.0*pow(mu,-0.5) - 11.0/240.0*pow(mu,-1.5);
         cumulants[1] = 1.0 + 5/36.0/mu + 109/720.0/mu/mu;
         cumulants[2] = -229/1080.0*pow(mu,-1.5);
         cumulants[3] = 509/1080.0/mu/mu;
         cumulants[4] = 0.0;
         cumulants[5] = 0.0;
     }
 
     TEST(likelihoodStatisticCumulants_mix_4)
     {
         const double eps = 1.0e-12;
         const unsigned order = 4;
 
         const double expectations[N_STATISTICS][N_CUMULANTS] = {
             {-0.08533849172695833, 0.6213017751479290, -0.8960541631330624, -4.568047337278107, -5.397659601730114, -32.50887573964497},
             {-0.3280198275754961, 1.527366863905325, -1.816109777064332, 4.877958579881657, -12.22473893988678, 49.50295857988166},
             {-0.3644664750838845, 1.499572649572650, -1.794775154132592, 5.116617357001972, -12.22473893988678, 49.50295857988166},
             {0.0, 1.0, 0.8320502943378437, 0.1538461538461538, -0.4906963274300104, -0.7337278106508876},
             {-0.1237408130040896, 1.001676528599606, -0.05813684748899036, 0.01681459566074951, 0.0, 0.0}
         };
 
         for (unsigned stat=0; stat<N_STATISTICS; ++stat)
         {
             const vector<double>& cumulants = mixtureModelCumulants(
                 static_cast<LikelihoodStatisticType>(stat),
                 order, N, gamma, rho, xi, zeta);
             CHECK(cumulants.size() == N_CUMULANTS);
             for (unsigned i=0; i<N_CUMULANTS; ++i)
                 CHECK_CLOSE(expectations[stat][i], cumulants[i], eps);
         }
     }
 
     TEST(likelihoodStatisticCumulants_poiss_4)
     {
         const double eps = 1.0e-12;
         const unsigned order = 4;
 
-        const double w_expectations[N_CUMULANTS] = {-0.08533849172695833, 0.6213017751479290, -0.7040425567474062,
-                                                    -4.514792899408284, -4.117582225825739, -29.40236686390533};
-        const double score_expectations[N_CUMULANTS] = {0., 1., 0.8320502943378437, 0.3846153846153846,
-                                                        0.1493423605221771, 0.06508875739644970};
-        const double slr_expectations[N_CUMULANTS] = {-0.1277410548037907, 0.9976084812623274,
-                                                      -0.05813684748899036, 0.01681459566074951, 0., 0.};
-
-        vector<double> cumulants = poissonProcessCumulants(
-            L_STAT_R_W, order, N, gamma, rho, xi, zeta);
-        CHECK(cumulants.size() == N_CUMULANTS);
-        for (unsigned i=0; i<N_CUMULANTS; ++i)
-            CHECK_CLOSE(w_expectations[i], cumulants[i], eps);
-
-        cumulants = poissonProcessCumulants(
-            L_STAT_R_S, order, N, gamma, rho, xi, zeta);
-        CHECK(cumulants.size() == N_CUMULANTS);
-        for (unsigned i=0; i<N_CUMULANTS; ++i)
-            CHECK_CLOSE(score_expectations[i], cumulants[i], eps);
-
-        cumulants = poissonProcessCumulants(
-            L_STAT_R_LR, order, N, gamma, rho, xi, zeta);
-        CHECK(cumulants.size() == N_CUMULANTS);
-        for (unsigned i=0; i<N_CUMULANTS; ++i)
-            CHECK_CLOSE(slr_expectations[i], cumulants[i], eps);
+        const double expectations[N_STATISTICS][N_CUMULANTS] = {
+            {-0.08533849172695833, 0.6213017751479290, -0.7040425567474062,
+             -4.514792899408284, -4.117582225825739, -29.40236686390533},
+            {-0.3400205529745996, 1.552884615384615, -1.864112678660746,
+             5.227071005917160, -12.22473893988678, 50.03550295857988},
+            {-0.3764672004829880, 1.510297501643655, -1.842778055729006,
+             5.465729783037475, -12.22473893988678, 50.03550295857988},
+            {0., 1., 0.8320502943378437, 0.3846153846153846,
+             0.1493423605221771, 0.06508875739644970},
+            {-0.1277410548037907, 0.9976084812623274,
+             -0.05813684748899036, 0.01681459566074951, 0., 0.}
+        };
+
+        for (unsigned stat=0; stat<N_STATISTICS; ++stat)
+        {
+            const vector<double>& cumulants = poissonProcessCumulants(
+                static_cast<LikelihoodStatisticType>(stat),
+                order, N, gamma, rho, xi, zeta);
+            CHECK(cumulants.size() == N_CUMULANTS);
+            for (unsigned i=0; i<N_CUMULANTS; ++i)
+                CHECK_CLOSE(expectations[stat][i], cumulants[i], eps);
+        }
     }
 
     TEST(likelihoodStatisticCumulants_poiss_counts)
     {
         const double eps = 1.0e-12;
         const double muvalues[] = {7.0, 13.0, 29.0};
         const unsigned nMu = sizeof(muvalues)/sizeof(muvalues[0]);
 
         for (unsigned imu=0; imu<nMu; ++imu)
         {
             const double mu = muvalues[imu];
 
             const vector<double>& cumulants = poissonProcessCumulants(
                 L_STAT_R_LR, 4U, mu, 1.0, 1.0, 1.0, 1.0);
             CHECK(cumulants.size() >= 6);
 
             double expectations[6];
             slr_cumulants_of_poisson_counts(mu, expectations);
 
             CHECK_CLOSE(expectations[0], cumulants[0], eps);
             CHECK_CLOSE(expectations[1], cumulants[1], eps);
             CHECK_CLOSE(expectations[2], cumulants[2], eps);
             CHECK_CLOSE(expectations[3], cumulants[3], eps);
             CHECK_CLOSE(expectations[4], cumulants[4], eps);
             CHECK_CLOSE(expectations[5], cumulants[5], eps);
         }
     }
 }
Index: trunk/npstat/stat/likelihoodStatisticCumulants.cc
===================================================================
--- trunk/npstat/stat/likelihoodStatisticCumulants.cc	(revision 602)
+++ trunk/npstat/stat/likelihoodStatisticCumulants.cc	(revision 603)
@@ -1,921 +1,1154 @@
 #include <cmath>
 #include <cassert>
 #include <stdexcept>
 
 #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<double>::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_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*(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*(41*gamma4/4 - 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*(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*(243*gamma4/8 - 132*gamma2*rho + 72*rho2 + 45*gamma*xi + 88*zeta);
+                }
+            }
+        }
+        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 - 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*((-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*(170*gamma4/27 - 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*(107*gamma*rho/12 - 5*gamma3/24 - 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*(163*gamma4/8 - 748*gamma2*rho/9 + 130*rho2/3 - 4*gamma*xi + 124*zeta);
+                }
+            }
+        }
+        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 - 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;
             }
         }
         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;
 
     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 - gamma2);
 
                 if (order >= 4U)
                 {
                     need_parameter(xi);
                     need_parameter(zeta);
                     const double gamma4 = gamma2*gamma2;
                     const double rho2 = rho*rho;
                     cum += fourth*(10*zeta - 10*gamma4 + 30*gamma2*rho - 8*rho2 - 22*gamma*xi);
                 }
             }
         }
         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 - 12*gamma3 - 9*xi);
                 }
             }
         }
         break;
 
         case 4U:
         {
             if (order >= 2U)
             {
                 need_parameter(rho);
                 cum += second*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*(66*zeta - 114*gamma4 + 288*gamma2*rho - 60*rho2 - 180*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*(50*gamma*rho - 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*(151*zeta - 420*gamma4 + 930*gamma2*rho - 150*rho2 - 510*gamma*xi);
             }
         }
         break;
 
         default:
             assert(0);
         }
     }
     break;
 
     default:
         assert(0);
     }
 
     return cum;
 }
 
 namespace npstat {
     std::vector<double> mixtureModelCumulants(
         const LikelihoodStatisticType t,
         const unsigned order, const unsigned N,
         const double gamma, const double rho,
         const double xi, const double zeta)
     {
         std::vector<double> 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<cumulantsNeeded; ++i)
             cumulants[i] = mixtureCumulantLS(i+1, t, order, N, gamma, rho, xi, zeta);
 
         return cumulants;
     }
 
     std::vector<double> poissonProcessCumulants(
         const LikelihoodStatisticType t,
         const unsigned order, const double mu,
         const double gamma, const double rho,
         const double xi, const double zeta)
     {
         std::vector<double> 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 || t == L_STAT_R_W))
-            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<cumulantsNeeded; ++i)
             cumulants[i] = poissonCumulantLS(i+1, t, order, mu, gamma, rho, xi, zeta);
 
         return cumulants;
     }
 }