Vamos ao Casino?

Nas aulas abordámos o método de Monte Carlo e vimos como nos pode auxiliar a calcular o valor aproximado de pi. Eis o código:

def monte_carlo_pi(num_dardos):
    """
    Calcula o valor de pi pelo método de Monte Carlo.
    """
    # define e inicializa acumulador
    conta_dardos_in = 0
    for i in range(num_dardos):
        # gera posição dardo i
        x= random.random()
        y= random.random()
        # dentro ou fora?
        d = (x**2 + y**2)**0.5
        if d <= 1:
            conta_dardos_in = conta_dardos_in + 1
    res_pi = 4 * (conta_dardos_in/num_dardos)
    return res_pi

A ideia consiste em ter um quarto de círculo de raio um (cuja área é igual a pi/4) inscrito num quadrado de lado um. Simulamos de seguida o lançamento de dardos cuja posição de impacto é escolhida com base numa distribuição uniforme (todos os pontos são equi-prováveis). A percentagem de dardos que cai dentro do círculo, multiplicado por 4 dá-nos o valor pretendido.

Este método pode ser aplicado nas mais diversas situações e problemas. Por exemplo, determinar a probabilidade de um dardo lançado aleatoriamente cair dentro de uma certa zona. Procurámos resolver esta questão no caso das áreas 1 e 3 da figura abaixo.

A ideia para a solução não é diferente da anterior: contar a percentagem dos dardos que caem dentro de uma das áreas pretendidas:

def monte_carlo_reg_2(num_dardos):
    """
    Calcula o valor de uma área  pelo método de Monte Carlo.
    """
    # define e inicializa acumulador 
    conta_dardos_in = 0
    for i in range(num_dardos):
        # gera posição dardo i
        x= random.uniform(0,2)
        y= random.uniform(0,2)
        # dentro ou fora
        zona_1 = (x <= 1 and y >=1)
        zona_3 = (x >= 1 and y <= 1) or (x>= 1 and y <= x)
        if zona_1 or zona_3:
            conta_dardos_in = conta_dardos_in + 1
    res = conta_dardos_in/num_dardos
    return res

Com estes dois exemplos o leitor está preparado para dar o salto para uma generalização da abordagem. Imaginemos que queremos calcular um dado integral:

Sabemos que o valor do integral é dado pela área debaixo de f(x), para x entre a e b. Então conhecido f(x) podemos encontrar uma figura geométrica de área conhecida e que envolva f(x). Feito isto só temos que calcular a percentagem de dados que ficam debaixo da curva e multiplicar pela área. E o programa que o faz é o seguinte:

def monte_carlo_int(funcao, min_x, max_x, min_y, max_y, num_dardos):
    """
    Calcula o valor de um integral  pelo método de Monte Carlo.
    """
    # define e inicializa acumulador 
    conta_dardos_in = 0
    for i in range(num_dardos):
        # gera posição dardo i
        x= random.uniform(min_x,max_x)
        y= random.uniform(min_y,max_y)
        # dentro ou fora
        zona_limite = funcao(x) 
        if y < zona_limite:
            conta_dardos_in = conta_dardos_in + 1
    area = (max_x - min_x) * (max_y -min_y)
    res = area * (conta_dardos_in/num_dardos ) 
    return res

Podemos testar o programa para um caso concreto, o da função y = e^x. A imagem mostra o seu comportamento entre x= 0 e x = 1, e ainda alguns dardos que foram “lançados”:

Como chamamos o programa? Fácil. Definimos a função, os respectivos valores máximos e mínimos para x e y e o número de dardos:

def my_exp(x):
    return pow(math.e, x)

minx = 0
maxx = 1
miny = 0 
maxy = math.e

print(monte_carlo_int(my_exp, minx,maxx,miny,maxy,1000000))