| 1 | #include <boost/math/tools/roots.hpp>
|
|---|
| 2 | #include <boost/math/distributions/gamma.hpp>
|
|---|
| 3 |
|
|---|
| 4 | typedef boost::math::policies::policy<boost::math::policies::digits10<16> > acc_policy;
|
|---|
| 5 | typedef boost::math::gamma_distribution<double, acc_policy> gamma_dist;
|
|---|
| 6 | typedef std::pair<double, double> result_t;
|
|---|
| 7 |
|
|---|
| 8 | int main() {
|
|---|
| 9 | const long unsigned int u = 55;
|
|---|
| 10 | const long unsigned int n = 2373;
|
|---|
| 11 | const double nu_a = 6456.947055550609;
|
|---|
| 12 | const double alpha = 27.21005923114458;
|
|---|
| 13 | const double nu_m = ((double)u/n)*nu_a;
|
|---|
| 14 |
|
|---|
| 15 | const gamma_dist g_m(nu_m, 1.0/alpha);
|
|---|
| 16 | const gamma_dist g_a(nu_a, 1.0/alpha);
|
|---|
| 17 | auto zero_func = [&] (const double& x) {
|
|---|
| 18 | return boost::math::pdf(g_m, x) - boost::math::pdf(g_a, x);
|
|---|
| 19 | };
|
|---|
| 20 | const double lower = boost::math::mean(g_m);
|
|---|
| 21 | const double upper = 237.3;
|
|---|
| 22 | boost::uintmax_t max_iter = 100;
|
|---|
| 23 | boost::math::tools::eps_tolerance<double> tol(16); // 10 bit precition on result
|
|---|
| 24 | result_t res = boost::math::tools::toms748_solve(zero_func, lower, upper, tol, max_iter);
|
|---|
| 25 | //result_t res = boost::math::tools::bisect(zero_func, lower, upper, tol, max_iter); // TODO: This works as expected
|
|---|
| 26 | return 0;
|
|---|
| 27 | }
|
|---|
