from scipy import * from numpy.random import * import math import matplotlib.pyplot as pyplot import numpy as np from scipy.integrate import tplquad def expf3d(x1, x2, x3): return 1.0*(np.exp(x1**2+x2**2+x3**2)) def mc3d(f, N_list): res = [] for N in N_list: res.append(sum(f(uniform(0,1,N), uniform(0,1,N), uniform(0,1,N)))/N) return array(res) def avg(f, N): res = f() for i in range(1,N): res += f() return res/N ########################################## pyplot.figure() print mc3d(expf3d, [1000000]) N_range = logspace(0, 6, 10) #print N_range #print mc2d(f2d, N_range) #val, abserr = tplquad(ff, bounds_x[0], bounds_x[1], lambda y: bounds_y[0], lambda y: bounds_y[1], \ # lambda y, z: bounds_z[0], lambda y, z: bounds_z[1]) val3d, abserr = tplquad(lambda x,y,z: expf3d(x, y, z), 0, 1, lambda y: 0, lambda y: 1, lambda y, z: 0, lambda y, z: 1) #print mc3d(expf3d, N_range) pyplot.xscale("log") pyplot.yscale("log") pyplot.plot(N_range, avg(lambda: abs(mc3d(expf3d, N_range) - val3d), 20), "rx", markersize=3) #pyplot.plot(N_range, avg(lambda: abs(mc2d(f2d, N_range) - pi), 20), "bo", markersize=3) pyplot.show() #figure.savefig("pi.pdf")