.. -*- coding: utf-8 -*- Monte Carlo =========== Os métodos de Monte Carlo são uma coleção de técnicas computacionais para a solução aproximada de problemas. Esses métodos fazem uso fundamental de amostras aleatórias. Duas classes de problemas estatísticos são mais comumente abordadas dentro desta estrutura: integração e otimização. A metodologia de Monte Carlo também é amplamente utilizada na simulação de sistemas físicos, químicos e biológicos. .. image:: ./imgs/Pi_30K.* :width: 500 :height: 500 :scale: 50 :alt: By nicoguaro - Own work, CC BY 3.0, :align: center :target: https://commons.wikimedia.org/w/index.php?curid=14609430 Muitos modelos interessantes têm estruturas extremamente complexas e não podem ser facilmente tratados usando técnicas tradicionais. Nesses casos o método de Monte Carlo é um arcabouço computacional para realizarmos inferência estatística. .. Dentro do paradigma bayesiano, todas as informações nas quais a inferência pode ser baseada são codificadas dentro da distribuição de probabilidade posterior. Usando os métodos de Monte Carlo, podemos caracterizar essas distribuições e calcular as expectativas sob elas: a técnica inferencial primária. Um Exemplo Simples Os métodos de Monte Carlo invertem o problema usual da estatística: em vez de estimar quantidades aleatórias de maneira determinística, quantidades aleatórias são empregadas para fornecer estimativas de quantidades determinísticas. Por exemplo, um experimento simples de Monte Carlo considera a chuva que cai uniformemente ao acaso (ou seja, a localização de qualquer gota de chuva pode ser interpretada como a realização de uma variável aleatória uniformemente distribuída) sobre uma região quadrada do espaço, e um círculo inscrito nesse quadrado . Sem referência a qualquer teoria de probabilidade formal, é intuitivo que a probabilidade de uma gota de chuva uniforme cair em qualquer região dentro do quadrado deve ser proporcional à área dessa região e independente de sua localização. Conseqüentemente, a probabilidade, p, de que uma gota de chuva esteja dentro do círculo inscrito pode ser expressa em termos de suas áreas. Se o quadrado tem lados de comprimento 2r, o círculo deve ter raio r e p=πr22r2=π4 Em si, isso pode não parecer particularmente interessante. Porém, expressando π em função dessa probabilidade, seus estimadores podem ser usados ​​para aproximar π. Proceder analiticamente não é possível: a obtenção da probabilidade requer o conhecimento de π. Intuitivamente, é possível estimar essa probabilidade contando a proporção de gotas de chuva que estão dentro do círculo: se n gotas de chuva são observadas e m delas estão dentro do círculo, então pode-se estimar p, usando pˆ=m/n. A Figura 1 mostra uma simulação computacional em que 500 gotas de chuva foram distribuídas uniformemente sobre um quadrado – neste caso, pˆ=383/500 e, definindo πˆ pela relação entre p e π, πˆ=383/125=3,06. Uma estimativa ruim de π, considerando o esforço computacional usado, mas ainda assim é uma estimativa. De fato, m é uma realização de uma variável aleatória binomial (n, p). Objetivos de aprendizagem ------------------------- Até aqui, procuramos desenvolver o pensamento computacional aplicado à resolução de problemas. Neste capítulo vamos introduzir o Método de Monte Carlo para realizar simulações computacionais. A simulação computacional pode ser um passo importante no entendimento e modelagem de um problema, antes de chegarmos a uma solução. Procuraremos entender o processo da simulação para adquirirmos a habilidade para aplicarmos a técnica para possivelmente estimar valores relativos a fenômenos mais complexos para os quais não são conhecidas soluções. Simularemos alguns fenômenos simples para os quais já conhecemos uma solução, mesmo sendo capazes de chegar a uma solução analiticamente. Nesses casos a simulação pode ser entendida com uma ferramenta importante para comprovar a validade das soluções. A execução de experimentos aleatórios nas simulações é semelhante à maneira que funcionam jogos de azar, como dados, roletas, ... `Monte Carlo `_ é uma cidade mundialmente famosa por seus cassinos com seus jogos de azar. Vem daí a razão do nome de batismo do método. Método de Monte Carlo --------------------- O método de **Monte Carlo**, ou experimentos de Monte Carlo, são uma ampla classe de algoritmos computacionais que dependem de amostragens aleatórias repetidas para obter resultados numéricos. Usaremos esse método para estimarmos alguns valores. Da página do curso de `Introdução à Computação `_ de Princeton: Em 1953, Enrico Fermi, John Pasta e Stanslaw Ulam criaram o primeiro "experimento computacional" para estudar uma estrutura atômica vibratória. O sistema não linear não podia ser analisado pela matemática clássica. A **simulação de Monte Carlo** usa valores gerados aleatoriamente para variáveis incertas. O método se apoia na **lei dos grandes números**. Essa lei diz que a média aritmética dos resultados da realização da mesma experiência repetidas vezes tende a se aproximar do valor esperado à medida que mais tentativas se sucederem. Quanto mais tentativas são realizadas, mais a probabilidade da média aritmética dos resultados observados irá se aproximar da probabilidade real. O método é usado quando outras técnicas falham. Normalmente, a precisão é proporcional à raiz quadrada do número de repetições. Tais técnicas são amplamente aplicadas em vários domínios, incluindo: * projeto de reatores nucleares, * prever a evolução das estrelas, * prever o mercado de ações etc. Estimando a área de um círculo ------------------------------ Talvez você esteja pensando "mas para que estimar a área de um círculo se eu já sei calcular a área de um círculo"? Sim, é verdade, sabemos que a área de um círculo de raio :math:`\mathtt{r}` é :math:`\mathtt{\pi r^2}`. Mas, por hora, imagine que você não conhece a fórmula para calcular a área de um círculo. O que você faria para estimar essa área? Nesse primeiro exemplo, vamos ver como aplicar o `Método Científico `_, junto com Monte Carlo, para chegar a uma estimativa da área de um círculo e quão boa ou ruim é essa estimativa. O método científico está presente em tudo que fazemos cotidianamente. A ideia consiste em * **observar** algum aspecto da natureza * **hipotetizar** um modelo consistente com as observações * **prever** eventos usando as hipóteses * **verificar** as predições fazendo mais observações * **validar** repetindo até que hipóteses e observações estejam de acordo Nosso experimento de Monte Carlo considera a chuva que cai uniformemente ao acaso. A localização de qualquer gota de chuva será interpretada como uma variável aleatória uniformemente distribuída sobre uma região quadrada do espaço com um quadrante do círculo inscrito nesse quadrado, como na imagem a seguir. .. image:: ./imgs/area.* :width: 403 :height: 361 :scale: 100 :alt: Quadrante de um círculo :align: center Sem referência a qualquer teoria de probabilidade, é intuitivo que a probabilidade de uma gota de chuva cair em qualquer região dentro do quadrado deve ser proporcional à área dessa região e independe de sua localização. Consequentemente, a probabilidade :math:`\mathtt{p}` de que uma gota de chuva esteja dentro do círculo inscrito na imagem acima pode ser expressa como uma relação das áreas do quadrado e do círculo. No nosso experimentos sortearemos coordenadas ``(x,y)`` de pontos que correspondem a posições no quadrante onde uma gota de chuva caiu. Esse sorteio e a quantidade de gotas que caíram na região do círculo pode ser calculada fazendo .. code:: python x = random.random()*r y = random.random()*r if x*x + y*y < r*r: cont += 1 No caso, ``(x,y)`` representa a coordenada do ponto e ``cont`` será a variável responsável por contabilizar o número total de gotas que caíram no círculo. Uma breve descrição de algumas funções do módulo ``random`` do Python está na seção `Aleatoriedade em Python <./23-monte-carlo.html#aleatoriedade-em-python>`_ que está mais adiante. Cada **experimento** consistirá no sorteio de :math:`\mathtt{n}` pontos. A área do círculo será aproximada pela fração :math:`\mathtt{cont/n}` de pontos sorteados que caíram no interior do círculo. Sabemos que um quadrado de lados de comprimento :math:`\mathtt{r}` tem área :math:`\mathtt{r^2}`. Assim, se sortearmos :math:`\mathtt{n}` pontos, teremos que uma estimativa para a área do círculo pode ser computada fazendo-se :math:`\mathtt{r^2 \times (cont/n)}`. .. image:: ./imgs/MonteCarloIntegrationCircle.svg.* :width: 440 :height: 427 :scale: 50 :alt: By Yoderj, Mysid - Own work, CC0 :align: center :target: https://commons.wikimedia.org/w/index.php?curid=35608043 A classe ``Area`` a seguir será a implementação do laboratório do nosso experimento. Na classe * ``n`` é o número de pontos a serem sorteados; * ``t`` é o número de experimentos * ``p`` é a média aritmética das áreas obtidas pelos ``t`` experimentos * ``mean()`` é o método que retorna a aproximação da área do circulo * ``experimento()`` é o método responsável pela realização de cada experimento. .. code:: Python import random # número pré-determinado de experimentos T = 4096 class Area: #------------------------------------------ def __init__(self, n, t = T, r = 1): self.n = n # número de pontos usados em cada experimento self.t = t # número de experimentos self.raio = r # raio do círculo area_semicirculo = 0 for i in range(t): area_semicirculo += self.experimento() self.p = 4*area_semicirculo/t #------------------------------------------ def mean(self): '''(Area) -> float RETORNA a estimativa experimental da área de um círculo de raio self.r ''' return self.p #----------------------------------------- def experimento(self): '''(Area) -> float SIMULA um experimento sorteando self.n pontos no quadrante de lado self.r . RETORNA um estimativa da área do circulo de raio self.r ''' n = self.n r = self.raio if n == 0: return 0 cont = 0 for i in range(n): x = random.random()*r y = random.random()*r if x*x + y*y < r*r: cont += 1 return r*r*(cont/n) # probabilidade de cair no semicírculo. Realizando alguns experimentos com a classe ``Area`` chegamos aos valores mostrados na tabela a seguir. Todos os experimentos foram realizados com :math:`\mathtt{n=100}` e :math:`\mathtt{t=4096}`. .. math:: \begin{array}{clnt} \texttt{raio} & \texttt{área experimental} \\ \hline 1 & 3.141718750000016 \\ 2 & 12.575898437500069 \\ 4 & 50.2770312500003 \\ 8 & 200.752500000001 \\ 16 & 804.2175000000044 \end{array} Baseando-nos nesses resultados experimentais seria razoável propormos o seguinte. **Hipótese.** A área do círculo tem a forma :math:`\mathtt{a r^b}` para constantes :math:`\mathtt{a}` e :math:`\mathtt{b}`. Fazendo o gráfico de :math:`\mathtt{\log \mbox{área}}` por :math:`\mathtt{\log r}`, o chamado bilog, se obtivermos uma reta isso indicará que, de fato, a área é uma função da forma :math:`\mathtt{a r^b}` e que o coeficiente angular da reta é o valor aproximado para :math:`\mathtt{b}`. No nosso caso, obtivemos uma reta de coeficiente angular aproximadamente 2. .. image:: ./imgs/grafico_area.* :width: 640 :height: 480 :scale: 70 :alt: Gráfico da área do círculo :align: center A medida que o raio dobra a área do círculo é multiplicada por 4. Esse fato por si só sugere que o valor de :math:`\mathtt{b}` pode ser algo como 2. Tendo o valor do coeficiente :math:`\mathtt{b}` podemos estimar o valor da constante :math:`\mathtt{a}` usando a nossa hipótese e os resultados experimentais. Se olharmos na tabela o valor correspondente a :math:`\mathtt{r = 8}` obtemos que .. math:: \begin{align*} \mathtt{a} & = \mathtt{\text{área} / r^2} \\ & \mathtt{= 804.2175000000044 / 16^2} \\ & \mathtt{= 3.141474609375017} \end{align*} Como sabemos, esse valor para :math:`\mathtt{a}` é na verdade uma aproximação de :math:`\mathtt{\pi}`. O erro relativo da nossa estimativa obtida através do método de Monte Carlo foi portanto de aproximadamente: .. math:: \mathtt{\frac{|3.141474609375017-3.141592653589793\ldots|}{3.141592653589793\ldots}} \approx .0000375\ldots Nada mal, certo?! |:+1_tone5:| |:+1_tone1:| |:+1_tone3:| |:+1_tone2:| Isso indica que o método tem futuro. |:slight_smile:| Paradoxo do Aniversário ----------------------- A questão que trataremos agora, que é conhecida como `Paradoxo do Aniversário `_, é a seguinte. Em uma sala com :math:`\mathtt{n}` pessoas, qual é a probabilidade de termos duas ou mais que fazem aniversário em uma mesma data? Qual a probabilidade desse evento ocorrer quando :math:`\mathtt{n} = 20`? E quando :math:`\mathtt{n} = 40`. Essa probabilidade deve depender do número :math:`\mathtt{n}` de pessoas, por isso representaremos essa probabilidade como uma função :math:`\mathtt{p(n)}`. É evidente que, se :math:`\mathtt{n}` é um número maior ou igual ao número de dias em um ano, teremos que :math:`\mathtt{p(n)=1}`. |:sleeping:| Vejamos como o método de Monte Carlo pode nos ajudar a estimar o valor de :math:`\mathtt{p(n)}`. Podemos facilmente realizar um experimento com um grupo de amigas e amigos que se encontram para uma festa. |:moon_cake:| |:guitar:| |:tropical_drink:| |:champagne_glass:| |:cupcake:| |:champagne:| |:cup_with_straw:| À medida que o pessoal vai chegando podemos perguntar a data de aniversário de cada uma e de cada um e interromper o experimento, não a festa, assim que descobrirmos duas pessoas que fazem aniversário em uma mesma data. Será que precisaremos perguntar a data de aniversário a muitas pessoas para encontrar duas que fazem aniversário em uma mesma data? Simularemos esse experimento computacionalmente, já que é mais fácil e rápido e não atrapalha a festa ou encontro algum. Por simplicidade, na nossa modelagem, ignoraremos pequenas variações na distribuição das datas, tais como anos bissextos, gêmeos, variações sazonais ou semanais. Suporemos que há 365 datas de aniversários possíveis e que são todas igualmente prováveis. Representaremos cada data por um inteiro no intervalo ``range(365)``. Um **experimento** será a simulação do processo de :math:`\mathtt{n}` pessoas entrarem na sala. O **resultado do experimento** será ``SUCESSO`` se houver ao menos duas pessoas na sala com uma **mesma data** de aniversário na sala, e será ``FRACASSO`` em caso contrário. Para realizarmos esse experimento de forma computacional podemos representar a sala através de uma lista ou conjunto. Inicialmente essa lista está vazia. Para simular o processo de uma nova pessoa entrar na sala: * sorteamos a data do seu aniversário; * verificamos se a mesma data já foi sorteada antes e, se isso ocorreu, paramos com ``SUCESSO``; * se a data sorteada ainda não ocorreu, acrescentamos a data à lista ou conjunto. Suponha que ``DATAS`` é um apelido para o valor ``365``. Para sortearmos uma data basta escrevermos .. code:: python data = random.randrange(DATAS) O método de Monte Carlo consiste em executar esse experimento computacional um dado número :math:`\mathtt{t}` vezes para deixar a lei dos grandes números fazer a sua mágica. A estimativa da probabilidade :math:`\mathtt{p(n)}` será então o número de experimentos que terminaram em ``SUCESSO`` dividido por :math:`\mathtt{t}`. Por exemplo, imagine que desejamos estimar :math:`\mathtt{p(20)}`. Podemos executar :math:`\mathtt{t = 1000}` experimentos e contar o número :math:`\mathtt{nS}` de ``SUCESSO``. A estimativa para :math:`\mathtt{p(20)}` será a razão :math:`\mathtt{nS/t}`. Agora, como exercício, complete a implementação da classe ``Aniversario`` escrevendo o método ``experimento()``. Um possível função ``main()`` *cliente* dessa classe é a seguinte. A função ``prob_aniversario()`` é descrita e apresentada na próxima seção. .. code:: python def main(): ''' Unidade de teste da classe Aniversario ''' print("Paradoxo do Aniversário") print("probabilidade calculada") p = prob_aniversario(23) print(f"p({23}) = {p}") print("probabilidade estimada") a = Aniversario(23, 1000) print(f"{a}") print(f"erro relativo = {abs(a.mean() - p)/p:.2f}") Descubra quão longe ou quão perto as probabilidades obtidas pelo método de Monte Carlo estão das probabilidade reais. .. raw:: html Aniversário sem simulação ------------------------- É possível determinarmos :math:`\mathtt{p(n)}` calculando a probabilidade :math:`\mathtt{q(n)}` de que todos os :math:`\mathtt{n}` aniversários sejam diferentes. Tendo :math:`\mathtt{q(n)}` teremos que :math:`\mathtt{p(n) = 1 -q(n)}` Se :math:`\mathtt{n > 365}`, pelo *Princípio da Casa dos Pombos*, temos que :math:`\mathtt{q(n) = 0}`. Por outro lado, se :math:`\mathtt{n \leq 365}`, ela é dada por .. math:: \begin{align*} \mathtt{q(n)} & \mathtt{= 1 \times (1 - 1/365 ) \times (1 - 2/365 ) \cdots (1-(n-1)/365)} \\ & \mathtt{= \frac{365 \times 364 \cdots (365-n+1)}{365^n}}\\ & \mathtt{= \frac{365!}{365^n(365-n)!}}. \end{align*} A seguinte função ``prob_aniversario(n)`` que se baseia nesses cálculos pode ser utilizada para obtermos alguns valores de :math:`\mathtt{p(n)}`. .. code:: python def prob_aniversario(n): '''(int) -> float RECEBE um inteiro n. RETORNA a probabilidade de em um grupo de n pessoas haja pelo menos 2 que fazem aniversário em um mesmo dia. Implementação alternativa: import math num = math.factorial(DATAS) den = math.factorial(DATAS-n) prob_c = num/(DATAS**n * den) ''' prob_c = 1 for i in range(n): prob_c *= (1 - i/DATAS) return 1 - prob_c Admiremos alguns valores de :math:`\mathtt{p(n) = 1 - q(n)}`: .. math:: \begin{array}{r|r} \mathtt{n} & \mathtt{p(n)} \\ \hline 8 & 0.07 \\ 10 & 0.11 \\ 12 & 0.16 \\ 16 & 0.28 \\ 18 & 0.34 \\ 20 & 0.41 \\ 22 & 0.47 \\ 24 & 0.53 \\ 30 & 0.70 \\ 40 & 0.89 \end{array} Um fato curioso mostrado nessa tabela é o de que funções injetoras são raras. Se considerarmos todas as funções de ``range(24)`` a ``range(365)`` os cálculos que acabamos de fazer mostram que mais da metade delas **não são injetoras**. Podemos usar computação para explorar esse problema e obtermos uma aproximação aceitável para :math:`\mathtt{p(n)}`? Aqui é que entra o Método de Monte Carlo. O interessante de usarmos o método em um problema que obtivemos a resposta analiticamente é verificar quão boas ou ruins são as aproximações que obtemos com o método. Colecionador de Figurinhas -------------------------- Trataremos agora de um problema bem atual, o problema do(a) colecionador(a) de figurinhas (`Coupon collector's problem `_). Suponha que temos um álbum com :math:`\mathtt{n}` figurinhas que são numeradas de 0 a :math:`\mathtt{n-1}`. Uma forma de completar o álbum é frequentar lugares onde aficionadas e aficionados por completar o álbum se encontram para trocar figurinhas repetidas por figurinhas que estão faltando no nosso álbum. Suponha agora que uma colecionadora curiosa deseja saber quantas figurinhas ela precisa comprar para completar seu álbum **sem trocar** figurinhas repetidas, apenas **comprando** figurinhas até que o álbum esteja completo. A questão que desejamos tratar é Quantas figurinhas a colecionadora deve comprar para completar o álbum? Suporemos que compramos uma figurinha por vez e que cada uma é **igualmente provável** de ser comprada. Assim, não consideraremos as figurinhas brilhantes como sendo *mais difíceis*. Escreveremos uma classe ``Colecionadora`` que utiliza o o método de Monte Carlo para estimar esse número de figurinhas. A classe ``Colecionadora`` recebe dois número inteiros :math:`\mathtt{n}` e :math:`\mathtt{t}` e realiza :math:`\mathtt{t}` experimentos aleatórios. Em cada experimento a ``Colecionadora`` conta o número de figurinhas que teve comprar para preecher um álbum de :math:`\mathtt{n}` figurinhas. Com esses números em mãos a ``Colecionadora`` estima o número de figurinhas que devem ser compradas para completar o álbum. A cada instante uma figurinha dentre as :math:`\mathtt{n}` possíveis é comprada com a mesma probabilidade. O adjetivo *aleatório* dos experimentos é devido a essa tal probabilidade. A compra de cada figurinha se assemelha ao lançamento de um dado *honesto* de :math:`\mathtt{n}` faces. Hãã?!. |:dizzy_face:| A face que fica para cima após o lançamento é o número da figurinha comprada. Por exemplo, suponha que tenhamos um álbum de :math:`\mathtt{n = 5}` figurinhas e que a ``Colecionadora`` deve realizar :math:`\mathtt{t = 3}` experimentos. A sequência de figurinhas compradas no primeiro experimento pode ser .. raw:: html 
   4   0   1   1   1   1   4   3   0   1   1   0   2
em que as novas figurinhas adquiridas estão em *vermelho*. Nesse primeiro experimento a ``Colecionadora`` teve que comprar :math:`\mathtt{f_0 = 13}` figurinhas até completar o álbum de :math:`\mathtt{n = 5}` figurinhas. As figurinhas têm número 0, 1, 2, 3, e 4, certo?! No segundo experimento as figurinhas compradas poderiam ter sido .. raw:: html 
   4   3   0   2   1
Uau! |:fireworks:| A ``Colecionadora`` comprou exatamente :math:`\mathtt{f_1 = 5}` figurinhas para completar o álbum. Finalmente, já que :math:`\mathtt{t = 3}`, a ``Colecionadora`` precisa realizar apenas um último experimento. Suponha que desta vez a ordem das figurinhas compradas seja .. raw:: html 
   3   3   3   3   4   3   2   0   4   0   1
Desta vez foram compradas :math:`\mathtt{f_2 = 11}` figurinhas para completar o álbum. Antes de terminar o seu serviço e poder ir para casa a ``Colecionadora`` calcula o número médio de figurinhas compradas nos :math:`\mathtt{t = 3}` experimentos. Esse número é :math:`\mathtt{F = (f_0 + f_1 + f_2)/3 = (13 + 5 + 11)/3 = 9{,}66...}`. A ``Colecionadora`` declara que :math:`\mathtt{9.66...}` é o número estimado de figurinhas que devem ser compradas para preencher um álbum de 5 figurinhas, fecha o laboratório e vai para casa descansar. Não se aborreça com o fato de :math:`\mathtt{9.666}` não ser um número inteiro. Ele é apenas um estimador fazendo o seu trabalho de estimar. Um exemplo de cliente da classe ``Colecionadora`` é a seguinte função ``main()``. A função ``no_figurinhas()`` é descrita e apresentada na próxima seção. .. code:: python def main(): ''' Unidade de teste da classe Colecionadora ''' print("Colecionadora de Figurinhas") print("número de figurinhas calculado") n = no_figurinhas(400) print(f"f({400}) = {n}") print("número de figurinhas estimado") colecionadora = Colecionadora(400, 100) print(f"{colecionadora}") print(f"erro relativo = {abs(colecionadora.mean() - n)/n:.2f}") Na classe ``Colecionadora`` que está a seguir ficou como exercício implementar o método ``experimento()``. .. raw:: html Figurinhas sem simulação ------------------------ É possível determinarmos diretamente o número esperado de figurinhas que a colecionadora deve comprar para completar o album. Seja :math:`\mathtt{T}` o número de figurinhas a serem compradas para completarmos um álbum. Seja :math:`\mathtt{t_i}` o número de figurinhas que compramos entre obtermos a :math:`\mathtt{(i-1)}`-ésima e a :math:`\mathtt{i}`-ésima figurinha distinta. :math:`\mathtt{T}` e :math:`\mathtt{t_i}` são variáveis aleatórias e :math:`\mathtt{T = t_1 + t_2 + \cdots + t_n}`. A probabilidade :math:`\mathtt{p_i}` de conseguirmos uma figurinha distinta, dado que já conseguimos :math:`\mathtt{i-1}` figurinhas distintas é :math:`\mathtt{(n - (i-1))/n}`. Seja :math:`\mathtt{q_i = 1 - p_i = 1/(i-1)}`. Assim, o valor esperado de :math:`\mathtt{t_i}` é dado por .. math:: \begin{align*} \mathtt{E(t_i)} & \mathtt{= 1 p_i + 2 p_i q_i + 3 p_i q_i^2 + \cdots} \\ & \mathtt{= p_i (1 + 2 q_i + 3 q_i^2 + 4 q_3^3 + \cdots)} \\ & \mathtt{= p_i 1/(1-q_i)^2 = 1/p_i.} \end{align*} Desta forma, temos que .. math:: \begin{align} \mathtt{E(T)} & \mathtt{= E(t_1) + E(t_2) + \cdots + E(t_n)} \\ & \mathtt{= n/n + n/(n-1) + n/(n-2)+ \cdots n/1} \\ & \mathtt{= n \times H_n} \\ & \mathtt{< n \ln (n+1)} \end{align} em que :math:`\mathtt{H_n}` é o :math:`\mathtt{n}`-ésimo `número harmônico `_. Portanto, temos que comprar aproximadamente :math:`\mathtt{n \ln n}` figurinhas para completarmos um album de :math:`\mathtt{n}`, supondo que a probabilidade ... Para um album de 400 figurinhas é esperado que tenhamos que comprar cerca de 2400 figurinhas. A seguinte função ``no_figurinhas()`` calcula e retorna o número esperado de figurinhas a serem compradas. .. code:: python def no_figurinhas(n): '''(int) -> float RECEBE um inteiro. RETORNA o número esperado de figurinhas que devem ser compradas para completarmos um álbum de n figurinhas. ''' hn = 0 for i in range(n,0,-1): hn += 1/i return n*hn Urnas eletrônicas ----------------- Suponha que em uma população de 100 milhões de eleitoras e eleitores, 51% votaram para a(o) candidata(o) ``A`` e 49% votaram para a(o) candidata(o) ``B``. Na votação foi utilizada uma urna eletrônica que é propensa a contabilizar 5% dos votos de maneira errada. Quando erra a urna troca a(o) candidata(o) que recebeu o voto. Se um voto foi para ``A``, a urna contabiliza para ``B`` e se o voto foi para ``B``, ela contabiliza para ``A``. **Problema.** Supondo que os erros ocorrem **uniformemente ao acaso**, 5% é uma porcentagem de erro suficiente para que seja provável, digamos que com mais de 50%, 25%, 10% ou 5% de chance, que o resultado da eleição tenha sido alterado? Qual porcentagem de erro poderia ser tolerada para que, digamos com mais de 95% de certeza, não haja alteração no resultado? Há certamente várias questões aqui que despertam paixões e que são associadas a problemas similares, afinal, 1% de incerteza ainda é incerteza. Nos concentraremos em aspectos que podem ser investigados através do Método de Monte Carlo. Está seção é baseada no exercício 5 do `Capítulo 2 `_ do livro `Computer Science: An Interdisciplinary Approach `_ . Modelagem é a procura por uma representação que aproxime de maneira *satisfatória* um fenômeno real que é difícil de observar diretamente. Modelos são usados para explicar e prever o comportamento de sistemas reais e são aplicados em uma variedade de atividades do dia a dia. Embora a modelagem seja um componente central da ciência, os modelos, na melhor das hipóteses, são **aproximações** dos sistemas que representam, não são réplicas exatas. Todas e todos estão constantemente trabalhando para melhorar e tornar os modelos mais satisfatórios e úteis. Na nossa modelagem da Eleição um componente central será a **aleatoriedade**: suporemos que os erros na contabilização dos votos ocorrerão *uniformemente ao acaso*. Isso significa que cada voto poderá ser contabilizado erroneamente com a mesma probabilidade. Por exemplo, suponha que nossa eleição tenha 20 votos e esses votos sejam .. raw:: html 
   A   B   B   A   B   A   A   B   A   B   A   A   B   A   B   B   A   A   A   B
Nessa eleição ``A`` venceu já que obteve 11 votos ante 9 votos recebidos por ``B``. Se as urnas eletrônicas 20% das vezes contabilizam erroneamente um voto, então nessa eleição de 20 votos teremos :math:`\mathtt{20 \times 0.2 = 4}` votos que serão contabilizados erroneamente. Assim, uma contabilização possível dos votos com erros poderia ser .. raw:: html 
   A   B   B   B   B   A   A   B   A   A   A   A   B   A   A   B   A   A   A   A
em que os votos em *vermelho* são aqueles contabilizados de forma errada. Como podemos ver, nesse caso o resultado da eleição *não foi alterado*, foram contabilizados 13 votos para ``A`` e 7 votos para ``B``. Uma outra contabilização possível poderia ser .. raw:: html 
   A   B   B   A   B   B   A   B   A   B   B   A   B   B   B   B   A   B   A   A
e nessa o resultado da eleição *foi alterado*, ``B`` passou a vencer por 12 votos a 8. A fim de obter uma estimativa para a probabilidade do resultado da eleição ser alterado escreveremos uma classe ``Eleicao`` que realiza experimentos de Monte Carlo. A estrutura da classe ``Eleicao`` é semelhante a das classes ``Area``, ``Aniversario`` e ``Colecionadora`` que fizemos. O construtor da classe **recebe** como argumentos: * ``no_votos``: um inteiro com o número total de votos em um da(o)s dois candidata(o)s; * ``votos_A`` : um inteiro que é número de votos recebidos pela(o) condidada(o) ``A``; * ``erro`` : um inteiro com porcentagem de votos contabilizados erroneamente; * ``t`` : um inteiro que indica o número de experimentos que devem ser realizados A classe supõe que a(o) candidata(o) ``A`` foi a(o) vencedora(or) da eleição, ``votos_A > no_votos/2``. Para estimar a probabilidade, um objeto da classe ``Eleicao`` realizará ``t`` experimentos. Se ``k`` é o número de experimentos em que o resultado da eleição foi alterado, a(o) candidata(o) ``B`` passou a ser a(o) vencedor(a) ou chegamos a um empate, então o objeto deve retornar através de um método ``prob()`` a razão :math:`\mathtt{k/t}`. Abaixo está um exemplo de *cliente* da classe ``Eleicao``: .. sourcecode:: python # constantes para a eleição NO_VOTOS = 100000000 # número total de votos VOTOS_A = 51000000 # número de votos em A VOTOS_B = NO_VOTOS - VOTOS_A # número de votos em B ERRO = 5 # porcentagem de erro T = 10 # número de experimentos ## ================================================================== # def main(): ''' Testes das classes Eleicao ''' # teste eleição print("Experimentos da classe Eleicao") eleicao = Eleicao(NO_VOTOS, VOTOS_A, 2, T) # 2% de erro print(f"Eleicao(v={NO_VOTOS},a={VOTOS_A},e={2},t={T})={eleicao.prob()}") print(f"{Eleicao(no_votos=NO_VOTOS, votos_A=VOTOS_A, erro=ERRO, t=T)}") # 5% de erro print(f"{Eleicao(no_votos=NO_VOTOS, votos_A=VOTOS_A, erro=10, t=T)}") # 10% de erro print(f"{Eleicao(no_votos=NO_VOTOS, votos_A=VOTOS_A, erro=20, t=T)}") # 20% de erro print(f"{Eleicao(no_votos=NO_VOTOS, votos_A=VOTOS_A, erro=30, t=T)}") # 30% de erro A saída esperada produzida por essa função ``main()`` terá a forma: :: Experimentos da classe Eleicao Eleicao(v=100000000,a=51000000,e=2,t=10)=?.? Eleicao(v=100000000,a=51000000,e=5,t=10)=?.? Eleicao(v=100000000,a=51000000,e=10,t=10)=?.? Eleicao(v=100000000,a=51000000,e=20,t=10)=?.? Eleicao(v=100000000,a=51000000,e=30,t=10)=?.? Em que as ``?.?`` representam as probabilidades estimadas pelos experimentos. A sua implementação da classe ``Eleicao`` será deixada como exercício. Depois de realizar os experimentos, procure verificar se a estimativa da probabilidade do resultado da eleição ser alterado parece razoável. Procure produzir uma argumentação que justifique o comportamento observado. Se não for possível produzir uma tal argumentação talvez o experimento esteja errado. De qualquer forma, fazer mais experimentos pode dar uma luz sobre o fenômeno e sobre o que esperar. É desejável que os experimentos possam ser realizados várias vezes. O número de vezes que um experimento é realizado está por conta do parâmetro ``t`` dos objetos da classe ``Eleição``. Se um objeto da classe ``Eleicao`` estiver fazendo o seu serviço muito lentamente há algumas coisas que podem ser tentadas para acelerá-lo. Como fizemos em `14. Operações com imagens `_ podemos apelar para NumPy. Frequentemente os comandos de repetições são chamados simplemente de **laços** ou **loops**. Objetos ``numpy.ndarrays`` permitem que troquemos o uso de laços que manipulam seus elementos um a um chamadas operações *vetoriais*. Essa técnica é conhecida pelo nome de **vetorização**. Seguem alguns exemplos em que laços como ``while ...``, ``for ...`` são substituídos por operações vetoriais que são muito mais eficiente. No que segue ``ar1`` é um objeto ``numpy.ndarray``, um array. Por isso quando escrevemos ``ar1 > 4`` teremos um novo ``numpy.ndarray`` de ``dtype=bool`` e com o mesmo ``shape`` de ``ar1``. Nesse novo ``numpy.ndarray`` temos que cada posição de índice ``i`` é ``True`` se ``ar1[i] > 4`` é ``False`` se ``ar1[i] <= 4``: .. sourcecode:: python In [4]: import numpy as np In [5]: lst = [ 1, 6, 7, 2, 6, 5, 3] In [6]: ar1 = np.array(lst) In [7]: ar1 Out[7]: array([1, 6, 7, 2, 6, 5, 3]) In [8]: ar1 > 4 # isto um exemplo de vetorização Out[8]: array([False, True, True, False, True, True, False]) Se desejarmos saber o número de ``True`` no array ``ar1 > 4`` podemos fazer .. sourcecode:: python In [9]: ar2 = ar1 > 4 In [10]: np.sum(ar2) # o resultado é o número de True em ar2 Out[10]: 4 Na adição feita pela função ``np.sum()`` tudo se passa como se ``True`` fosse 1 e ``False`` fosse 0. Para gerarmos um ``numpy.ndarray`` com uma amostra aleatória de 3 números inteiros em ``range(5)`` podemos usar a função ``np.random.choice()`` e escrever .. sourcecode:: python In [13]: np.random.choice(5, 3, replace=False) # gera uma amostra aleatória de np.arange(5) sem repetições Out[13]: array([2, 0, 4]) O argumento ``replace=False`` indica que a amostra será criada sem repetições. Sem esse argumento a amostra poderá ter elementos repetidos: .. sourcecode:: python In [15]: np.random.choice(5, 3) # gera uma amostra aleatória de range(5) com possíveis repetições Out[15]: array([4, 2, 2]) In [16]: np.random.choice(5, 3) # gera uma amostra aleatória de range(5) com possíveis repetições Out[16]: array([3, 1, 0]) Para simularmos ``6`` lançamentos de uma moeda em que a probabilidade de cara é ``p=0.3`` e de coroa ``0.7`` podemos utilizar a função ``np.random.binomial()`` e escrevemos: .. sourcecode:: python In [56]: np.random.binomial(n=1, p=0.3, size=6) Out[56]: array([1, 1, 1, 0, 0, 1]) Neste exemplo estamos interpretando ``1`` como cara e ``0`` como coroa. Para termos o número de caras ou de coroas na amostra podemos fazer: .. sourcecode:: python In [57]: lances = np.random.binomial(n=1, p=0.3, size=6) In [58]: lances Out[58]: array([0, 0, 1, 0, 0, 0]) In [59]: caras = np.sum(lances == 1) In [60]: coroas = np.sum(lances==0) In [61]: caras Out[61]: 1 In [62]: coroas Out[62]: 5 O número de caras em ``20000`` lançamentos de uma moeda em que a probabilidade de cara é ``0.1`` pode ser simulado com .. sourcecode:: python In [84]: np.sum(np.random.binomial(1, 0.1, 20000) == 1) Out[84]: 2057 Aleatoriedade em Python ----------------------- Em Python temos a classe ``random`` com vários métodos e funções para gerarmos objetos aleatoriamente. Veja a página `Generate pseudo-random numbers `_. Entre esses métodos e funções temos alguns que já usamos: * ``random.seed(a=None, version=2)``: inicializa o gerador de números aleatórios * ``random.randrange(start, stop[, step])``: retorna um valor escolhido uniformemente ao acaso em ``range(start, stop[,step])`` * ``random.choice(lst)``: retorna um item escolhido uniformemente ao acaso da lista ``lst``. * ``random.shuffle(lst[, random])``: embaralha *in place* os itens na lista ``lst``. * ``random.sample(populacao, k)``: retorna uma lista com ``k`` elementos da ``polulacao`` que é um objeto da classe ``list`` ou ``set`` .. code:: python In [1]: import random In [2]: random.seed("oi") In [3]: random.randrange(27) Out[3]: 7 In [4]: random.seed("oi") In [5]: random.randrange(27) Out[5]: 7 In [6]: random.seed("Eduardo") In [6]: random.randrange(27) Out[6]: 14 In [7]: lst = ["joão", "pedro", "juliana", "giovana", "julia"] In [8]: random.choice(lst) Out[8]: 'giovana' In [9]: random.choice(lst) Out[9]: 'julia' In [10]: random.choice(lst) Out[10]: 'julia' In [11]: random.choice(lst) Out[11]: 'pedro' In [12]: lst Out[12]: ['joão', 'pedro', 'juliana', 'giovana', 'julia'] In [13]: random.shuffle(lst) In [14]: lst Out[14]: ['julia', 'juliana', 'pedro', 'giovana', 'joão'] In [15]: random.shuffle(lst) In [16]: lst Out[16]: ['giovana', 'joão', 'pedro', 'juliana', 'julia'] In [17]: random.sample(lst,2) Out[17]: ['giovana', 'julia'] In [18]: random.sample(lst,2) Out[18]: ['juliana', 'pedro'] In [19]: random.sample(lst,2) Out[19]: ['joão', 'julia'] Eficiência de ``set`` --------------------- Python possui o tipo nativo ``set`` que permite a criação e manipulação de conjuntos. Um conjunto é uma coleção não ordenada sem elementos duplicados. Os usos básicos incluem testes de adesão e eliminando entradas duplicadas. Objetos da classe conjunto admitem operações como união, intersecção, diferença e diferença simétrica. Internamente, para o Python, conjuntos são dicionários em que o valor associado a cada chave é irrelevante, apenas as chaves, elementos do conjunto, importam. Por essa razão, de maneira semelhante ao que ocorre com dicionário, os elementos de um conjunto em Python devem ser objetos imutáveis como números, strings e tuples. Para se criar um conjunto vazio podemos escrever .. code:: python In [2]: s = {} Chaves ou a função ``set()`` podem ser usadas para criar conjuntos. Para criar um conjunto vazio devemos usar ``set()`` e não ``{}``, já que o par de chaves ``{}`` é usado para criar um dicionário vazio. .. code:: python In [1]: s = set('abracadabra') In [2]: t = set('alacazam') In [3]: s Out[3]: {'a', 'r', 'b', 'c', 'd'} In [4]: s - t # letras em s mas não em t Out[4]: {'r', 'd', 'b'} In [5]: s | t # letras em s ou em t ou em ambos Out[5]: {'a', 'c', 'r', 'd', 'b', 'm', 'z', 'l'} In [6]: s & t # letras em ambos s e t Out[6]: {'a', 'c'} In [7]: s ^ t # letras em s ou em t Out[7]: {'r', 'd', 'b', 'm', 'z', 'l'} O comportamento básico de objeto do tipo ``set`` também ilustrado pelo trecho de código abaixo. .. code:: python def main(): # crie um conjunto vazio conj0 = set() print(f' 1: conj0: {conj0} tem tamanho {len(conj0)}') # crie um conjunto e obeserve que não há ordem entre os elementos de um conjunto conj1 = { 'set', "e'", (2,2), 'da hora', '!!' } print(f' 2: conj1: {conj1} tem tamanho {len(conj1)}') # use o método add() para inserir elementos em um conjunto for i in range(5): conj0.add( (i, 2) ) print(' 3: conj0:', conj0) # use o método union() para unir conjuntos print(' 4: uniao:', conj0.union(conj1)) # use o método intersection() para determinar a intersecção de conjuntos print(' 5: intersecção:', conj1.intersection(conj0)) # uso de 'for ... in ...' para percorrer os elementos de um conjunto. for elemento in conj0: # testa se elemento é um elemento de conj1 if elemento in conj1: print(f" >>> {elemento} pertence a {conj1}") else: print(f" x {elemento} não pertence a {conj1}") A função ``main()`` produz a saída a seguir. Examinar a função ``main()`` e o seu comporatmento pode ser instrutivo para compreender o funcionamento de objetos da classe ``set``. :: 1: conj0: set() tem tamanho 0 2: conj1: {'!!', 'set', (2, 2), 'da hora', "e'"} tem tamanho 5 3: conj0: {(1, 2), (4, 2), (0, 2), (2, 2), (3, 2)} 4: uniao: {(1, 2), (4, 2), (0, 2), '!!', 'set', (2, 2), 'da hora', (3, 2), "e'"} 5: intersecção: {(2, 2)} x (1, 2) não pertence a {'!!', 'set', (2, 2), 'da hora', "e'"} x (4, 2) não pertence a {'!!', 'set', (2, 2), 'da hora', "e'"} x (0, 2) não pertence a {'!!', 'set', (2, 2), 'da hora', "e'"} >>> (2, 2) pertence a {'!!', 'set', (2, 2), 'da hora', "e'"} x (3, 2) não pertence a {'!!', 'set', (2, 2), 'da hora', "e'"} .. raw:: html Agora vamos aos consumos de tempo. Suponha que o dicionário ``s`` e ``t`` são conjuntos ``x`` é um objeto imutável qualquer. A função ``len()`` quando usada com um objeto da classe ``set`` retorna o número de elementos no conjunto. A tabela foi copiada de `TimeComplexity `_. .. math:: \begin{array}{lc} \text{operação} & \text{consumo de tempo médio} \\ \hline \texttt{len(s)}& \mathtt{\textup{O}(1)}\\ \texttt{x in s}& \mathtt{\textup{O}(1)}\\ \text{união: }\texttt{s | t} & \mathtt{\textup{O}(len(s)+len(t))}\\ \text{intersecção: }\texttt{s & t} & \mathtt{\textup{O}(min(len(s), lent(t))}\\ \text{diferença : }\texttt{s - t}& \mathtt{\textup{O}(len(s))}\\ \text{diferença simétrica: }\texttt{s^t}& \mathtt{\textup{O}(len(s))} \end{array} Lições ------ .. image:: https://online.anyflip.com/bybf/cgdb/files/mobile/1.jpg?1575726812 :width: 500 :height: 700 :scale: 10 :alt: This work is licensed under a Creative Commons Attribution-NonCommercial 2.5 License. :align: left :target: https://online.anyflip.com/bybf/cgdb/files/mobile/1.jpg?1575726812 Uma das lições que aprendemos com nesses experimentos é que vale a pena trocar figurinhas repetidas para completar álbuns. Vale a pena significa que apenas comprar figurinhas para completar um album sai muito caro. Se bem que isso já era intuitivo. Não era necessário o Método de Monte Carlo para chegarmos a essa conclusão. O ponto é que, frequentemente, devido à complexidade envolvida na obtenção de resultados exatos através de modelos analíticos, as simulações computacionais do Método de Monte Carlo são ferramentas muito úteis para estimarmos algum valor de interesse. A simulação deve ser minimamente eficiente para que obtenhamos estimadores sem termos que esperar até o final dos tempos. Para tornar as simulações efetivas entram em campo várias atrizes e atores como biologia, computação, estatística, engenharia, física, matemática,... todas e todos trabalhando lado a lado. Os programas que consideramos ao longo do ano são um ponto de partida, não um estudo completo. Esta disciplina também é um ponto de partida para o seus estudos posteriores em ciências ou engenharias. A abordagem de programação (*computational thinking*) e as ferramentas que você aprendeu devem auxiliá-la(o) a resolver vários problemas computacionais que serão encontrados no restante de seu curso e na vida profissional. Exercícios ---------- Para saber mais --------------- * `Conjuntos `_. * `Generate pseudo-random numbers `_. * `Mathematical statistics functions `_. * `Monte Carlo Simulation `_. * `Random sampling (numpy.random) `_. * `Resampling: The New Statistics `_.