Replace exp(lbeta(a,b)) with beta(a,b)

[andrjohns]

Description

Several functions are calling exp(lbeta(a, b)) (i.e. the exponential of the log of the beta function), which could be replaced with the introduction of the beta function.

Example

Can be seen in prim/scal/fun/grad_reg_inc_beta and prim/scal/prob/beta_binomial_cdf (among others

Current Version:

v2.19.1

[andrjohns]

Doing some benchmarking comparing exp(lbeta(x,y)) with boost::math::beta(x,y), and it looks like stan's implementation is an order of magnitude faster than Boost's:

andrew@Lappy:~/perf-math$ ./beta_func 
2019-04-30 16:28:36
Running ./beta_func
Run on (8 X 4000 MHz CPU s)
CPU Caches:
  L1 Data 32K (x4)
  L1 Instruction 32K (x4)
  L2 Unified 256K (x4)
  L3 Unified 8192K (x1)
--------------------------------------------------
Benchmark           Time           CPU Iterations
--------------------------------------------------
Beta_Stan      171114 ns     171043 ns       4139
Beta_Boost    1571291 ns    1570855 ns        437

With benchmark:

#include <benchmark/benchmark.h>
#include <stan/math/prim/scal.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <Eigen/Dense>

static void Beta_Stan(benchmark::State& state) {
  using Eigen::Matrix;
  using Eigen::VectorXd;
  using namespace stan::math;
  VectorXd v1 = VectorXd::Random(1000).array().exp();
  VectorXd v2 = VectorXd::Random(1000).array().exp();

  for (auto _ : state) {
    VectorXd out(1000);
    for(int i = 0; i < 1000; i++)
      out[i] = exp(lbeta(v1[i],v2[i]));
  }
}
BENCHMARK(Beta_Stan);

static void Beta_Boost(benchmark::State& state) {
  using Eigen::Matrix;
  using Eigen::VectorXd;
  using namespace stan::math;
  VectorXd v1 = VectorXd::Random(1000).array().exp();
  VectorXd v2 = VectorXd::Random(1000).array().exp();

  for (auto _ : state) {
    VectorXd out(1000);
    for(int i = 0; i < 1000; i++)
      out[i] = boost::math::beta(v1[i],v2[i]);
  }
}
BENCHMARK(Beta_Boost);

BENCHMARK_MAIN();

So it may just be worth defining stan::math::beta using exp(lbeta()) and adding gradients for that

[wds15]

Surprising. What about the accuracy of the two functions?

[andrjohns]

Identical, as far as I can tell. They're both defined up to the same extreme with the same result:

#include <stan/math/prim/scal.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <cmath>
#include <iostream>

double beta_stan(double x) {
  using std::exp;
  return exp(stan::math::lbeta(x, x));
}
double beta_boost(double x) {
  return boost::math::beta(x, x);
}

int main() {
  std::cout.precision(18);
  std::cout << std::scientific;
  for (double x = 536; x < 537; x += 0.01)
    std::cout << "x = " << x << "  1: " << beta_stan(x)
              << "  2: " << beta_boost(x) << std::endl;
}

Returns:

x = 5.360000000000000000e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360099999999999909e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360199999999999818e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360299999999999727e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360399999999999636e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360499999999999545e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360599999999999454e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360699999999999363e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360799999999999272e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360899999999999181e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.360999999999999091e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.361099999999999000e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.361199999999998909e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.361299999999998818e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.361399999999998727e+02  1: 4.940656458412465442e-324  2: 4.940656458412465442e-324
x = 5.361499999999998636e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.361599999999998545e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.361699999999998454e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.361799999999998363e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.361899999999998272e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.361999999999998181e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.362099999999998090e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.362199999999997999e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.362299999999997908e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00
x = 5.362399999999997817e+02  1: 0.000000000000000000e+00  2: 0.000000000000000000e+00

[bob-carpenter]

Stan just calls the standard library functions for double values. That is, std::exp is the same as stan::math::exp and std::lgamma is the same as stan::math::lgamma.

If you look at Boost's doc for beta, you'll see that it's templated. The doc says:

There are effectively two versions of this function internally: a fully generic version that is slow, but reasonably accurate, and a much more efficient approximation that is used where the number of digits in the significand of T correspond to a certain Lanczos approximation. In practice any built-in floating-point type you will encounter has an appropriate Lanczos approximation defined for it. It is also possible, given enough machine time, to generate further Lanczos approximation's using the program libs/math/tools/lanczos_generator.cpp.

I have no idea what that number of digits thing is, but here's their implementation if anyone wants to dive in.

[andrjohns]

Closed via #1228