C++ (Qt)using bounds_t = std::vector<std::pair<double, double>>; using grid_t = std::vector<size_t>; class result_t{public: result_t() = default; result_t(const double & mean, const double & stddev) : m_mean(mean), m_stddev(stddev) {} double mean() const { return m_mean; } double stddev() const { return m_stddev; } double confidence_interval(const double & level = 0.997) const { using namespace boost::math; normal_distribution<double> dist; double t = quantile(complement(dist, (1.0-level)/2.0)); return t*m_stddev; } private: double m_mean; double m_stddev;}; template <class F, class Tol>inline result_t mc_integrate(const F & f, const bounds_t & bounds, const Tol & tol, size_t pps = 100){ static thread_local std::random_device rd; static thread_local std::mt19937 gen(rd()); size_t n = 0; double mean = 0.0; double variance = 0.0; std::vector<double> x(bounds.size()); double V = 1.0; for (const auto & range : bounds) V *= (range.second - range.first); do { for (size_t i = 0; i < pps; ++i, ++n) { for (size_t j = 0; j < bounds.size(); ++j) { x[j] = std::uniform_real_distribution<double>(bounds[j].first, bounds[j].second)(gen); } double d = f(x) - mean; mean += d/(n+1); variance += d*d*n/(n+1); } } while (!tol(V*sqrt(variance/(n*n)), n)); return result_t(V*mean, V*sqrt(variance/(n*n)));}
C++ (Qt)#include <iostream>#include <cmath>#include <vector> #include <specmath/monte.h> inline double func(const std::vector<double> & x){ return exp(-x[0]*x[0] - x[1]*x[1] - x[2]*x[2] - x[3]*x[3]);} int main(){ size_t point_per_step = 1000; size_t max_calls = 10000000; specmath::monte::bounds_t bounds{ {-10.0, 10.0}, {-10.0, 10.0}, {-10.0, 10.0}, {-10.0, 10.0} }; auto res = specmath::monte::mc_integrate(func, bounds, [&](const double & err, size_t n) { return (err < 1.e-2) || (n > max_calls); }, point_per_step); // We use the 3 sigma rule (p = 0.997) std::cout << res.mean() << " +- " << res.confidence_interval(0.997) << std::endl; std::cout << "exact solution = " << M_PI*M_PI << std::endl; return 0;}
C++ (Qt)#include <iostream>#include <cmath>#include <vector>#include <chrono> #include <specmath/monte.h> inline double func(const std::vector<double> & x){ return exp(-x[0]*x[0] - x[1]*x[1] - x[2]*x[2] - x[3]*x[3]);} int main(){ size_t point_per_step = 1000; size_t max_calls = 100000000000; specmath::monte::bounds_t bounds{ {-10.0, 10.0}, {-10.0, 10.0}, {-10.0, 10.0}, {-10.0, 10.0} }; auto start = std::chrono::high_resolution_clock::now(); auto res = specmath::monte::mc_integrate(func, bounds, [&](const double & err, size_t n) { return (err < 1.e-2) || (n > max_calls); }, point_per_step); auto stop = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop-start).count(); std::cout << "duration (ms) = " << duration << std::endl; // We use the 3 sigma rule (p = 0.997) std::cout << res.mean() << " +- " << res.confidence_interval(0.997) << std::endl; std::cout << "exact solution = " << M_PI*M_PI << std::endl; specmath::monte::grid_t grid{ 100, 100, 100, 100 }; specmath::monte::mc_integrator mc(100*100); start = std::chrono::high_resolution_clock::now(); res = mc.integrate(func, bounds, grid, 1.e-2, 100); stop = std::chrono::high_resolution_clock::now(); duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop-start).count(); std::cout << "duration (ms) = " << duration << std::endl; std::cout << res.mean() << " +- " << res.confidence_interval(0.997) << std::endl; std::cout << "exact solution = " << M_PI*M_PI << std::endl; return 0;}
