diff --git a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp index 247f415af9f..a464ae3780d 100644 --- a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp @@ -3,17 +3,16 @@ #include #include +#include #include -#include -#include +#include #include +#include #include #include -#include #include #include #include -#include #include namespace stan { @@ -52,13 +51,20 @@ return_type_t neg_binomial_2_log_lpmf( size_t size_phi = stan::math::size(phi); size_t size_eta_phi = max_size(eta, phi); size_t size_n_phi = max_size(n, phi); - size_t max_size_seq_view = max_size(n, eta, phi); + size_t size_all = max_size(n, eta, phi); VectorBuilder eta_val(size_eta); for (size_t i = 0; i < size_eta; ++i) { eta_val[i] = value_of(eta_vec[i]); } + VectorBuilder phi_val(size_phi); + VectorBuilder log_phi(size_phi); + for (size_t i = 0; i < size_phi; ++i) { + phi_val[i] = value_of(phi_vec[i]); + log_phi[i] = log(phi_val[i]); + } + VectorBuilder::value, T_partials_return, T_log_location> exp_eta(size_eta); @@ -68,17 +74,19 @@ return_type_t neg_binomial_2_log_lpmf( } } - VectorBuilder phi_val(size_phi); - VectorBuilder log_phi(size_phi); - for (size_t i = 0; i < size_phi; ++i) { - phi_val[i] = value_of(phi_vec[i]); - log_phi[i] = log(phi_val[i]); + VectorBuilder::value, + T_partials_return, T_log_location, T_precision> + exp_eta_over_exp_eta_phi(size_eta_phi); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_eta_phi; ++i) { + exp_eta_over_exp_eta_phi[i] = inv(phi_val[i] / exp_eta[i] + 1); + } } VectorBuilder - logsumexp_eta_logphi(size_eta_phi); + log1p_exp_eta_m_logphi(size_eta_phi); for (size_t i = 0; i < size_eta_phi; ++i) { - logsumexp_eta_logphi[i] = log_sum_exp(eta_val[i], log_phi[i]); + log1p_exp_eta_m_logphi[i] = log1p_exp(eta_val[i] - log_phi[i]); } VectorBuilder n_plus_phi( @@ -87,38 +95,25 @@ return_type_t neg_binomial_2_log_lpmf( n_plus_phi[i] = n_vec[i] + phi_val[i]; } - for (size_t i = 0; i < max_size_seq_view; i++) { - if (phi_val[i] > 1e5) { - // TODO(martinmodrak) This is wrong (doesn't pass propto information), - // and inaccurate for n = 0, but shouldn't break most models. - // Also the 1e5 cutoff is way too low. - // Will be addressed better once PR #1497 is merged - logp += poisson_log_lpmf(n_vec[i], eta_val[i]); - } else { - if (include_summand::value) { - logp -= lgamma(n_vec[i] + 1.0); - } - if (include_summand::value) { - logp += multiply_log(phi_val[i], phi_val[i]) - lgamma(phi_val[i]); - } - if (include_summand::value) { - logp += n_vec[i] * eta_val[i]; - } - if (include_summand::value) { - logp += lgamma(n_plus_phi[i]); - } - logp -= (n_plus_phi[i]) * logsumexp_eta_logphi[i]; + for (size_t i = 0; i < size_all; i++) { + if (include_summand::value) { + logp += binomial_coefficient_log(n_plus_phi[i] - 1, n_vec[i]); + } + if (include_summand::value) { + logp += n_vec[i] * eta_val[i]; } + logp += -phi_val[i] * log1p_exp_eta_m_logphi[i] + - n_vec[i] * (log_phi[i] + log1p_exp_eta_m_logphi[i]); if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] - += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp_eta[i] + 1); + += n_vec[i] - n_plus_phi[i] * exp_eta_over_exp_eta_phi[i]; } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += 1.0 - n_plus_phi[i] / (exp_eta[i] + phi_val[i]) + log_phi[i] - - logsumexp_eta_logphi[i] - digamma(phi_val[i]) - + digamma(n_plus_phi[i]); + += exp_eta_over_exp_eta_phi[i] - n_vec[i] / (exp_eta[i] + phi_val[i]) + - log1p_exp_eta_m_logphi[i] + - (digamma(phi_val[i]) - digamma(n_plus_phi[i])); } } return ops_partials.build(logp); diff --git a/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp index a9887050958..7742732d86f 100644 --- a/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp +++ b/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -212,29 +213,21 @@ TEST(ProbNegBinomial2, log_matches_lpmf) { TEST(ProbDistributionsNegBinomial2Log, neg_binomial_2_log_grid_test) { std::vector mu_log_to_test = {-101, -27, -3, -1, -0.132, 0, 4, 10, 87}; - // TODO(martinmodrak) Reducing the span of the test, should be fixed - // along with #1495 - // std::vector phi_to_test = {2e-5, 0.36, 1, 10, 2.3e5, 1.8e10, 6e16}; - std::vector phi_to_test = {0.36, 1, 10}; + std::vector phi_to_test = {2e-5, 0.36, 1, 10, 2.3e5, 1.8e10, 6e16}; std::vector n_to_test = {0, 1, 10, 39, 101, 3048, 150054}; - // TODO(martinmdorak) Only weak tolerance for this quick fix - auto tolerance = [](double x) { return std::max(fabs(x * 1e-8), 1e-8); }; - for (double mu_log : mu_log_to_test) { for (double phi : phi_to_test) { for (int n : n_to_test) { double val_log = stan::math::neg_binomial_2_log_lpmf(n, mu_log, phi); - EXPECT_LE(val_log, 0) - << "neg_binomial_2_log_lpmf yields " << val_log - << " which si greater than 0 for n = " << n - << ", mu_log = " << mu_log << ", phi = " << phi << "."; + std::stringstream msg; double val_orig = stan::math::neg_binomial_2_lpmf(n, std::exp(mu_log), phi); - EXPECT_NEAR(val_log, val_orig, tolerance(val_orig)) + msg << std::setprecision(22) << "neg_binomial_2_log_lpmf yields different result (" << val_log << ") than neg_binomial_2_lpmf (" << val_orig << ") for n = " << n << ", mu_log = " << mu_log << ", phi = " << phi << "."; + stan::test::expect_near_rel(msg.str(), val_log, val_orig); } } } diff --git a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp index fede0d0850f..21da83beab2 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -26,394 +27,448 @@ struct TestValue { }; // Test data generated in Mathematica (Wolfram Cloud). The code can be re-ran at -// https://www.wolframcloud.com/obj/martin.modrak/Published/NegBinomial_2_Log_test.nb +// https://www.wolframcloud.com/obj/martin.modrak/Published/neg_binomial_2_log_lpmf.nb // but is also presented below for convenience: // -// nb2log[n_,eta_,phi_]:= LogGamma[n + phi] - LogGamma[n + 1] - -// LogGamma[phi ] + n * (eta - Log[Exp[eta] + phi]) + -// phi * (Log[phi] - Log[Exp[eta] + phi]); -// nb2logdeta[n_,eta_,phi_]= D[nb2log[n, eta, phi],eta]; -// nb2logdphi[n_,eta_,phi_]= D[nb2log[n, eta, phi],phi]; -// out = OpenWrite["nb2log_test.txt"] -// etas= {-130,-50,-1, -135*10^-3,0,98 *10^-2,5,368}; -// phis= {4*10^-4,65*10^-3,442*10^-2,800, 15324}; -// ns = {0,6,14,1525,10233}; -// WriteString[out, "std::vector testValues = {"]; -// Block[{$MaxPrecision = 80, $MinPrecision = 40}, { -// For[i = 1, i <= Length[etas], i++, { -// For[j = 1, j <= Length[phis], j++, { -// For[k = 1, k <= Length[ns], k++, { -// ceta = etas[[i]]; -// cphi = phis[[j]]; -// cn=ns[[k]]; -// val = N[nb2log[cn,ceta,cphi]]; -// ddeta= N[nb2logdeta[cn,ceta,cphi]]; -// ddphi= N[nb2logdphi[cn,ceta,cphi]]; -// WriteString[out," TestValue(",CForm[cn],",",CForm[ceta],",", -// CForm[cphi],",", CForm[val],","CForm[ddeta],",", -// CForm[ddphi],"),"] -// }] -// }] +// toCString[x_] := ToString[CForm[N[x, 24]]]; +// nb2log[n_,logmu_,phi_]:= LogGamma[n + phi] - LogGamma[n + 1] - +// LogGamma[phi] + n * (logmu - Log[Exp[logmu] + phi]) + +// phi * (Log[phi] - Log[Exp[logmu] + phi]) +// nb2logdmu[n_,logmu_,phi_]= D[nb2log[n, logmu, phi],logmu]; +// nb2logdphi[n_,logmu_,phi_]= D[nb2log[n, logmu, phi],phi]; +// logmus= SetPrecision[{-36,-26.1, -4.5,-0.2,2.0,13.4, 84.3}, Infinity]; +// phis= SetPrecision[{0.0001,0.084,3.76,461.0, 142311.0}, Infinity]; +// ns = {0,7,11,2338,9611}; +// out = "std::vector testValues = {\n"; +// For[k = 1, k <= Length[ns], k++, { +// For[i = 1, i <= Length[logmus], i++, { +// For[j = 1, j <= Length[phis], j++, { +// clogmu = logmus[[i]]; +// cphi = phis[[j]]; +// cn=ns[[k]]; +// val = nb2log[cn,clogmu,cphi]; +// ddmu= nb2logdmu[cn,clogmu,cphi]; +// ddphi= nb2logdphi[cn,clogmu,cphi]; +// out = StringJoin[out," {",ToString[cn],",",toCString[clogmu],",", +// toCString[cphi],",", toCString[val],",",toCString[ddmu],",", +// toCString[ddphi],"},\n"]; +// }] // }] -// }]; -// WriteString[out,"};"]; -// Close[out]; -// FilePrint[%] +// }] +// out = StringJoin[out,"};\n"]; +// out std::vector testValues = { - TestValue(0, -130, 0.0004, 0., -3.4811068399043105e-57, 0.), - TestValue(6, -130, 0.0004, -742.670616198677, 6., -12497.717251921473), - TestValue(14, -130, 0.0004, -1720.9251872606271, 14., -32496.820494410564), - TestValue(1525, -130, 0.0004, -186333.4804666206, 1525., - -3.80999209402007e6), - TestValue(10233, -130, 0.0004, -1.250243590665905e6, 10233., - -2.557999019011787e7), - TestValue(0, -130, 0.065, 0., -3.4811068399043105e-57, 0.), - TestValue(6, -130, 0.065, -767.9794906909146, 6., -74.73014818539053), - TestValue(14, -130, 0.065, -1786.9017778439902, 14., -196.917187021929), - TestValue(1525, -130, 0.065, -194091.16634050777, 1525., - -23438.34928792469), - TestValue(10233, -130, 0.065, -1.3023307775882862e6, 10233., - -157405.6761910148), - TestValue(0, -130, 4.42, 0., -3.4811068399043105e-57, 0.), - TestValue(6, -130, 4.42, -784.0828117852828, 6., -0.4312637869430638), - TestValue(14, -130, 4.42, -1833.6283513006856, 14., -1.6501470360217263), - TestValue(1525, -130, 4.42, -200493.63376876694, 1525., -339.0590808438865), - TestValue(10233, -130, 4.42, -1.3454684320637877e6, 10233., - -2307.29338827833), - TestValue(0, -130, 800, 0., -3.4811068399043105e-57, 0.), - TestValue(6, -130, 800, -786.5605440348701, 6., -0.000023330624470518056), - TestValue(14, -130, 800, -1845.078105689095, 14., -0.00014060783924468162), - TestValue(1525, -130, 800, -206952.5277963584, 1525., -0.8389763486176387), - TestValue(10233, -130, 800, -1.395830983586415e6, 10233., - -10.166635855973585), - TestValue(0, -130, 15324, 0., -3.4811068399043105e-57, 0.), - TestValue(6, -130, 15324, -786.5782724723963, 6., -6.386208176026916e-8), - TestValue(14, -130, 15324, -1845.1852845285289, 14., -3.872952154126844e-7), - TestValue(1525, -130, 15324, -207834.0167357739, 1525., - -0.004643063079575183), - TestValue(10233, -130, 15324, -1.4117087218942088e6, 10233., - -0.15627194802547817), - TestValue(0, -50, 0.0004, 0., -1.9287498479639178e-22, 0.), - TestValue(6, -50, 0.0004, -262.6706161986769, 6., -12497.717251921473), - TestValue(14, -50, 0.0004, -600.9251872606274, 14., -32496.820494410564), - TestValue(1525, -50, 0.0004, -64333.4804666206, 1525., -3.80999209402007e6), - TestValue(10233, -50, 0.0004, -431603.59066590504, 10233., - -2.557999019011787e7), - TestValue(0, -50, 0.065, 0., -1.9287498479639178e-22, 0.), - TestValue(6, -50, 0.065, -287.9794906909146, 6., -74.73014818539053), - TestValue(14, -50, 0.065, -666.9017778439904, 14., -196.917187021929), - TestValue(1525, -50, 0.065, -72091.16634050779, 1525., -23438.34928792469), - TestValue(10233, -50, 0.065, -483690.77758828626, 10233., - -157405.6761910148), - TestValue(0, -50, 4.42, 0., -1.9287498479639178e-22, 0.), - TestValue(6, -50, 4.42, -304.08281178528296, 6., -0.4312637869430638), - TestValue(14, -50, 4.42, -713.6283513006858, 14., -1.6501470360217263), - TestValue(1525, -50, 4.42, -78493.63376876694, 1525., -339.0590808438865), - TestValue(10233, -50, 4.42, -526828.4320637881, 10233., -2307.29338827833), - TestValue(0, -50, 800, 0., -1.9287498479639178e-22, 0.), - TestValue(6, -50, 800, -306.5605440348702, 6., -0.000023330624470518056), - TestValue(14, -50, 800, -725.0781056890954, 14., -0.00014060783924468162), - TestValue(1525, -50, 800, -84952.52779635842, 1525., -0.8389763486176387), - TestValue(10233, -50, 800, -577190.983586415, 10233., -10.166635855973585), - TestValue(0, -50, 15324, 0., -1.9287498479639178e-22, 0.), - TestValue(6, -50, 15324, -306.5782724723962, 6., -6.386208176026916e-8), - TestValue(14, -50, 15324, -725.1852845285287, 14., -3.872952154126844e-7), - TestValue(1525, -50, 15324, -85834.01673577388, 1525., - -0.004643063079575183), - TestValue(10233, -50, 15324, -593068.7218942087, 10233., - -0.15627194802547817), - TestValue(0, -1, 0.0004, -0.002730053093156532, -0.00039956554729340564, - -5.8262188646578155), - TestValue(6, -1, 0.0004, -9.624142649118, 0.00611722505162245, - 2480.164552716579), - TestValue(14, -1, 0.0004, -10.479775574199106, 0.014806279183510258, - 2459.33867489777), - TestValue(1525, -1, 0.0004, -16.81061433296054, 1.65595137834382, - -1638.7975986621), - TestValue(10233, -1, 0.0004, -28.176651635192684, 11.113986800903698, - -25281.98225286028), - TestValue(0, -1, 0.065, -0.12324467951140204, -0.05523977671804774, - -1.0462292737439123), - TestValue(6, -1, 0.065, -5.4791673253246715, 0.8457039108467774, - 2.6706427321759207), - TestValue(14, -1, 0.065, -7.57003041826529, 2.046962160933211, - -14.9203691826154), - TestValue(1525, -1, 0.065, -257.79937372403765, 228.93461414600833, - -3500.7778852403826), - TestValue(10233, -1, 0.065, -1676.4055320453917, 1536.504219365091, - -23615.329484008667), - TestValue(0, -1, 4.42, -0.35336977067787206, -0.33961321498519975, - -0.0031123429168939953), - TestValue(6, -1, 4.42, -10.915869027469208, 5.199373371007742, - -0.3300744524827952), - TestValue(14, -1, 4.42, -29.10099183821663, 12.58468881899833, - -1.409888798391908), - TestValue(1525, -1, 4.42, -3890.907704212679, 1407.4861440582208, - -312.5521835201079), - TestValue(10233, -1, 4.42, -26229.892416216317, 9446.402009195977, - -2129.4099898544955), - TestValue(0, -1, 800, -0.36779488254126136, -0.3677103498239346, - -1.0566589647709179e-7), - TestValue(6, -1, 800, -12.931097379030476, 5.629531822552386, - -0.000019989005837535087), - TestValue(14, -1, 800, -39.45233698208112, 13.625854719054148, - -0.00013266984123889644), - TestValue(1525, -1, 800, -10228.596700235812, 1523.9313417958242, - -0.8381002694655955), - TestValue(10233, -1, 800, -75779.05593758884, 10227.928814637991, - -10.160756617874249), - TestValue(0, -1, 15324, -0.36787502543801764, -0.36787060979358055, - -2.8815172470331163e-10), - TestValue(6, -1, 15324, -12.946291536606935, 5.631985353162647, - -5.475079230526231e-8), - TestValue(14, -1, 15324, -39.55349564443637, 13.631793303770952, - -3.656513370973221e-7), - TestValue(1525, -1, 15324, -11109.421220654047, 1524.5955199749144, - -0.004640674342975881), - TestValue(10233, -1, 15324, -91652.33542736099, 10232.386474212051, - -0.1562559175659044), - TestValue(0, -0.135, 0.0004, -0.0030758014883219406, -0.0003998169579138641, - -6.689961326020191), - TestValue(6, -0.135, 0.0004, -9.620714324994289, 0.0023458143341247776, - 2488.7287085224075), - TestValue(14, -0.135, 0.0004, -10.471315153383557, 0.006006656056842966, - 2480.473361726521), - TestValue(1525, -0.135, 0.0004, -15.851716649474838, 0.6974481364352408, - 756.5961351210399), - TestValue(10233, -0.135, 0.0004, -21.74031670264958, 4.68227435161399, - -9203.56550062657), - TestValue(0, -0.135, 0.065, -0.17355816530480334, -0.060499170784905056, - -1.7393691464599734), - TestValue(6, -0.135, 0.065, -4.983802576662806, 0.3549619875315516, - 9.446464847896287), - TestValue(14, -0.135, 0.065, -6.347094690329737, 0.9089101986201605, - 1.814068917686857), - TestValue(1525, -0.135, 0.065, -119.1564692857828, 105.53587856798117, - -1603.1098530521738), - TestValue(10233, -0.135, 0.065, -745.8016166677262, 708.5085063379321, - -10877.707952603063), - TestValue(0, -0.135, 4.42, -0.7972828882914406, -0.7295110644556051, - -0.015332991818062247), - TestValue(6, -0.135, 4.42, -6.77237913731839, 4.280202422708732, - -0.22255005640011594), - TestValue(14, -0.135, 4.42, -18.840964604379963, 10.959820405594515, - -1.142704342330766), - TestValue(1525, -0.135, 4.42, -2725.386686190179, 1272.5726669231467, - -282.1292052356146), - TestValue(10233, -0.135, 4.42, -18406.520499591774, 8543.336841294322, - -1925.1970362834456), - TestValue(0, -0.135, 800, -0.8732391466033107, -0.8727627283344418, - -5.955228363774268e-7), - TestValue(6, -0.135, 800, -8.250332475072994, 5.12069155120305, - -0.00001574399672842758), - TestValue(14, -0.135, 800, -27.85662652076425, 13.111963923919705, - -0.00012211167739840934), - TestValue(1525, -0.135, 800, -8910.940647628244, 1522.463533320778, - -0.8368973142018654), - TestValue(10233, -0.135, 800, -66934.48164579559, 10220.963511022857, - -10.15268179368541), - TestValue(0, -0.135, 15324, -0.8736910046546811, -0.873666098586244, - -1.625299006491332e-9), - TestValue(6, -0.135, 15324, -8.262305564366748, 5.125991823849705, - -4.3164385630243066e-8), - TestValue(14, -0.135, 15324, -27.949773736920378, 13.125535720430971, - -3.3683352693505975e-7), - TestValue(1525, -0.135, 15324, -9790.852373971315, 1524.0393891872175, - -0.004637390943832059), - TestValue(10233, -0.135, 15324, -82801.63401513056, 10231.542920615924, - -0.15623387778406472), - TestValue(0, 0, 0.0004, -0.0031297783723510477, -0.00039984006397441024, - -6.824845770941593), - TestValue(6, 0, 0.0004, -9.620421562314927, 0.001999200319872051, - 2489.460301347968), - TestValue(14, 0, 0.0004, -10.470560071286346, 0.005197920831667333, - 2482.3602575793902), - TestValue(1525, 0, 0.0004, -15.763640987330291, 0.6093562574970012, - 976.6908902568914), - TestValue(10233, 0, 0.0004, -21.149006354171433, 4.091163534586165, - -7725.923400264901), - TestValue(0, 0, 0.065, -0.18176228253611273, -0.06103286384976526, - -1.8573756720976538), - TestValue(6, 0, 0.065, -4.939309822938063, 0.3051643192488263, - 10.086365633302695), - TestValue(14, 0, 0.065, -6.2323394419969205, 0.7934272300469484, - 3.464512784485386), - TestValue(1525, 0, 0.065, -105.77088536835072, 93.01408450704226, - -1410.5930846874387), - TestValue(10233, 0, 0.065, -655.9353073694074, 624.4882629107981, - -9585.215040143006), - TestValue(0, 0, 4.42, -0.9014860475798108, -0.8154981549815498, - -0.019454274343498), - TestValue(6, 0, 4.42, -6.20803454903441, 4.077490774907749, - -0.2002630680488462), - TestValue(14, 0, 4.42, -17.385223019332848, 10.601476014760147, - -1.0852063261438887), - TestValue(1525, 0, 4.42, -2555.5683368414957, 1242.819188191882, - -275.42122433697745), - TestValue(10233, 0, 4.42, -17266.416519266422, 8344.177121771218, - -1880.16185158575), - TestValue(0, 0, 800, -0.9993755203460353, -0.9987515605493134, - -7.799497456772997e-7), - TestValue(6, 0, 800, -7.567414871619192, 4.9937578027465666, - -0.00001474727833628009), - TestValue(14, 0, 800, -26.094970281047694, 12.983770287141073, - -0.00011954009860382087), - TestValue(1525, 0, 800, -8705.432231464434, 1522.0973782771537, - -0.8365972908645132), - TestValue(10233, 0, 800, -65554.76622405997, 10219.225967540575, - -10.150667534799734), - TestValue(0, 0, 15324, -0.9999673728476992, -0.9999347471451876, - -2.129059595290528e-9), - TestValue(6, 0, 15324, -7.578631375147566, 4.999673735725938, - -4.0441863191631455e-8), - TestValue(14, 0, 15324, -26.18616547115994, 12.99915171288744, - -3.298092945414055e-7), - TestValue(1525, 0, 15324, -9585.116216997136, 1523.9005546492658, - -0.004636571433929149), - TestValue(10233, 0, 15324, -81420.38961583155, 10231.33233278956, - -0.15622837586171023), - TestValue(0, 0.98, 0.0004, -0.003521678449611314, -0.0003999399592377696, - -7.8043462259338625), - TestValue(6, 0.98, 0.0004, -9.619314621296192, 0.0005006714742186634, - 2492.2268732689504), - TestValue(14, 0.98, 0.0004, -10.46745467547295, 0.0017014867188272407, - 2490.1215926683385), - TestValue(1525, 0.98, 0.0004, -15.383077442174908, 0.2285054660442723, - 1927.8381186954143), - TestValue(10233, 0.98, 0.0004, -18.59312476501509, 1.5355928598007087, - -1337.9764634954095), - TestValue(0, 0.98, 0.065, -0.2429355772781735, -0.06345207263809681, - -2.76128468677041), - TestValue(6, 0.98, 0.065, -4.767248786178039, 0.0794335299991208, - 12.618019394958765), - TestValue(14, 0.98, 0.065, -5.74929929656745, 0.26994766684874427, - 10.576916914579899), - TestValue(1525, 0.98, 0.065, -46.55166607298088, 36.25330526432138, - -538.291454718534), - TestValue(10233, 0.98, 0.065, -258.21532828733325, 243.62794322513656, - -3726.7666341288596), - TestValue(0, 0.98, 4.42, -2.085194344338019, -1.6623571643545987, - -0.09566452035824002), - TestValue(6, 0.98, 4.42, -3.118586687545845, 2.0810493998699275, - -0.01638680599463549), - TestValue(14, 0.98, 4.42, -8.598233613515248, 7.072258152169295, - -0.5545480533310743), - TestValue(1525, 0.98, 4.42, -1470.6581882505307, 949.7868112427125, - -209.3921137821333), - TestValue(10233, 0.98, 4.42, -9979.732399673303, 6382.717538120574, - -1436.6605223201655), - TestValue(0, 0.98, 800, -2.6600289899732843, -2.655611540996233, - -5.521811221598227e-6), - TestValue(6, 0.98, 800, -3.3605232422682434, 3.324471372446295, - -3.956077495042773e-6), - TestValue(14, 0.98, 800, -14.064685186393206, 11.297915257036331, - -0.0000880381480067513), - TestValue(1525, 0.98, 800, -7215.758505610536, 1517.2821289589797, - -0.8326540460538301), - TestValue(10233, 0.98, 800, -55549.328711222974, 10196.375797335235, - -10.124180638880098), - TestValue(0, 0.98, 15324, -2.6642246279711728, -2.663993040879676, - -1.5112705753494993e-8), - TestValue(6, 0.98, 15324, -3.363540258013437, 3.33496389215052, - -1.0907247577351953e-8), - TestValue(14, 0.98, 15324, -14.13194319100731, 11.33357313619078, - -2.435836623743224e-7), - TestValue(1525, 0.98, 15324, -8092.446096303554, 1522.0708941042951, - -0.004625777692748342), - TestValue(10233, 0.98, 15324, -71394.82522420201, 10228.557056242118, - -0.15615587395017738), - TestValue(0, 5, 0.0004, -0.005129619482412584, -0.0003999989219313857, - -11.824051401202995), - TestValue(6, 5, 0.0004, -9.620038054348033, -0.0003838278927170233, - 2490.418269104288), - TestValue(14, 5, 0.0004, -10.466998764550176, -0.00036226652043120673, - 2491.261123184485), - TestValue(1525, 5, 0.0004, -15.159872938045737, 0.003710137670052396, - 2485.8065870493465), - TestValue(10233, 5, 0.0004, -17.08620434453769, 0.027179691403163724, - 2429.0366049190375), - TestValue(0, 5, 0.065, -0.5026973821845353, -0.06497154463642642, - -6.734243654586291), - TestValue(6, 5, 0.065, -4.885023351671648, -0.062344895691172866, - 10.80289048394233), - TestValue(14, 5, 0.065, -5.677757542844056, -0.058842697097501465, - 11.638894745962721), - TestValue(1525, 5, 0.065, -10.72300452714444, 0.6026350622821846, - 6.184059083508757), - TestValue(10233, 5, 0.065, -16.31585327378707, 4.414778231493505, - -50.560431225223354), - TestValue(0, 5, 4.42, -15.66097567784475, -4.292171719051572, - -2.5721275924871447), - TestValue(6, 5, 4.42, -11.003030464726871, -4.1186491657279145, - -1.685183812308864), - TestValue(14, 5, 4.42, -8.894227315595337, -3.887285761296371, - -1.146456971892401), - TestValue(1525, 5, 4.42, -37.685674017920974, 39.81147725071138, - -6.586785126365316), - TestValue(10233, 5, 4.42, -286.73197869332216, 291.6505429744463, - -61.66251014536436), - TestValue(0, 5, 800, -136.14280118600206, -125.18861230732865, - -0.013692736098342095), - TestValue(6, 5, 800, -113.72441622976723, -120.12752689963362, - -0.012542423482431009), - TestValue(14, 5, 800, -93.60340589585249, -113.3794130227069, - -0.011094843043363584), - TestValue(1525, 5, 800, -1473.1928123052458, 1161.170595481826, - -0.5543680944524236), - TestValue(10233, 5, 800, -16253.562993271455, 8506.492550516554, - -8.17868004560178), - TestValue(0, 5, 15324, -147.69907423135155, -146.9895630824659, - -0.000046300649236385993), - TestValue(6, 5, 15324, -124.33517719483288, -141.0471157696556, - -0.00004260878910443758), - TestValue(14, 5, 15324, -103.0192965724122, -133.1238526859085, - -0.000037924592620959194), - TestValue(1525, 5, 15324, -2121.414393156032, 1363.3824622568213, - -0.0037347843330450686), - TestValue(10233, 5, 15324, -30500.05087098574, 9987.854328915519, - -0.14991286444070076), - TestValue(0, 368, 0.0004, -0.15032961840434253, -0.0004000000000000001, - -374.8240460108563), - TestValue(6, 368, 0.0004, -9.765221882218956, -0.0004000000000000001, - 2127.4587020676695), - TestValue(14, 368, 0.0004, -10.61216103101976, -0.0004000000000000001, - 2128.35545957858), - TestValue(1525, 368, 0.0004, -15.300962794837687, -0.0004000000000000001, - 2133.081933919416), - TestValue(10233, 368, 0.0004, -17.20382461596455, -0.0004000000000000001, - 2134.9858361205343), - TestValue(0, 368, 0.065, -24.097668920590625, -0.065, -369.7333680090865), - TestValue(6, 368, 0.065, -28.477367666024247, -0.065, -352.1558238867848), - TestValue(14, 368, 0.065, -29.266598891792, -0.065, -351.2659396464002), - TestValue(1525, 368, 0.065, -33.65022328528903, -0.065, -346.5441943953198), - TestValue(10233, 368, 0.065, -35.430094188966905, -0.065, - -344.6403282546427), - TestValue(0, 368, 4.42, -1619.991262543284, -4.420000000000001, - -365.5138603039104), - TestValue(6, 368, 4.42, -1615.1572361520293, -4.420000000000001, - -364.58765802750503), - TestValue(14, 368, 4.42, -1612.8136580987152, -4.420000000000001, - -363.99658652545247), - TestValue(1525, 368, 4.42, -1597.2619947735784, -4.420000000000001, - -359.55031671340777), - TestValue(10233, 368, 4.42, -1590.7558162463538, -4.420000000000001, - -357.6488775415163), - TestValue(0, 368, 800, -289052.31061786565, -800., -360.3153882723321), - TestValue(6, 368, 800, -289018.7634915345, -800., -360.30791160295655), - TestValue(14, 368, 800, -288983.8041593674, -800., -360.29802888017133), - TestValue(1525, 368, 800, -287560.8055295305, -800., -359.2481146209497), - TestValue(10233, 368, 800, -286189.66239505477, -800., -357.6907741283057), - TestValue(0, 368, 15324, -5.491551922548806e6, -15324., - -357.36282449417945), - TestValue(6, 368, 15324, -5.4915006777682435e6, -15324., - -357.36243301536337), - TestValue(14, 368, 15324, -5.491442187376253e6, -15324., - -357.36191128189233), - TestValue(1525, 368, 15324, -5.486439246638204e6, -15324., - -357.2679504598954), - TestValue(10233, 368, 15324, -5.474353427491954e6, -15324., - -356.85132040461684), + {0, -36., 0.000100000000000000004792174, -2.31952283024087929523226e-16, + -2.31952283023818920215225e-16, -2.69009308000224930599517e-24}, + {0, -36., 0.0840000000000000052180482, -2.3195228302435661858205e-16, + -2.31952283024356298332874e-16, -3.81249019275872869665102e-30}, + {0, -36., 3.75999999999999978683718, -2.31952283024356931676723e-16, + -2.31952283024356924522221e-16, -1.9027933171192913158222e-33}, + {0, -36., 461., -2.31952283024356938772873e-16, + -2.3195228302435693871452e-16, -1.26580106437037714160219e-37}, + {0, -36., 142311., -2.31952283024356938831037e-16, + -2.31952283024356938830848e-16, -1.32828224194512041655272e-42}, + {0, -26.1000000000000014210855, 0.000100000000000000004792174, + -4.62289481781287824178924e-12, -4.62289471095709740559354e-12, + -1.06855780836195688701657e-15}, + {0, -26.1000000000000014210855, 0.0840000000000000052180482, + -4.62289492454145310046865e-12, -4.62289492441424382973675e-12, + -1.51439608014165778339631e-21}, + {0, -26.1000000000000014210855, 3.75999999999999978683718, + -4.62289492466582046196526e-12, -4.6228949246629785527253e-12, + -7.55826925521030645308565e-25}, + {0, -26.1000000000000014210855, 461., -4.62289492466863919207562e-12, + -4.62289492466861601294603e-12, -5.02801075764912450865117e-29}, + {0, -26.1000000000000014210855, 142311., -4.62289492466866229611911e-12, + -4.62289492466866222103301e-12, -5.27619828240266701257698e-34}, + {0, -4.5, 0.000100000000000000004792174, -0.000471930181119562098823563, + -0.0000991078594800272766579705, -3.728223216395348042993}, + {0, -4.5, 0.0840000000000000052180482, -0.0104333684327369739264569, + -0.0098114347031002672745465, -0.00740397297186555491995328}, + {0, -4.5, 3.75999999999999978683718, -0.0110926179127638047017042, + -0.0110762714687203261403217, -4.34745852220174529502355e-6}, + {0, -4.5, 461., -0.0111088626902796869123984, -0.0111087288444673011605548, + -2.90337987821587513149328e-10}, + {0, -4.5, 142311., -0.0111089961046503927849115, + -0.0111089956710585016382076, -3.04679112048052440606817e-15}, + {0, -0.20000000000000001110223, 0.000100000000000000004792174, + -0.000901046250479348295492963, -0.000099987787464060911106451, + -8.010584630152873459984}, + {0, -0.20000000000000001110223, 0.0840000000000000052180482, + -0.199467033424133669556274, -0.0761837159352978552652188, + -1.46765854153375960277075}, + {0, -0.20000000000000001110223, 3.75999999999999978683718, + -0.740730806872720638510927, -0.672332093234306416324947, + -0.0181911472442591026765043}, + {0, -0.20000000000000001110223, 461., -0.818004584479470198072977, + -0.81727927438862128024214, -1.57334076106055928597981e-6}, + {0, -0.20000000000000001110223, 142311., -0.818728397963215257316515, + -0.818726042857481443695725, -1.65490069904197201165846e-11}, + {0, 2., 0.000100000000000000004792174, -0.00112103539054129295184384, + -0.000099998646665483029681996, -10.2103674387580987323199}, + {0, 2., 0.0840000000000000052180482, -0.377012371003043676678775, + -0.0830558079711177822673545, -3.49948289323721281132084}, + {0, 2., 3.75999999999999978683718, -4.08687891675795882940335, + -2.49194646483517166950394, -0.424184162745422141042474}, + {0, 2., 461., -7.33046427302122464456484, -7.27249028826061546955612, + -0.000125757016834293221277066}, + {0, 2., 142311., -7.38886427869112561814738, -7.38867246509109332505865, + -1.34784802321881715908971e-9}, + {0, 13.4000000000000003552714, 0.000100000000000000004792174, + -0.00226103403721276985381512, -0.000099999999984848563673037, + -21.6103403722792118658158}, + {0, 13.4000000000000003552714, 0.0840000000000000052180482, + -1.33366284302251745111466, -0.0839999893091445105829689, + -14.876938734683010272656}, + {0, 13.4000000000000003552714, 3.75999999999999978683718, + -45.4042061406096035346291, -3.75997857962063060636234, + -11.0755924364332381088221}, + {0, 13.4000000000000003552714, 461., -3350.22538971724911576917, + -460.678224812138394142983, -6.26799818851433996014357}, + {0, 13.4000000000000003552714, 142311., -246124.815117870661913638, + -117068.49513614471120162, -0.906861170125471331885925}, + {0, 84.2999999999999971578291, 0.000100000000000000004792174, + -0.00935103403719761843271571, -0.000100000000000000004792174, + -92.5103403719761798459793}, + {0, 84.2999999999999971578291, 0.0840000000000000052180482, + -7.28926283233166137753023, -0.0840000000000000052180482, + -85.7769384801388205324283}, + {0, 84.2999999999999971578291, 3.75999999999999978683718, + -311.988184720169192006433, -3.75999999999999978683718, + -81.9755810425981940674846}, + {0, 84.2999999999999971578291, 461., -36034.8035021785437227894, -461., + -77.1666019570033486394563}, + {0, 84.2999999999999971578291, 142311., -1.03081876937799361823842e7, + -142311., -71.4342299174339030881957}, + {7, -36., 0.000100000000000000004792174, -198.683622924671001775538, + 6.99999999998376310823605, -59997.5501489646175798213}, + {7, -36., 0.0840000000000000052180482, -238.883518369596162432722, + 6.9999999999999804386908, -69.0960424923444586909594}, + {7, -36., 3.75999999999999978683718, -256.751996432901803890657, + 6.99999999999999933622166, -0.718640295127228617069051}, + {7, -36., 461., -260.47982082426917795879, 6.99999999999999976452566, + -0.0000978945853445233788324339}, + {7, -36., 142311., -260.525013799174310121032, 6.99999999999999976803631, + -1.03688150994232402700446e-9}, + {7, -26.1000000000000014210855, 0.000100000000000000004792174, + -129.383623248262034973091, 6.99999967639274733829224, + -59997.5469131006865090386}, + {7, -26.1000000000000014210855, 0.0840000000000000052180482, + -169.583518369986016957648, 6.9999999996101358613744, + -69.0960424877584835201958}, + {7, -26.1000000000000014210855, 3.75999999999999978683718, + -187.451996432915042522719, 6.99999999998677065175815, + -0.718640295124939781566091}, + {7, -26.1000000000000014210855, 461., -191.17982082427388076164, + 6.99999999999530690927387, -0.000097894585344371117922657}, + {7, -26.1000000000000014210855, 142311., -191.225013799178942958982, + 6.99999999999537687768416, -1.03688150994072626043362e-9}, + {7, -4.5, 0.000100000000000000004792174, -11.2192075332052297804221, + 0.0623507285386109391248531, 9374.22326367570748540949}, + {7, -4.5, 0.0840000000000000052180482, -19.2633991074236276185347, + 6.17256900670521084397013, -59.3698802916059729529247}, + {7, -4.5, 3.75999999999999978683718, -36.2837402011840741159725, + 6.96830301037142800168473, -0.713160409032598892878726}, + {7, -4.5, 461., -39.9910983681716742384461, 6.9887225919756817853728, + -0.0000975289772446567319934955}, + {7, -4.5, 142311., -40.0361233417087705744711, 6.98889045789915253328171, + -1.03304486908480829056641e-9}, + {7, -0.20000000000000001110223, 0.000100000000000000004792174, + -11.1577615044596219717133, 0.000754889728272497048508259, + 9985.89049108549401489409}, + {7, -0.20000000000000001110223, 0.0840000000000000052180482, + -5.10523818836474804913996, 0.575173289456547933675185, + 5.01538223526623988873864}, + {7, -0.20000000000000001110223, 3.75999999999999978683718, + -8.27174735895246138021499, 5.07598581829756986530093, + -0.403937269906561294876212}, + {7, -0.20000000000000001110223, 461., -10.7102463026561453773833, + 6.1703108450892087653941, -0.0000725484455802695111298855}, + {7, -0.20000000000000001110223, 142311., -10.7437824687875434412923, + 6.18123368560834338804528, -7.70447935296972460390227e-10}, + {7, 2., 0.000100000000000000004792174, -11.1572212979358090452568, + -5.26523047722532179308571e-6, 9991.29214927237179177473}, + {7, 2., 0.0840000000000000052180482, -4.67822832418615765300115, + -0.00437313889759920793061993, 9.8011095064000150019666}, + {7, 2., 3.75999999999999978683718, -2.44742652447511137215548, + -0.131208500432565998079953, 0.0910218281904859865055242}, + {7, 2., 461., -1.92159366542738408699561, -0.382918557279757136122044, + 0.0000158891115286490464064998}, + {7, 2., 142311., -1.91424152166072074784582, -0.38903589945144240315715, + 1.69074180881189252797788e-10}, + {7, 13.4000000000000003552714, 0.000100000000000000004792174, + -11.1582665635858573490992, -0.0000999989393839702240949489, + 9980.83949989472447227895}, + {7, 13.4000000000000003552714, 0.0840000000000000052180482, + -5.55575146449510519129741, -0.083999098404519957659751, + -0.639658499701806072156294}, + {7, 13.4000000000000003552714, 3.75999999999999978683718, + -40.8853097501782241578175, -3.75993870125478372781648, + -9.93254120984925579863814}, + {7, 13.4000000000000003552714, 461., -3315.77631190381925371792, + -460.673338854838977134308, -6.25292229992679818636882}, + {7, 13.4000000000000003552714, 142311., -246051.646530372853518603, + -117067.253506656850199859, -0.906820707877630374849315}, + {7, 84.2999999999999971578291, 0.000100000000000000004792174, + -11.1653565625252413192581, -0.000100000000000000004792174, + 9909.93951050103628769457}, + {7, 84.2999999999999971578291, 0.0840000000000000052180482, + -11.5113505628995678711385, -0.0840000000000000052180482, + -71.5396476391501811780814}, + {7, 84.2999999999999971578291, 3.75999999999999978683718, + -307.469248451258373600901, -3.75999999999999978683718, + -80.8325192100658482257719}, + {7, 84.2999999999999971578291, 461., -36000.3495367018363608841, -461., + -77.151515469809951297333}, + {7, 84.2999999999999971578291, 142311., -1.03081131584031573940356e7, + -142311., -71.4341807304248851244502}, + {11, -36., 0.000100000000000000004792174, -306.294198663985101325787, + 10.9999999999744850169151, -99997.0711864556783212622}, + {11, -36., 0.0840000000000000052180482, -373.387720780818404406357, + 10.999999999999969393344, -116.240975127892379549188}, + {11, -36., 3.75999999999999978683718, -405.018198139691527464642, + 10.9999999999999990894639, -1.45345494122108727163354}, + {11, -36., 461., -413.383897627979795451303, 10.9999999999999997625131, + -0.000254934049399102214755753}, + {11, -36., 142311., -413.501921377875076957036, 10.9999999999999997680298, + -2.71559116633421423577824e-9}, + {11, -26.1000000000000014210855, 0.000100000000000000004792174, + -197.394199172482654828876, 10.9999994914769589000084, + -99997.0661015266437808503}, + {11, -26.1000000000000014210855, 0.0840000000000000052180482, + -264.487720781428391423827, 10.9999999993899980078309, + -116.240975120685847137987}, + {11, -26.1000000000000014210855, 3.75999999999999978683718, + -296.118198139709689507612, 10.999999999981852678434, + -1.45345494121749053012889}, + {11, -26.1000000000000014210855, 461., -304.483897627984544048369, + 10.9999999999952667973873, -0.000254934049398862947611818}, + {11, -26.1000000000000014210855, 142311., -304.601921377879715609259, + 10.9999999999953767477464, -2.71559116633170345973834e-9}, + {11, -4.5, 0.000100000000000000004792174, -11.706990517292534287664, + 0.0980363493375200627828593, 9017.8460181027726944107}, + {11, -4.5, 0.0840000000000000052180482, -28.2644285868714287029341, + 9.70535783036710290753852, -100.952775113605121643751}, + {11, -4.5, 3.75999999999999978683718, -58.5617425653278014914233, + 10.9565197428515127604419, -1.44484120738179930405144}, + {11, -4.5, 461., -66.8952715611464156307853, 10.9886262038729098348204, + -0.0002543593564776044431017}, + {11, -4.5, 142311., -67.0130312326551432919643, 10.9888901456535588389502, + -2.70956041824890637536164e-9}, + {11, -0.20000000000000001110223, 0.000100000000000000004792174, + -11.6101872629383735675807, 0.00124339116583624445400238, + 9981.48443912601296955335}, + {11, -0.20000000000000001110223, 0.0840000000000000052180482, + -5.90787076264096269638241, 0.947377292537602670212559, + 1.05849718208658084906379}, + {11, -0.20000000000000001110223, 3.75999999999999978683718, + -14.1259605624152918587705, 8.36073891060149916908714, + -0.948526674591890563801079}, + {11, -0.20000000000000001110223, 461., -20.4214207600281899855949, + 10.1632194847908259343291, -0.000214205349334668674488331}, + {11, -0.20000000000000001110223, 142311., -20.520713059859749302951, + 10.1812106733031004347544, -2.2874532593255647910321e-9}, + {11, 2., 0.000100000000000000004792174, -11.6092126588923464290597, + 0.0000488681502017790827148631, 9991.22977788173817623571}, + {11, 2., 0.0840000000000000052180482, -5.13540078317237206980856, + 0.0405883862872685488332284, 9.73996823769877214616492}, + {11, 2., 3.75999999999999978683718, -3.06137175973074859416703, + 1.21778462208320867130519, 0.0612622027041340431507753}, + {11, 2., 461., -2.88927536521627684556178, 3.55397957470930476869743, + -4.26994469309227722196226e-6}, + {11, 2., 142311., -2.89135678253022200706033, 3.61075642377121526650085, + -5.03190674865121337520085e-11}, + {11, 13.4000000000000003552714, 0.000100000000000000004792174, + -11.6102037914014672054394, -0.0000999983333263254586217559, + 9981.31845625030445310433}, + {11, 13.4000000000000003552714, 0.0840000000000000052180482, + -5.96770830536101890173932, -0.0839985893161630702750551, + -0.165549576777919563836645}, + {11, 13.4000000000000003552714, 3.75999999999999978683718, + -39.8538584150632717115069, -3.759915913617156940076, + -9.60353212925099392334374}, + {11, 13.4000000000000003552714, 461., -3300.14958948598756276581, + -460.670546879239310272208, -6.24440860614920345254555}, + {11, 13.4000000000000003552714, 142311., -246013.941380067564168339, + -117066.544004092358198853, -0.906797587679355764396099}, + {11, 84.2999999999999971578291, 0.000100000000000000004792174, + -11.617293789734793530787, -0.000100000000000000004792174, + 9910.41847291719271617468}, + {11, 84.2999999999999971578291, 0.0840000000000000052180482, + -11.9233068946770922978236, -0.0840000000000000052180482, + -71.065532655650617438992}, + {11, 84.2999999999999971578291, 3.75999999999999978683718, + -306.43777432844088456675, -3.75999999999999978683718, + -80.503504068925664332461}, + {11, 84.2999999999999971578291, 461., -35984.7200213335603843011, -461., + -77.1429957196861533815963}, + {11, 84.2999999999999971578291, 142311., -1.03080746722304058304262e7, + -142311., -71.4341526246487950815918}, + {2338, -36., 0.000100000000000000004792174, -62651.1907684420448777443, + 2337.99999999457695539095, -2.33699916660572733525724e7}, + {2338, -36., 0.0840000000000000052180482, -78386.4573396010374044213, + 2337.99999999999354376284, -27813.2247691071302770408}, + {2338, -36., 3.75999999999999978683718, -87244.5784991555256038266, + 2337.99999999999985553814, -615.235652274809944062375}, + {2338, -36., 461., -97261.3541997290761812182, 2337.99999999999999859168, + -3.26705772262828347042341}, + {2338, -36., 142311., -99951.6902594622698079577, 2337.99999999999999976424, + -0.000133435966326975414974749}, + {2338, -26.1000000000000014210855, 0.000100000000000000004792174, + -39504.990876519910618993, 2337.99989191671703492836, + -2.33699905852787203745752e7}, + {2338, -26.1000000000000014210855, 0.0840000000000000052180482, + -55240.2573397297094689766, 2337.99999987132480170888, + -27813.2247675754145700052}, + {2338, -26.1000000000000014210855, 3.75999999999999978683718, + -64098.3784991584079601655, 2337.99999999712082169713, + -615.235652274045473004386}, + {2338, -26.1000000000000014210855, 461., -74115.1541997291075706003, + 2337.99999999997193170739, -3.26705772262823261527954}, + {2338, -26.1000000000000014210855, 142311., -76805.4902594622778290634, + 2337.99999999999530115643, -0.000133435966326974881320714}, + {2338, -4.5, 0.000100000000000000004792174, -37.9188749421642944369079, + 20.858146249102902750828, -198577.847904345129887008}, + {2338, -4.5, 0.0840000000000000052180482, -5029.86319434730944140515, + 2064.90525599567281088845, -24562.2210710608448107488}, + {2338, -4.5, 3.75999999999999978683718, -13604.4870759968535893564, + 2331.10160386314082116743, -613.403922615515723001624}, + {2338, -4.5, 461., -23614.4216481166468801215, 2337.93255242508532760098, + -3.26693551284037806579056}, + {2338, -4.5, 142311., -26304.7015509659310958491, 2337.98870849677942716161, + -0.000133434683874347561598764}, + {2338, -0.20000000000000001110223, 0.000100000000000000004792174, + -17.2530056900344775246763, 0.285429102468546297600207, + 7145.03240130482692470614}, + {2338, -0.20000000000000001110223, 0.0840000000000000052180482, + -238.089236939508565582574, 217.47705608494119565083, + -2571.27861575441500369895}, + {2338, -0.20000000000000001110223, 3.75999999999999978683718, + -4005.5119497678293099549, 1919.26585035841237164672, + -504.067189818768777258691}, + {2338, -0.20000000000000001110223, 461., -13565.92078287865980031, + 2333.03782063120661396226, -3.25806818947358951236853}, + {2338, -0.20000000000000001110223, 142311., -16252.1224385913391334784, + 2337.1678232647280123578, -0.000133341466693716057784701}, + {2338, 2., 0.000100000000000000004792174, -16.9993203274403325176464, + 0.0315409623602125914052141, 9681.71391098754441853037}, + {2338, 2., 0.0840000000000000052180482, -36.345344890082340900417, + 26.196955662584086046202, -296.248198078484197433633}, + {2338, 2., 3.75999999999999978683718, -945.921470460610089695847, + 785.994533645635122586109, -203.555176892467693501323}, + {2338, 2., 461., -8461.86172575984929863877, 2293.84446785934606789742, + -3.18717688126682423770834}, + {2338, 2., 142311., -11115.2005139685862037106, 2330.48994045855231459004, + -0.000132584343734508495425318}, + {2338, 13.4000000000000003552714, 0.000100000000000000004792174, + -16.9688195052148641926364, -0.0000996457592914831445917036, + 9986.7200057158703803923}, + {2338, 13.4000000000000003552714, 0.0840000000000000052180482, + -10.8734665707860518042194, -0.0837024271645438342281827, + 5.22808308495811843394431}, + {2338, 13.4000000000000003552714, 3.75999999999999978683718, + -25.5045023028520732667892, -3.74665920542777317204743, + -4.50627645969873246867534}, + {2338, 13.4000000000000003552714, 461., -2105.32744445484601087282, + -459.046315074133113245574, -4.46701233139860018126931}, + {2338, 13.4000000000000003552714, 142311., -234622.84254414070827656, + -116653.790887199136613496, -0.893479868994494336336804}, + {2338, 84.2999999999999971578291, 0.000100000000000000004792174, + -16.9759091509741556489599, -0.000100000000000000004792174, + 9915.82354812310706660278}, + {2338, 84.2999999999999971578291, 0.0840000000000000052180482, + -16.8287689979316593747652, -0.0840000000000000052180482, + -65.668374254014350440872}, + {2338, 84.2999999999999971578291, 3.75999999999999978683718, + -292.075161470279026145921, -3.75999999999999978683718, + -75.4027226791102688967341}, + {2338, 84.2999999999999971578291, 461., -34788.2730773814556666438, -461., + -75.3620761655318490321712}, + {2338, 84.2999999999999971578291, 142311., -1.02962292135863589242571e7, + -142311., -71.4179345460698058518492}, + {9611, -36., 0.000100000000000000004792174, -257493.798714111504824811, + 9610.99999997770706584663, -9.60999902521143419375755e7}, + {9611, -36., 0.0840000000000000052180482, -322200.978630068417237347, + 9610.999999999973460561, -114395.144355373590193239}, + {9611, -36., 3.75999999999999978683718, -358701.177684561134391123, + 9610.99999999999940687087, -2548.13160475324762913113}, + {9611, -36., 461., -403079.072217515117383779, 9610.99999999999999493227, + -17.7630043020635701626676}, + {9611, -36., 142311., -424212.312866130535985746, 9610.99999999999999975238, + -0.00218253709318427116222024}, + {9611, -26.1000000000000014210855, 0.000100000000000000004792174, + -162344.899158395651111601, 9610.99955569358470701868, + -9.60999858092731645759248e7}, + {9611, -26.1000000000000014210855, 0.0840000000000000052180482, + -227052.078630597345206472, 9610.99999947105914950335, + -114395.144349077046283779}, + {9611, -26.1000000000000014210855, 3.75999999999999978683718, + -263552.277684572968739345, 9610.99999998817871670057, + -2548.13160475010505798557}, + {9611, -26.1000000000000014210855, 461., -307930.172217515232038494, + 9610.99999999989899826967, -17.7630043020633611084385}, + {9611, -26.1000000000000014210855, 142311., -329063.412866130554578654, + 9610.999999999995064897, -0.00218253709318426896848674}, + {9611, -4.5, 0.000100000000000000004792174, -104.508893420504154556264, + 85.743526266719416841998, -847430.234306281205498662}, + {9611, -4.5, 0.0840000000000000052180482, -20648.2403016158121274787, + 8488.39853461890805547167, -101030.965402842249733567}, + {9611, -4.5, 3.75999999999999978683718, -55983.0428066363797886683, + 9582.67667769505493377765, -2540.60175643222855702327}, + {9611, -4.5, 461., -100332.814925682181363541, 9610.75729475722022850913, + -17.7625019237987339653397}, + {9611, -4.5, 142311., -121465.824725374770167655, 9610.98814075622869244341, + -0.00218253182129617637581483}, + {9611, -0.20000000000000001110223, 0.000100000000000000004792174, + -19.5547487056230065372451, 1.17364684131883001763989, + -1735.73121296897564419702}, + {9611, -0.20000000000000001110223, 0.0840000000000000052180482, + -949.731171379774212295606, 894.23698468706897035991, + -10626.5306853797828683123}, + {9611, -0.20000000000000001110223, 3.75999999999999978683718, + -16521.513038999314726487, 7891.768160440031828256, + -2091.0860971061479068882}, + {9611, -0.20000000000000001110223, 461., -59023.1441094345908559896, + 9593.14395476867203137831, -17.7260454286430745149716}, + {9611, -0.20000000000000001110223, 142311., -80139.3868875039743196636, + 9610.1259811407200125367, -0.00218214857464869223866328}, + {9611, 2., 0.000100000000000000004792174, -18.5112200482628704439815, + 0.129968981779812349901792, 8698.84749102058108657733}, + {9611, 2., 0.0840000000000000052180482, -119.854424704305068624943, + 107.948248829969884782069, -1268.06413157597183571962}, + {9611, 2., 3.75999999999999978683718, -3933.80532649936636751687, + 3238.80127865994241519564, -854.484838151306300912966}, + {9611, 2., 461., -38021.2293458402145037005, 9452.10949634845787633542, + -17.4342406591600541232224}, + {9611, 2., 142311., -59002.2007387401537465017, 9603.1123321581496224457, + -0.0021790320685329142613448}, + {9611, 13.4000000000000003552714, 0.000100000000000000004792174, + -18.3822916425459253100299, -0.0000985437949788883229584548, + 9988.12276030177863802081}, + {9611, 13.4000000000000003552714, 0.0840000000000000052180482, + -12.1692487377720201796787, -0.0827767772596333470047606, + 6.63081050986194646155585}, + {9611, 13.4000000000000003552714, 3.75999999999999978683718, + -21.6460442652836866611064, -3.70522558331286636290617, + -3.10473788018826412150356}, + {9611, 13.4000000000000003552714, 461., -1491.91977767158488399654, + -453.969805440038841232377, -3.19739819420502413093782}, + {9611, 13.4000000000000003552714, 142311., -212175.813403243944037729, + -115363.73784931155578397, -0.853487617294708440971187}, + {9611, 84.2999999999999971578291, 0.000100000000000000004792174, + -18.3893801863409040880501, -0.000100000000000000004792174, + 9917.23732235214127244709}, + {9611, 84.2999999999999971578291, 0.0840000000000000052180482, + -18.1236255149538125593602, -0.0840000000000000052180482, + -64.2546271873853975660471}, + {9611, 84.2999999999999971578291, 3.75999999999999978683718, + -288.17526969257348869956, -3.75999999999999978683718, + -73.9901645192500912911873}, + {9611, 84.2999999999999971578291, 461., -34169.7871284528721914203, -461., + -74.0814500768543374677113}, + {9611, 84.2999999999999971578291, 142311., -1.02723620903825239882665e7, + -142311., -71.3688772675071099858592}, }; } // namespace neg_binomial_2_log_test_internal @@ -457,16 +512,13 @@ TEST(ProbDistributionsNegBinomial2Log, derivativesPrecomputed) { TEST(ProbDistributionsNegBinomial2Log, derivativesComplexStep) { using boost::math::differentiation::complex_step_derivative; using stan::math::is_nan; + using stan::math::log1p_exp; + using stan::math::log_sum_exp; using stan::math::neg_binomial_2_log_lpmf; using stan::math::var; - // TODO(martinmodrak) Reduced span of test, as the quick fix is not stable - // enough should be resolved once PR #1497 is merged - // std::vector n_to_test = {0, 7, 100, 835, 14238, 385000, 1000000}; - // std::vector eta_to_test = - // {-1531,-831, -124.5, -13, -2, 0, 0.536844, 1.26845, 11, 850, 2423}; - std::vector n_to_test = {0, 7, 100, 835}; - std::vector eta_to_test = {-2, 0, 0.536844, 1.26845, 11}; + std::vector n_to_test = {0, 7, 100, 835, 14238, 385000, 1000000}; + std::vector eta_to_test = {-124.5, -13, -2, 0, 0.536844, 1.26845, 11}; auto nb2_log_for_test = [](int n, const std::complex& eta, const std::complex& phi) { @@ -548,15 +600,88 @@ TEST(ProbDistributionsNegBinomial2Log, derivativesComplexStep) { } } +TEST(ProbDistributionsNegBinomial2Log, derivativesZeroOne) { + using stan::math::log1p_exp; + using stan::math::log_diff_exp; + using stan::math::log_sum_exp; + using stan::math::var; + using stan::test::expect_near_rel; + using stan::test::relative_tolerance; + + std::vector eta_to_test = {-943, -130, -15, -1, 0, 0.138, 14, 611}; + double phi_start = 1e-8; + double phi_max = 1e20; + for (double eta_dbl : eta_to_test) { + for (double phi_dbl = phi_start; phi_dbl < phi_max; + phi_dbl *= stan::math::pi()) { + std::stringstream msg; + msg << std::setprecision(20) << ", eta = " << eta_dbl + << ", phi = " << phi_dbl; + + var eta0(eta_dbl); + var phi0(phi_dbl); + var val0 = neg_binomial_2_log_lpmf(0, eta0, phi0); + + std::vector x0; + x0.push_back(eta0); + x0.push_back(phi0); + + std::vector gradients0; + val0.grad(x0, gradients0); + + var eta1(eta_dbl); + var phi1(phi_dbl); + var val1 = neg_binomial_2_log_lpmf(1, eta1, phi1); + + std::vector x1; + x1.push_back(eta1); + x1.push_back(phi1); + + std::vector gradients1; + val1.grad(x1, gradients1); + + double expected_value_0 = phi_dbl * (-log1p_exp(eta_dbl - log(phi_dbl))); + expect_near_rel("value, n = 0 " + msg.str(), val0.val(), + expected_value_0); + + double expected_deta_0 = -phi_dbl / (1 + phi_dbl / exp(eta_dbl)); + expect_near_rel("deta, n = 0 " + msg.str(), gradients0[0], + expected_deta_0); + + double expected_dphi_0 = 1.0 / (1.0 + phi_dbl / exp(eta_dbl)) + - log1p_exp(eta_dbl - log(phi_dbl)); + expect_near_rel("dphi, n = 0 " + msg.str(), gradients0[1], + expected_dphi_0); + + double expected_value_1 + = (phi_dbl + 1) * (-log1p_exp(eta_dbl - log(phi_dbl))) + eta_dbl; + expect_near_rel("value, n = 1 " + msg.str(), val1.val(), + expected_value_1); + + double expected_deta_1 + = exp(log(phi_dbl) - log_sum_exp(eta_dbl, log(phi_dbl))) + + expected_deta_0; + expect_near_rel("deta, n = 1 " + msg.str(), gradients1[0], + expected_deta_1); + + double expected_dphi_1 + = (1 + phi_dbl) / (phi_dbl + (phi_dbl * phi_dbl) / exp(eta_dbl)) + - log1p_exp(eta_dbl - log(phi_dbl)); + expect_near_rel("dphi, n = 1 " + msg.str(), gradients1[1], + expected_dphi_1); + } + } +} + TEST(ProbDistributionsNegBinomial2Log, derivatives_diff_sizes) { using stan::math::neg_binomial_2_log_lpmf; using stan::math::var; int N = 100; - double mu_dbl = 1.5; + double eta_dbl = 1.5; std::vector phi_dbl{2, 4, 6, 8}; - var mu(mu_dbl); + var mu(eta_dbl); std::vector phi; for (double i : phi_dbl) { phi.push_back(var(i)); @@ -568,8 +693,8 @@ TEST(ProbDistributionsNegBinomial2Log, derivatives_diff_sizes) { val.grad(x, gradients); double eps = 1e-6; - double grad_diff = (neg_binomial_2_log_lpmf(N, mu_dbl + eps, phi_dbl) - - neg_binomial_2_log_lpmf(N, mu_dbl - eps, phi_dbl)) + double grad_diff = (neg_binomial_2_log_lpmf(N, eta_dbl + eps, phi_dbl) + - neg_binomial_2_log_lpmf(N, eta_dbl - eps, phi_dbl)) / (2 * eps); EXPECT_FLOAT_EQ(grad_diff, gradients[0]); }