cálculo de integrais em pyton

Como posso realizar um programa capaz de calcular uma integral pelo método de monte carlo?

Iniciei o programa dessa maneira:

from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

#imagem
a = -1
b = 1

#função
xr = np.linspace(a,b,1000)
yr = xr

c = np.min(yr)
d = np.max(xr)

n = int(input("Insira o número de pontos que deseja gerar aleatoriamente: "))

#número de rodadas
r = int(input("Insira o número de rodadas: "))

xa = np.random.random(n)*(b-a)+a
xb = np.random.random(n)*(d-c)+c

if #adicionar a condição para o cálculo da integral

print("O número total de pontos gerados é: ",(len(xa)))
print("O número de pontos contidos na área desejada é: ",(j))

plt.plot(xr,yr)
plt.plot(xa,xb,"+")
plt.show()

preciso fazer de modo com que o programa se repita um numero r de vezes para poder calcular a integral dada pela média das estimativas, o desvio padrão das estimativas, o desvio padrão da média para a estimativa, lembrando que a integral nada mais é que a área desejada da função. Como posso criar a condição "if" para que ele realize o cálculo e repita o programa r vezes?