Fix issue #591 neg_binomial_rng

#### Summary:

Fix neg_binomial_rng behaviour

#### Intended Effect:

Fix issue #591 

#### How to Verify:

Tests should pass.

Now, to check that the test are correctly written, one can do the following:

In BDA (edition 2 at least), section standard probability distributions, we have "(...) under alternative parametrization, p = beta/(beta+1) (...)", which implies: beta = p/(p-1)

So, with BDA parametrization

neg_binomial_rng(alpha, beta)

should give the same as:

neg_binomial_rng(alpha, p/(1-p))

And this is exactly what is done in the tests now:

    TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) {
      
      ...
      
      boost::math::negative_binomial_distribution<>dist (5,0.6); //alpha=5, p=0.6
      boost::math::chi_squared mydist(K-1);

      ...

      while (count < N) {
        int a = stan::prob::neg_binomial_rng(5, 0.6/(1-0.6), rng); //alpha=5, beta=0.6/(1-0.6)=p/(1-p)
        
        ...
        
       }

       ...
    }

#### Side Effects:

None

#### Documentation:

None

#### Reviewer Suggestions:

Daniel Lee