diff options
| -rw-r--r-- | python/distribution-integrals.py | 74 |
1 files changed, 74 insertions, 0 deletions
diff --git a/python/distribution-integrals.py b/python/distribution-integrals.py new file mode 100644 index 0000000..95287fa --- /dev/null +++ b/python/distribution-integrals.py @@ -0,0 +1,74 @@ +# nsmc --- n-sphere Monte Carlo method +# Copyright © 2022 Arun I <arunisaac@systemreboot.net> +# Copyright © 2022 Murugesan Venkatapathi <murugesh@iisc.ac.in> +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see +# <https://www.gnu.org/licenses/>. + +import nsmc +from numpy import linspace +from statistics import mean +import matplotlib.pyplot as plt + +def mean_thunk(f, trials): + return mean([f() for i in range(trials)]) + +def samples_for_dimension_rtol(integrand, oracle, true_integral, n, rtol): + print(n, end=' ', flush=True) + integral, samples = nsmc.integral(integrand, oracle, true_integral(n), n, rtol) + return samples + +def average_samples_for_rtol(integrand, oracle, true_integral, rtol, trials, dimension): + for n in dimension: + yield mean_thunk(lambda: samples_for_dimension_rtol(integrand, oracle, true_integral, n, rtol), trials) + +def experiment(integrand, oracle, true_integral, rtols, trials, dimension, filename): + print(filename) + for rtol in rtols: + print(f'rtol = {rtol}:', end=' ', flush=True) + samples = list(average_samples_for_rtol(integrand, oracle, true_integral, rtol, trials, dimension)) + plt.plot(dimension, samples, '-x', label=f'rtol = {rtol}') + print() + plt.yscale('log') + plt.xlabel('Dimension') + plt.ylabel('Average number of samples') + plt.legend() + plt.savefig(filename) + plt.close() + +# Number of trials to average each dimension, rtol point over +trials = 2 +# Dimensions and relative error tolerances to plot +dimension = linspace(10, 100, 10, dtype=int) +rtols = [0.2, 0.1, 0.05] + +# uniform(0,1) domain +max_extent = 1.0 +oracle = nsmc.make_uniform_extent_oracle(0, max_extent) + +# Gaussian integrand +integrand = nsmc.gaussian_integrand() +true_integral = lambda n: nsmc.gaussian_true_integral(n, max_extent) +experiment(integrand, oracle, true_integral, rtols, trials, dimension, 'gaussian.png') + +# Polynomial integrand +coefficients = [1, -1.5, 0.6875, -0.09375] +integrand = nsmc.polynomial_integrand(coefficients) +true_integral = lambda n: nsmc.polynomial_true_integral(n, max_extent, coefficients) +experiment(integrand, oracle, true_integral, rtols, trials, dimension, 'polynomial.png') + +# x-coordinate integrand +integrand = nsmc.x_coordinate_integrand() +true_integral = lambda n: nsmc.x_coordinate_true_integral(n, max_extent) +experiment(integrand, oracle, true_integral, rtols, trials, dimension, 'x-coordinate.png') |