domingo, 6 de novembro de 2016

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:
01.def monte_carlo_pi(num_dardos):
02.    """
03.    Calcula o valor de pi pelo método de Monte Carlo.
04.    """
05.    # define e inicializa acumulador
06.    conta_dardos_in = 0
07.    for i in range(num_dardos):
08.        # gera posição dardo i
09.        x= random.random()
10.        y= random.random()
11.        # dentro ou fora?
12.        d = (x**2 + y**2)**0.5
13.        if d <= 1:
14.            conta_dardos_in = conta_dardos_in + 1
15.    res_pi = 4 * (conta_dardos_in/num_dardos)
16.    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:
01.def monte_carlo_reg_2(num_dardos):
02.    """
03.    Calcula o valor de uma área  pelo método de Monte Carlo.
04.    """
05.    # define e inicializa acumulador
06.    conta_dardos_in = 0
07.    for i in range(num_dardos):
08.        # gera posição dardo i
09.        x= random.uniform(0,2)
10.        y= random.uniform(0,2)
11.        # dentro ou fora
12.        zona_1 = (x <= 1 and y >=1)
13.        zona_3 = (x >= 1 and y <= 1) or (x>= 1 and y <= x)
14.        if zona_1 or zona_3:
15.            conta_dardos_in = conta_dardos_in + 1
16.    res = conta_dardos_in/num_dardos
17.    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:
01.def monte_carlo_int(funcao, min_x, max_x, min_y, max_y, num_dardos):
02.    """
03.    Calcula o valor de um integral  pelo método de Monte Carlo.
04.    """
05.    # define e inicializa acumulador
06.    conta_dardos_in = 0
07.    for i in range(num_dardos):
08.        # gera posição dardo i
09.        x= random.uniform(min_x,max_x)
10.        y= random.uniform(min_y,max_y)
11.        # dentro ou fora
12.        zona_limite = funcao(x)
13.        if y < zona_limite:
14.            conta_dardos_in = conta_dardos_in + 1
15.    area = (max_x - min_x) * (max_y -min_y)
16.    res = area * (conta_dardos_in/num_dardos )
17.    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:
01.def my_exp(x):
02.    return pow(math.e, x)
03. 
04.minx = 0
05.maxx = 1
06.miny = 0
07.maxy = math.e
08. 
09.print(monte_carlo_int(my_exp, minx,maxx,miny,maxy,1000000))

Sem comentários:

Enviar um comentário