diff options
| -rw-r--r-- | python/distribution-volumes.py | 77 |
1 files changed, 77 insertions, 0 deletions
diff --git a/python/distribution-volumes.py b/python/distribution-volumes.py new file mode 100644 index 0000000..57e6c03 --- /dev/null +++ b/python/distribution-volumes.py @@ -0,0 +1,77 @@ +# 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(oracle, true_volume, n, rtol): + print(n, end=' ', flush=True) + volume, samples = nsmc.volume(oracle(n), true_volume(n), n, rtol) + return samples + +def average_samples_for_rtol(oracle, true_volume, rtol, trials, dimension): + for n in dimension: + yield mean_thunk(lambda: samples_for_dimension_rtol(oracle, true_volume, n, rtol), trials) + +def experiment(oracle, true_volume, rtols, trials, dimension, filename): + print(filename) + for rtol in rtols: + print(f'rtol = {rtol}:', end=' ', flush=True) + samples = list(average_samples_for_rtol(oracle, true_volume, 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) +oracle = lambda n: nsmc.make_uniform_extent_oracle(0, 1) +true_volume = lambda n: nsmc.uniform_true_volume(0, 1, n) +experiment(oracle, true_volume, rtols, trials, dimension, 'uniform.png') + +# beta(2,2) +oracle = lambda n: nsmc.make_beta_extent_oracle(2, 2) +true_volume = lambda n: nsmc.beta_true_volume(2, 2, n) +experiment(oracle, true_volume, rtols, trials, dimension, 'beta.png') + +# arcsine +oracle = lambda n: nsmc.make_beta_extent_oracle(0.5, 0.5) +true_volume = lambda n: nsmc.beta_true_volume(0.5, 0.5, n) +experiment(oracle, true_volume, rtols, trials, dimension, 'arcsine.png') + +# cube (pretty slow for large dimensions and tight tolerances) +oracle = lambda n: nsmc.make_cube_extent_oracle(n, 1.0) +true_volume = lambda n: nsmc.cube_true_volume(n, 1.0) +dimension = linspace(10, 40, 4, dtype=int) +rtols = [0.2, 0.1] +experiment(oracle, true_volume, rtols, trials, dimension, 'cube.png') |