ClearAll[i, j, x, y, xv, yv, myint, mycount, maxvals, f]; 
    f[x_, y_] := (1 + Sin[x^2 + y])/8; 
    myint = Integrate[f[x, y], {x, -1, 1}, {y, -1, 1}] // 
        N

Plot3D[f[x, y], {x, -1, 1}, {y, -1, 1}]

Print["Value of the integral is ", 
        myint];
 For[j = 0, 
      j < 9, j++, 
mycount = 
        0; 
maxvals = 10^j; 
For[i = 0, 
        i < maxvals, i++, 
xv = 
          Random[Real, {-1, 1}]; 
yv = 
          Random[Real, {-1, 1}]; 
If[
          Random[Real, {0, 1/4}] < 
            f[xv, yv], mycount++]]; 
approx = 
        N[mycount/maxvals]; 
Print[approx, "\t",  
        approx - myint];
]