double f = 1.0 - boost::math::ibetac(_r, _k + 1, _p);

NegativeBinomial dist(r, p);
double prob = dist.cdf(k);

struct NegativeBinomialSystem
{
   NegativeBinomialSystem(int k, double r, double p) : _k(k), _r(r), _p(p) { }

   // Return both f(p) and dot(f)(p).
   std::pair<double, double> operator()(double const& mu)
   {
        double f = 1.0 - boost::math::ibetac(_r, _k + 1, _p);

        double df = std::pow(1.0 - _p, _k) * std::pow(_p, _r - 1.0);
        df /= boost::math::beta(_r, _k + 1);

        return std::make_pair(f, df);
   }
   const int _k;
   const double _r;
   const double _p;
};

double solveNegativeBinomialParameter(int k, double r, double prob, double guess)
{
   if (k == 0)
      return -std::log(prob);
   else
   {
      return boost::math::tools::newton_raphson_iterate(
         NegativeBinomialSystem(k, r, prob), guess, 0.0001, 10.0, 17);
   }
}

#include <iostream>
#include <boost/test/unit_test.hpp>
#include <boost/test/output_test_stream.hpp>
#include "math/distributions/negative_binomial.h"

#include <cmath>
#include <boost/math/distributions/negative_binomial.hpp>
#include <boost/math/tools/roots.hpp>
#include "math/engine_math.h"
#include <boost/math/special_functions/beta.hpp>

BOOST_AUTO_TEST_SUITE(negative_binomial_distribution_test)

void solveNegativeBinomialParameterFromProbability_verifyCorrect(double r, double p, int k)
{
    NegativeBinomial dist(r, p);
    double prob = dist.cdf(k);

    double p2 = solveNegativeBinomialParameter(k, r, prob, 1.0);
    BOOST_CHECK_CLOSE_FRACTION(p, p2, 0.0001);
}

BOOST_AUTO_TEST_CASE(htSolveNegativeBinomialParameterFromProbability_multipleCases_verifyCorrect)
{
    const double _R = 8.618515;
    solveNegativeBinomialParameterFromProbability_verifyCorrect(_R, 0.5, 1);
}
BOOST_AUTO_TEST_SUITE_END()