n_d_Monte_Carlo_integral.py

# function [value]  =  n_d_Monte_Carlo_integral (n, m, n_cases, f)
#  Accompanying program for the text
#
#     Classical and Modern Numerical Analysis:
#     Theory, Methods and Practice
#  by Azmy S. Ackleh, Edward J. Allen,
#     R. Baker Kearfott, and Padmanabhan Seshaiyer
#
#     (Taylor and Francis / CRC Press, 2009)
#
#  [value]  =  n_d_Monte_Carlo_integral (n, m, n_cases, f) returns a guess
#  for the integral of a function f of n variables over the unit n-cube
#  [0,1]^n, using m sample points.  The mean of these tries is returned in
#  value and the estimated standard deviation is returned in deviation.
#  This routine also prints the estimated value and estimated
#  for each try standard deviation.
#
#  The string "f" is the name of an m-file containing the function to be
#  integrated.  For example, if sum_of_squares.m contains
#
#     function [y]  =  sum_of_squares(x)
#     y  =  x'*x;
#
#  then issuing
#
#     value  =  n_d_Monte_Carlo_integral (10, 5000, 5, 'sum_of_squares')
#
#  will return an approximation to 10/3 in value by guessing at a
#  10-dimensional integral, using 5000 sample points in each of 5 tries,
#  resulting in a standard deviation of deviation.
#
#  The state of the random nuber generator is reset at each call of this
#  function, so the results are reproducible.

import numpy as np

def n_d_Monte_Carlo_integral (n, m, n_cases, f):

        np.random.seed(1)  # remove if you want different random numbers for each run

        value  =  0
        print('  try no.       m  est. mean   est. dev. ')
        for icase in range(n_cases):
                s1 = 0.0
                s2 = 0.0
                for i in range(m):
                        x = np.random.rand(n)
                        fval  =  f(x)
                        s1 += fval/m
                        s2 += fval**2/m
                estmean = s1
                value  += estmean
                eststdev = np.sqrt(s2-s1*s1)/np.sqrt(m)
                print( ' {:5.0f} {:10.0f} {:10.6f} {:10.6f}'.format( icase, m, estmean, eststdev ) )

        value  /= n_cases
        return value

if __name__ == '__main__':

        from sum_of_squares import sum_of_squares

        value = n_d_Monte_Carlo_integral (10, 5000, 5, sum_of_squares)
        print( value )


'''
  try no.       m  est. mean   est. dev.
     0       5000   3.329908   0.013245
     1       5000   3.329091   0.013095
     2       5000   3.338456   0.013254
     3       5000   3.332952   0.013198
     4       5000   3.330448   0.013404
3.3321708311896954
'''