RNG for Beta is Failing

[Peter Li]

Hi,

I'm finding that RNG for Beta is not passing the chi-squared goodness of fit test. I am using the following code:

    template <class RNG>

    inline double

    beta_rng(double alpha,

             double beta,

     RNG& rng1,

     RNG& rng2) {

      using boost::variate_generator;

      using boost::random::gamma_distribution;

      variate_generator<RNG&, gamma_distribution<> >

        rng_gamma_alpha(rng1, gamma_distribution<>(alpha, 1.0));

      variate_generator<RNG&, gamma_distribution<> >

        rng_gamma_beta(rng2, gamma_distribution<>(beta, 1.0));

      return rng_gamma_alpha() / (rng_gamma_alpha() + rng_gamma_beta());

    }

and test:

TEST(ProbDistributionsBeta, doesFit) {

  boost::random::mt19937 rng1, rng2;

  int N = 10000;

  int K = 5;

  boost::math::beta_distribution<>dist (2.0,2.0);

  boost::math::chi_squared mydist(K-1);

  double loc[4];

  for(int i = 1; i < K; i++)

    loc[i - 1] = quantile(dist, 0.2 * i);

  int count = 0;

  int bin [5] = {0, 0, 0, 0, 0};

  while (count < N) {

    double a = stan::prob::beta_rng(2.0,2.0,rng1,rng2);

    int i = 0;

    while (i < K-1 && a > loc[i]) 

++i;

    ++bin[i];

    count++;

   }

  double chi = 0;

  double expect [5] = {N / K, N / K, N / K, N / K, N / K};

  for(int j = 0; j < K; j++)

    chi += ((bin[j] - expect[j]) * (bin[j] - expect[j]) / expect[j]);

  EXPECT_TRUE(chi < quantile(complement(mydist, 1e-6)));

}

I'm not sure what is going wrong; I took the beta dist. rv generation from wiki (http://en.wikipedia.org/wiki/Beta_distribution#Generating_beta-distributed_random_variates).

When I run the test, I get that the value of chi is 2120.29 which is not less than 33.3768 (qunatile(complement(mydist,1e-6))). 

However, when I decrease N to 100 or less, the test passes..

Thanks,

Peter

[Peter Li]

And I just realized that the same thing is happening for double exponential in that the chi squared goodness of fit test fails for N>100.

[Bob Carpenter]

You should only have a single rng for an argument
and use it twice.

My guess is that you're making two copies and that's
why it's failing. If you copy an RNG in Boost, it'll
create another one that generates the same random
number sequence.

Otherwise, you should just break things down and see
what's going on step by step in the algorithm.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Bob Carpenter]

On 3/2/13 7:34 PM, Peter Li wrote:

> I'm finding that RNG for Beta is not passing the chi-squared goodness of fit test. I am using the following code:
>
>
> template <class RNG>
> inline double
> beta_rng(double alpha,
> double beta,
> RNG& rng1,
> RNG& rng2) {

You just wnat to use a single rng argument and pass it into
the constructor for rng_gamma_alpha and rng_gamma_beta.
The types look OK.

> using boost::variate_generator;
> using boost::random::gamma_distribution;
> variate_generator<RNG&, gamma_distribution<> >
> rng_gamma_alpha(rng1, gamma_distribution<>(alpha, 1.0));
> variate_generator<RNG&, gamma_distribution<> >
> rng_gamma_beta(rng2, gamma_distribution<>(beta, 1.0));
> return rng_gamma_alpha() / (rng_gamma_alpha() + rng_gamma_beta());

Also, instead of the last line, you need to do this:

double a = rng_gamma_alpha();
double b = rng_gamma_beta(0;

return a / (a + b);

You're drawing three random numbers here instead of two.

- Bob

[Ben Goodrich]

On Saturday, March 2, 2013 7:34:00 PM UTC-5, Peter Li wrote:

and test:

TEST(ProbDistributionsBeta, doesFit) {

  boost::random::mt19937 rng1, rng2;

  int N = 10000;

  int K = 5;

A general comment, probably unrelated to this particular bug, is that I don't know if K = 5 is so good. This link

http://www.itl.nist.gov/div898/handbook/prc/section2/prc211.htm

suggests setting K = round(2 * N^(2/5)), which would be 80 in this case.

Ben

[Peter Li]

I'll be sure to change that for all of the tests.

Thanks,

Peter

[Peter Li]

Ah. I see now. That explains why the other three distributions weren't working either.


Thanks,
Peter

[Bob Carpenter]

On 3/2/13 7:56 PM, Ben Goodrich wrote:

> A general comment, probably unrelated to this particular bug, is that I don't know if K = 5 is so good. This link
>
> http://www.itl.nist.gov/div898/handbook/prc/section2/prc211.htm
>
> suggests setting K = round(2 * N^(2/5)), which would be 80 in this case.

Random numerical algorithms are a tricky business!

Do you think at 79/80, we need to worry about using
complementary inverse CDFs? I guess we'll find out!

- Bob

[Ben Goodrich]

I would guess that empirical quantiles would be fine, but if the theoretical quantiles are already set up, then we might as well use them.

Ben
 

[Bob Carpenter]

I'm not sure what you mean by empirical quantiles.

We can compute the inverse CDFs and inverse complementary CDFs
for most of these distributions through Boost's library. My
question was whether 80 bins was enough to specialize the code
to use the complemnentary CDFs for p > 0.5.

- Bob

[Ben Goodrich]

On Sunday, March 3, 2013 4:18:00 PM UTC-5, Bob Carpenter wrote:

We can compute the inverse CDFs and inverse complementary CDFs
for most of these distributions through Boost's library.  My
question was whether 80 bins was enough to specialize the code
to use the complemnentary CDFs for p > 0.5.

Oh. Well, in R they are the same to at least 13 digits of precision

qchisq(ppoints(80), df = 79) - qchisq(1 - ppoints(80), df = 79, lower.tail = FALSE)

so I doubt it would make any difference to the tests.

Ben