montecarlo1

from math import cos, sqrt
from ROOT import gRandom, TF2

ax = 0.01
ay = 0
bx = 1
by = 1
m = 10
N = 10000
n = 1000

def f(x, y):
    return 1.0 / (x ** 3 * (y + (cos(x)) ** 2 / 3))
def fdraw(x, p):
    return f(x[0], x[1])
def integrate(f, xlim, ylim):
    xleft, xright = xlim
    yleft, yright = ylim
    x = gRandom.Uniform(xleft, xright)
    y = gRandom.Uniform(yleft, yright)
    return (xright - xleft) * (yright - yleft) * f(x,y)

res = []

for i in range(0,m):
    cres = []
    for j in range(0, n):
    r = integrate(f, (ax,bx), (ay, by))
    res.append(r)
    cres.append(r)
    avg = sum(res) / len(res)
    cavg = sum(cres) / len(cres)
    var = (sum(r ** 2 for r in res) - avg ** 2) / len(res)
    cvar = (sum(r ** 2 for r in cres) - cavg ** 2) / len(cres)
    print "iteration number", i + 1
    print "cavg", cavg, "crms",  sqrt(cvar / len(cres))
    print "avg", avg, "rms", sqrt(var / len(res))


t = TF2("f", fdraw, ax, bx, ay, by, 1)
t.Draw("surf")
raw_input()