# 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 '''