hypervolumens.py

"""
Approximate unit hypersphere volumes by Monte Carlo sampling.
"""

from __future__ import division
from random import random

# Area of a unit circle, just to get started.
def area(trials):
    # Sprinkle `trials` random points (x,y) where 0<=x<=1, 0<=y<=1.
    # Count how many are at distance <= 1 from the origin.
    # The occupancy is the fraction of the total #points. It's also the estimated area of
    # the part of the circle within the [0..1,0..1] quadrant.
    # (r2 is the squared distance; we don't need to bother taking the sqrt.)
    occupancy = 1/trials * sum(r2 <= 1
                               for _ in xrange(trials)
                               for r2 in [sum(random()**2 for _ in xrange(2))])
    # The whole circle has 4 quadrants. Above, we didn't count the quadrants
    # with one or more negative coordinates.
    nquadrants = 2**2
    return occupancy * nquadrants

## area(10000)
#. 3.1616

# Hypervolume of a unit hypersphere
def hypervolume(ndims, trials):
    occupancy = 1/trials * sum(r2 <= 1
                               for _ in xrange(trials)
                               for r2 in [sum(random()**2 for _ in xrange(ndims))])
    ncells = 2**ndims
    return occupancy * ncells

for n in range(1, 11): print n, hypervolume(n, 1000000)

#. 1 2.0
#. 2 3.13884
#. 3 4.188072
#. 4 4.925312
#. 5 5.255552
#. 6 5.147648
#. 7 4.711424
#. 8 4.054272
#. 9 3.265024
#. 10 2.450432