#!/net/python/bin/python
from numpy import *
import pylab as p

EXACT = 155.0/6.0

def error(measured,exact):
    return abs( (measured-exact)/exact )

def f(x):
    return (sum(x))**2

def SampleMean(f,N,xmin,xmax):
    assert len(xmin) == len(xmax)
    L = len(xmin)

    x = random.random((N,L))

    domain = exp( sum (log(xmax[i]-xmin[i]) for i in range(L) ) )
                #(x0M-x0m)(x1M-x1m)...(xnM-xnm)
    return domain * 1.0/N * sum(f(x[i]) for i in range(N))

if __name__ == "__main__":
    Nlist = [1]
    Elist = []
    while 2*Nlist[-1] < 200000:
        Nlist.append(2*Nlist[-1])

    xmin = [0,0,0,0,0,0,0,0,0,0]
    xmax = [1,1,1,1,1,1,1,1,1,1]
    
    for N in Nlist:
        s = SampleMean(f,N,xmin,xmax)
        Elist.append(error(s,EXACT))

    Nlist = array(Nlist)
    Elist = array(Elist)

    p.loglog( 1/sqrt(Nlist),Elist )
    p.xlabel('1/sqrt(N)')
    p.ylabel('abs(error)')
    p.title('error vs. N in 10D Sample Mean Monte Carlo Integration')
    p.show()