Estatística Computacional com Python
Estatística Computacional com Python
2024-08-11
ii
Table of contents
Prefácio 1
1 Introdução ao Python 3
1.1 Introdução aos Comandos e Objetos do Python . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 Operações Aritméticas Básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Variáveis Booleanas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Strings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Listas, Tuplas, Conjuntos e Dicionários . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.1 Listas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2 Tuplas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4.3 Conjuntos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4.4 Dicionários . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.5 Matrizes e Arranjos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 Arquivos de Dados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.7 Estruturas de Controle de Programação . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.8 Funções . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.9 Estatística Computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
1.10 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
iii
iv TABLE OF CONTENTS
6 Aproximação de Distribuições 85
6.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2 Modelos Probabilísticos Discretos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3 Modelos Probabilísticos Contínuos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.4 Quadraturas Gaussianas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.5 Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
6.6 Funções Pré-Existentes no Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.7 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
References 145
Prefácio
1
2 Prefácio
Chapter 1
Introdução ao Python
O programa Python foi escolhido para ministrar este curso por uma série de razões. Além de ser um
programa livre, no sentido de possuir livre distribuição e código fonte aberto, pode ser utilizado nas
plataformas Windows e Linux. Além do mais, o Python possui grande versatilidade no sentido de possuir
inúmeros pacotes já prontos e nos possibilitar criar novas rotinas e funções. O PyPi é o repositório oficial
do Python onde todos os pacotes são armazenados. Você pode pensar nele como um Github para os
pacotes do Python. O Python foi criado pelo holandês Guido van Rossum para ser uma linguagem de
programação simples e legível, além de ser muito produtiva. O Python evoluiu e se tornou em uma
linguagem muito atrativa e uma das principais escolhas para aplicações de desenvolvimento web, análise
de dados e inteligência artificial, entre outras. Por ser genuinamente um programa orientado por objeto
nos possibilita programar com muita eficiência e versatilidade, embora apresente algumas mudanças em
sua implementação em relação a outras linguagens orientadas por objetos. Outro aspecto que é bastante
atrativo no Python refere-se ao fato de o mesmo receber contribuições de pesquisadores de todo o mundo
na forma de pacotes. Essa é uma característica que faz com que haja grande desenvolvimento do programa
em relativamente curtos espaços de tempo e que nos possibilita encontrar soluções para quase todos os
problemas com os quais nos deparamos em situações reais. Para os problemas que não conseguimos
encontrar soluções, o ambiente de programação Python nos possibilita criar nossas próprias soluções.
Nestas notas de aulas pretendemos apresentar os conceitos básicos da estatística computacional de uma
forma bastante simples. Inicialmente obteremos nossas próprias soluções para um determinado método ou
técnica e em um segundo momento mostraremos que podemos ter a mesma solução pronta do Python
quando esta estiver disponível. Particularmente neste capítulo vamos apresentar algumas características
do ambiente e da linguagem para implementarmos nossas soluções. Nosso curso não pretende dar soluções
avançadas e de eficiência máxima para os problemas que abordaremos, mas propiciar aos alunos um
primeiro contato com a linguagem Python e com os problemas básicos da estatística computacional.
A desvantagem é que o Python não é um programa fácil de aprender. Alguns esforços iniciais são necessários
até que consigamos obter algum benefício. Não temos a intenção de apresentar neste curso os recursos
do Python para análises de modelos lineares de posto completo ou incompleto, de modelos não-lineares,
de modelos lineares generalizados ou de gráficos. Eventualmente poderemos utilizar algumas destas
funções como um passo intermediário da solução do problema que estaremos focando. Este material será
construído com uma breve e simplificada abordagem teórica do tópico e associará exemplificações práticas
dos recursos de programação Python para resolver algum problema formulado, em casos particulares da
teoria estudada.
Este material é apenas uma primeira versão que deverá ter muitos defeitos. Assim, o leitor que encontrá-los
ou tiver uma melhor solução para o problema poderá contribuir enviando um e-mail para danielff@[Link].
Visite minha homepage [Link]
3
4 CHAPTER 1. INTRODUÇÃO AO PYTHON
Uma vez instalado, podemos digitar os códigos e com Ctrl Enter executamos os códigos linha por linha
ou um bloco de linhas marcadas. É importante instalarmos algumas bibliotecas básicas, caso elas já não
estejam instaladas. A seguir, temos um código Python para esse propósito.
pip install numpy
pip install sympy
pip install PIL
1.1. INTRODUÇÃO AOS COMANDOS E OBJETOS DO PYTHON 5
Em seguida devemos importar as libraries que precisarmos. No script a seguir consideramos a importação
de todas as libraries. A library math não precisa ser instalada, pois já vem com as distribuições do Python.
import math
import numpy
import sympy
import PIL
import jupyter
import matplotlib
6
59056
1.2625
Podemos observar que o símbolo # é usado para inserirmos comentários no código Python e o operador /
faz divisão usando operadores reais. Para divisão de inteiros, podemos usar //, assim, 6 // 5 retorna 1
e 6 / 5 retorna 1,2. O resto da divisão por inteiro é obtido pelo operador %. Assim, 6 % 5 retorna 1.
Temos também que o operador ** é a função potência, ou seja, por exemplo, 310 é 3**10 em Python.
O Python também pode realizar operações com números complexos, que no caso, √ são representados por
a + bj, em que a é a parte real do número e bj, a parte imaginária, sendo j = −1. O j é representado
por i nos livros de matemática e de outras áreas. O programa a seguir ilustra uma operação com números
complexos dada por (6 − 4i)2 . Assim, temos
(6-4j)**2
(20-48j)
Apresentamos a seguir um script que faz uso da library math, cujo
√ primeiro comando foi para importá-la,
o penúltimo para obter o valor de π e o último comando calculou 2. Também fizemos mais uma operação
com números complexos.
import math
2 + 4 + 5.6
2 / 3 - 4
(3-4j)*(3+4j)
[Link]
[Link](2)
11.6
6 CHAPTER 1. INTRODUÇÃO AO PYTHON
-3.3333333333333335
(25+0j)
3.141592653589793
1.4142135623730951
As bibliotecas numpy e sympy são para diversos cálculos matemáticos, sendo que a última efetua cálculos
simbólicos.
import sympy
import numpy
numpy.set_printoptions(legacy='1.25')
[Link]([Link]/5)
[Link]([Link]/5)
type(1.5 + 2.1j) # tipo do objeto
√
s
5 5
−
8 8
0.5877852522924731
complex
Podemos realizar uma operação matemática básica como algumas das anteriormente apresentadas ou até
mesmo, mais complexa e armazenar o valor em uma variável, digamos x. Essa variável pode ser usada
para outras operações matemáticas e até simbólicas, se usarmo o sympy. Veja Knuth [1984] para discussão
sobre programação simbólica. O script a seguir ilustra alguns casos deste procedimento.
x = [Link]('y')
z = (1+x**2)**2
[Link](z)
y = 2.3 + 6.7**2
r = y**2 + 1 / 2**3
print('r =', r, 'y = ', y, 'e', 'z =', z)
2
y2 + 1
r = 2227.0211 y = 47.19 e z = (y**2 + 1)**2
Assim, atribuímos dados às variáveis Python. O Python diferencia maiúsculas de minúsculas e nomes
como X e x são diferentes. No Python as variáveis são ponteiros (pointers). Comandos como x = 2 cria um
objeto x e atribuí (armazena) o valor 2 nele em outras linguagens, mas no Python, há um objeto inteiro 2
e x é um ponteiro, apontando para ele. Veja as consequências disso a seguir, sendo que o comando \n
realiza uma quebra de linha. Não há maiores implicações em objetos escalares como este, mas quando se
trata, por exemplo, de listas, nosso próximo objeto, a questão já é bem diferente.
x = z = b = 1
b = 7
print('x is pointing to', x,
'\nz is pointing to', z, '\nb is pointing to', b)
x is pointing to 1
z is pointing to 1
b is pointing to 7
Para lidarmos com funções de números complexos a library cmath deve ser importada e as funções
trigonométricas de números complexos podem ser usadas com base nesta biblioteca e não na library math,
1.2. VARIÁVEIS BOOLEANAS 7
que é designada para números reais (float). Veja o script ilustrativo a seguir.
import cmath
([Link](0.1 - 0.4j) + [Link](0.2 + 0.6j))**2.5
(1.0143106793360148+2.4164569923716543j)
True
False
True
True
1.3 Strings
As strings (variáveis texto) são um importante tipo de objeto Python. Uma vez que temos um objeto
definido, os métodos e funções estão disponíveis para serem usados. As strings são denotadas por str em
Python. Veja alguns exemplos, em que as strings foram atribuídas ou não a objetos (variáveis).
'Esta é uma string'
mensagem = 'Universidade Federal '
type(mensagem)
mensagem
UFLA = mensagem + 'de Lavras, MG.'
UFLA
Se o usuário entrar com Daniel, o resultado será Daniel foi aprovado!. Esta versão de Markdown ainda
não suporta interatividade com o usuário. Portanto, o comando input não foi avaliado na saída deste
script. No segundo comando, se o usuário entrar com um número não inteiro, haverá uma mensagem de
erro do Python. Existem opções para lidar com erros deste tipo e de outras causas também.
1.4.1 Listas
As listas, lists são os primeiros blocos de construção para lidarmos como manipulação de dados. As
listas são vetores cujo primeiro elemento inicia-se no 0, mas cujos elementos de cada célula pode ser
diferentes tipos mistos, desde inteiros, booleanos, reais, complexos, strings, caracteres, conjuntos, tuplas
e outras listas. As listas fazem parte do quarteto list, tuple, set e dictionary. A biblioteca numpy
fornece ferramentas adicionais para lidarmos com grande coleções de dados.
x = [1, 2, 3, 4]
x
y = [1, 'Estat', 3.5, 4+5j]
y
type(x)
type(y)
y[3]
[1, 2, 3, 4]
[1, 'Estat', 3.5, (4+5j)]
1.4. LISTAS, TUPLAS, CONJUNTOS E DICIONÁRIOS 9
list
list
(4+5j)
A variável x é uma lista de inteiros com 4 elementos, que são indexados por 0, 1, 2, 3. Assim, x[1] aponta
para o valor 2 e x[0] para o valor 1. A variável y também é uma lista com 4 elementos de diferentes tipo,
sendo y[0] um inteiro, y[1] uma string y[2] um float e y[3], um número complexo. Para criar a lista,
simplesmente utilizamos as chaves [], com cada elemento da lista separado por uma vírgula. É possível
criar uma lista com elementos com valores repetidos e eles serão identificados como sendo diferentes, pois
a lista respeita as ordens de entradas dos valores e preserva a ordem. Podemos verificar se um elemento
pertence a lista com o comando in, como mostra o script a seguir, entre outros exemplos.
[2, 7] == [7, 2]
[5, 7] == [5, 7, 7]
x
2 in x
7 in x
y
4+5j in y
'Daniel' in y
False
False
[1, 2, 3, 4]
True
False
[1, 'Estat', 3.5, (4+5j)]
True
False
Podemos, como foi feito com as strings realizar algumas operações aritméticas com as listas, como
mostra o exemplo do seguinte script.
x+y
x+[[0,1],'teste',[1,0]]
y*2
y[0] # primeiro elemento da lista
y[-1] # último elemento da lista
# numeração -1,-2,-3,-4 para o índice
# acessa as posições, 3,2,1,0,
# respectivamente da lista y
Vejamos agora o problema dos ponteiros, por meio do exemplo do script apresentado na sequência.
x = z = b = [1,2,3]
b[1] = 7
print('x is pointing to', x,
'\nz is pointing to', z, '\nb is pointing to', b)
# todos os objetos foram alterados e não só b
# pois eles apontam para a mesma lista [1,2,3]
x is pointing to [1, 7, 3]
z is pointing to [1, 7, 3]
b is pointing to [1, 7, 3]
Observamos que se x, z e b apontarem para o mesmo objeto, então se alterarmos o valor b[1] de 2 para 7,
então todos os três objetos serão alterados na posição 1, que corresponde ao segundo valor da lista, pois
ela se inicia na posição 0. Entretanto, se em vez de b[1] = 7 tivéssemos usado a atribuição b = [7,9],
então os vetores x e z não seriam alterados, com a nova atribuição do vetor b. Nos exemplos anteriores,
vimos também que os elementos de uma lista são acessados pelo seu índice que varia de 0 a n-1, sendo n o
seu tamanho. Assim, a lista x=[1,2,3,4] tem seus elementos x[0] igual a 1, x[1] igual a 2, x[2]igual a
3 e x[3] igual a 4. Também podemos variar o índice de -1 a -n, sendo que -1 significa a última posição
da lista, ou seja, a posição n-1, -2 corresponde a posição n-2 e assim por diante até -n, que corresponde
a posição 0 da lista.
As listas são objetos e como tais podemos utilizar alguns métodos associados a eles. As listas são mutáveis
e dinâmicas (podemos alterar seus elementos), são ordenadas (cada elemento da lista possui uma ordem
definida na sua criação) e permitem elementos repetidos. Para usarmos um método ou uma função
deveremos considerar a diferença entre eles. Embora todos métodos sejam funções em Python, nem toda
função é um método. As funções recebem os objetos como entradas e não os modifica e os métodos agem
nos objetos. A seguir apresentamos uma relação de alguns métodos ou funções associados às listas:
• sort(): ordena a lista em ordem crescente.
• append(): adiciona um elemento ao final da lista.
• extend: adiciona múltiplos elementos à lista.
• index(): usado para encontrar o índice de um elemento na lista.
• max(list): retorna o valor máximo de uma lista.
• min(list): retorna o valor mínimo de uma lista.
• list(tuple): transforma uma tuple numa lista.
• len(list): retorna o tamanho da lista (número de elementos).
• filter(fun,list): filtra uma lista usando uma função fun Python.
Vamos ilustrar alguns destes métodos com exemplos particulares. Vamos considerar uma lista e aplicarmos
o método sort() para ordenarmos os seus valores. Neste exemplo a seguir, vamos ver a diferença de um
método e de uma função, observando como o método modifica o objeto que o chamou. Neste caso, a
chamada de um método é dada pelo nome da lista (objeto) seguida de um ponto e do nome do método:
[Link]().
x = [7.4, 5.8, 9.3, 3.2]
x # objeto x original
[Link]()
x # objeto x modificado pelo método sort() ordenado
3
Já as funções len, max() e min() atuam no objeto, passado como entrada da função, mas não o modificam.
Veja o exemplo na lista x dos exemplos anteriores o efeito destas duas funções.
len(x) # tamanho da lista x
len(mes) # tamanho da lista mes
b = max(x)
a = min(x)
a
b
mes
x # x e mes não modificados
6
4
1.6
11.6
['Janeiro', 'Fevereiro', 'Março', 'Abril']
[1.6, 3.2, 5.8, 7.4, 9.3, 11.6]
Para ilustrar a uso da função filter() vamos considerar uma função que retorna True ou False para
uma certa condição de interesse. Por exemplo, se quiséssemos saber quais números dos seis elementos da
lista x possui resto da divisão por 2 menor que 1,5. Esse resultado é obtido com a comparação a % 2 <=
1.5, que irá retornar verdadeiro ou falso para o número representado por a. Só que devemos fazer isso
para todos os elementos da lista x ou de outra lista qualquer. Devemos criar uma função para receber
12 CHAPTER 1. INTRODUÇÃO AO PYTHON
cada elemento da lista e verificar a condição, retornando True ou False e passar pela função filter()
para realizar a iteração nos elementos da lista x ou na lista de interesse. Vamos criar um primeira função
no exemplo a seguir e em seguida aplicar a função filter(). O tipo de objeto retornado desta função é
filter, logo, tem de ser transformado em lista antes de imprimir.
def resto(a):
if ((a % 2) <= 1.5):
return True
else:
return False
# aplicar a função filter
x_filtrado = filter(resto, x)
print(list(x_filtrado))
38.9
6.483333333333333
Outros métodos como o count() (conta o número de ocorrências de um dado valor), reverse() (ordena
a lista em ordem reversa a ordem original), clear() (limpa todos os dados da lista), copy() (copia todos
os dados da lista), insert() (insere um elemento em uma posição específica da lista) e pop() (remove
um elemento em uma posição específica) como apresentado no script a seguir.
L = [1, 2, 3, 5, 7, 2, 4.5]
[Link](2)
[Link](2.1)
L1 = [Link]()
[Link]()
L
[Link]()
L
L1
1.4. LISTAS, TUPLAS, CONJUNTOS E DICIONÁRIOS 13
2
0
[4.5, 2, 7, 5, 3, 2, 1]
[]
[1, 2, 3, 5, 7, 2, 4.5]
[1, 3, 2, 3, 5, 7, 2, 4.5]
3
[1, 2, 3, 5, 7, 2, 4.5]
Podemos construir uma matriz, haja vista que Python não possui um objeto matricial, usando uma lista.
Se criarmos uma lista de n componentes, com cada um dos componentes tendo m componentes, teremos
uma matriz n × m. Vejamos no script a seguir a construção da seguinte matriz:
1 4
A= 2 5 . (1.1)
3 6
[1, 2, 3, 4, 5, 6, 7]
Acessando posições em particulares, para impressão ou para atribuição:
L[0]
L[1] = 8
L
1
[1, 8, 3, 4, 5, 6, 7]
Acessando, posições com os índices negativos.
14 CHAPTER 1. INTRODUÇÃO AO PYTHON
7
1
Para blocos de elementos temos que o comando L{m:s:r] acessa os elementos nas posições m, m+1+r,
m+1+2r, ... até a (s-1)-ésima posição ou até a posição mais próxima de s-1 possível. Muito cuidado
deve ser tomado, pois o limite, superior não indica onde o subconjunto termina entre os índices válidos de
uma lista e sim,que ela termina na posição destacada subtraída de 1.
L[1:4]
L[1:23] # passa do limite len(L)
L[0:5:2]
L[-1:-8:-1]
L[4:]
L[:6]
[8, 3, 4]
[8, 3, 4, 5, 6, 7]
[1, 3, 5]
[7, 6, 5, 4, 3, 8, 1]
[5, 6, 7]
[1, 8, 3, 4, 5, 6]
Se omitirmos os limites inferior ou superior da sequência, então a lista selecionada será iniciada no índice
0 (valor inicial) ou terminará no último índice (valor final da lista), como nos dois últimos exemplos
apresentados.
1.4.2 Tuplas
As tuplas são objetos Python muito parecidos com as listas. Vários métodos e funções que se aplicam
às listas também se aplicam às tuplas. Ao contrário das listas, as tuplas são objetos imutáveis, ou seja,
uma vez criadas elas não podem ser modificadas. Assim, se criarmos uma tupla por t = (1,2,3) não
poderemos atribuir valor, por exemplo, deste jeito t[1] = 9. Elas podem conter mais de um valor idêntico
e são ordenadas, como as listas. A forma de criar a tupla em relação à lista é o uso dos parênteses no
lugar dos colchetes.
t = (1, 2, 'DFF', 3)
'DFF' in t
t[2]
t[0]
t[:3]
len(t)
True
'DFF'
1
(1, 2, 'DFF')
4
1.4. LISTAS, TUPLAS, CONJUNTOS E DICIONÁRIOS 15
Os métodos count() e index podem ser usados nas tuplas, como ilustrado a seguir. O método index()
tem a seguinte sintaxe, sendo que os dois últimos argumentos são opcionais: [Link](element,
start, end).
[Link](1) # número de ocorrência de 1
[Link]('DFF') # índice da posição de 'DFF'
1
2
A tupla pode ter qualquer tipo como sendo seus elementos, incluindo uma tupla ou uma lista.
t1 = ((1,2),'r', [2,3,4])
t1
print('Componente lista da tupla ',t1[2])
print('Elemento 0 do componente lista da tupla ',t1[2][0])
t1[2].append(5)
print('Modificando o componente lista da tupla ',t1)
3
8
A função any(t) retorna True se há algum item True na tupla e retorna False, caso contrário. Neste caso,
a tupla ou qualquer outro tipo apropriado poderá ter elementos 0 e 1 ou booleanos. Pode ser aplicada nas
listas, conjuntos e nos dicionários.
t = (True, False, False,False,True)
any(t)
True
Podemos usar ainda as funções min(), max(), sum() e sorted(), como ilustrado a seguir. A função
sorted() ordena a tupla e retorna uma lista ordenada como resultado. Veja que o método [Link]()
altera o objeto lista e não pode ser aplicado na tupla, pelo fato de a tupla ser imutável.
t = (2.3,4.5,1.2,1.1,9.7,5.3)
min(t)
max(t)
sum(t)
16 CHAPTER 1. INTRODUÇÃO AO PYTHON
sorted(t)
t
1.1
9.7
24.099999999999998
[1.1, 1.2, 2.3, 4.5, 5.3, 9.7]
(2.3, 4.5, 1.2, 1.1, 9.7, 5.3)
1.4.3 Conjuntos
Os conjuntos set em Python tem uma conotação muito próxima com a definição de conjuntos da
matemática. Esses são conjuntos que a ordem ou duplicação de seus elementos não mudam o conjunto.
Assim, são imutáveis, não ordenados e não pode ter mais de um elemento idêntico em suas ocorrências.
Podemos usar o construtor (função) set() para criar um conjunto ou usarmos as chaves{} para digitar
seus elementos separados por vírgula.
phi = set()
print('Conjunto vazio: ', phi)
# precisa ser uma lista ou tupla de argumento
A = set(['A','D','B','C','E'])
A
B = {1,2,3.4,5,6,6}
'A' in A
B # repare que o elemento
#repetido 6 aparece 1 vez apenas
B[0]
True
Algumas operações matemáticas com conjuntos estão disponíveis em Python, como união, interseção,
diferença (Ac ∩ B) e diferença simétrica ((Ac ∩ B) ∪ (A ∩ B c )), como ilustrados no exemplo a seguir.
A = {1,2,3,4,5,6}
B = {4,5,7,8,9,10}
[Link](B)
[Link](B)
[Link](B) # esta em A, mas não em B
[Link](A) # está em B, mas não em A
A.symmetric_difference(B) # está só em A ou só em B
A & B # intersecção
A | B # união
A - B # diferença
B - A # diferença
AˆB # diferença simétrica
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
{4, 5}
{1, 2, 3, 6}
{7, 8, 9, 10}
{1, 2, 3, 6, 7, 8, 9, 10}
{4, 5}
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
{1, 2, 3, 6}
{7, 8, 9, 10}
{1, 2, 3, 6, 7, 8, 9, 10}
1.4.4 Dicionários
Os dicionários dictionary() são objetos mutáveis e iteráveis. Os seus elementos vem sempre aos pares,
sendo o primeiro valor uma chave e o segundo elemento, o valor da chave. Tanto a chave quanto o seu
valor são objetos Python. A chave é imutável, mas seus valores associados são objetos mutáveis ou não.
D = {1: [1,2,3.4,5], 2: 3.7, 3: {1,2,3}}
D[1]
D[2]
type(D[3])
type(D[1])
[1, 2, 3.4, 5]
3.7
set
list
Temos um dicionário, com as chaves 1, 2 e 3. Para a chave 1 temos uma lista como seu valor, para a chave
2, temos um valor float e para a chave 3, criamos um objeto do tipo conjunto. A seguir, acrescentamos
18 CHAPTER 1. INTRODUÇÃO AO PYTHON
uma chave, nomeada 4 com um valor booleano associado. O método keys() recupera as chaves do objeto
dicionário D, transforma numa lista e imprime a lista e seu primeiro elemento C[0].
D[4] = True
print(D)
C = list([Link]())
print(C)
C[0]
Em seguida, criamos uma matriz com o uso da função array(). Vamos usar nossa lista A anterior, para
fazer isso.
B = [Link](A)
B
array([[4, 1],
[1, 2]])
Podemos criar também a partir de tuplas, em vez de listas, a matriz numpy. Além disso, existem funções
próprias do pacote para criarmos matrizes, como, por exemplo, a matriz de zeros 2 × 4 a seguir. Também
existem funções para criarmos arrays (vetores) unidimensionais, como o método arange() e o linspace().
O primeiro cria um vetor indo de n até o máximo m (inteiros) sem incluí-lo de 1 em 1, ou do mínimo n até
20 CHAPTER 1. INTRODUÇÃO AO PYTHON
[[0. 0. 0. 0.]
[0. 0. 0. 0.]]
array([2, 3, 4, 5, 6])
array([2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. , 6.5])
array([2. , 2.8, 3.6, 4.4, 5.2, 6. ])
Vamos ilustrar alguns cálculos simples com vetores.
x = [Link](2,3,0.2)
y = [Link](3,6,0.7)
print('', x,'\n', y)
x + y # adição dos vetores
x * y # produto elementwise
y ** x # potenciação elementwise
array([[4, 1],
[1, 2]])
array([[2, 1],
[1, 2]])
array([[9, 6],
[4, 5]])
array([[9, 6],
[4, 5]])
array([[4, 1],
1.5. MATRIZES E ARRANJOS 21
[1, 2]])
As inversas podem ser obtidas com a função [Link]() e as inversas generalizadas de Moore-
Penrose pelo método [Link](). O determinante pode ser obtido por [Link]() e os
autovalores e autovetores (decomposição espectral) por [Link](). É importante observar que os
autovalores do [Link]() não necessariamente estará ordenado do maior para o menor, como é
convencionalmente adotado em diferentes outros programas. Este método é bem geral e pode ser usado
em matrizes reais ou complexas quadradas. Alternativamente, para matrizes simétricas o numpy possui o
método [Link](), que retorna os autovalores em ordem crescente. Veja o exemplo a seguir.
[Link](B)
[Link](B) # igual a inversa (posto completo)
[Link](B)
L, P = [Link](B)
print('Autovalores: ',L,'\nAutovetores: ',P)
L1, P1 = [Link](B)
print('Autovalores: ',L1,'\nAutovetores: ',P1)
[-0.78116852 0.62413562]]
Vetor dos valores singulares: [11.19813546 5.88232625]
Vetores singulares à direita (transposto): [[-0.75238412 -0.65872463]
[ 0.65872463 -0.75238412]]
2
array([[11.19813546, 0. ],
[ 0. , 5.88232625]])
Para verificarmos que a decomposição realmente é adequada, temos o seguinte script:
[Link]([Link](U,[Link](L)),Vt)
[Link]([Link](L)).dot(Vt) # alternativo
U @ [Link](L) @ Vt # alternativo
array([[1., 4.],
[2., 7.],
[9., 3.]])
array([[1., 4.],
[2., 7.],
[9., 3.]])
array([[1., 4.],
[2., 7.],
[9., 3.]])
Muitas outras funções existem no pacote numpy para lidarmos ou operarmos matrizes. Se necessitarmos
de alguma outra função para alguma operação matricial, vamos apresentá-la nestas ocasiões.
X1 X2
0 13.4 14
1 14.6 15
2 13.5 19
3 15.0 23
4 14.6 17
5 14.0 20
6 16.4 21
7 14.8 16
8 15.2 27
9 15.5 34
10 15.2 26
11 16.9 28
12 14.8 24
13 16.2 26
14 14.7 23
15 14.7 9
16 16.5 18
17 15.4 28
18 15.1 17
19 14.2 14
Os DataFrames são estruturas de dados tabular, sendo que cada coluna possui constitui de uma sequência
de valores do mesmo tipo (booleano, float, strings, etc.), em que as diferentes colunas podem ser e
potencialmente são de diferentes tipos. Os DataFrames possuem uma coluna adicional chamada index
que no exemplo anterior, do objeto dados, variou de 0 a 19, pois nosso DataFrame possui 20 linhas e duas
variáveis X1 e X2, que no arquivo [Link] estavam identificadas na primeira linha física do arquivo que
foi lido pelo método read_csv(). O index mapeia as linhas do DataFrame com os labels mencionados.
Vamos mostrar como construir um DataFrame diretamente a partir de um objeto dict para o construtor
DataFrame() do pandas. Vamos criar um DataFrame de um delineamento inteiramente casualizado, com
2 tratamentos e 3 repetições de cada um, com os respectivas produtividades avaliadas nas 6 parcelas
experimentais.
dic = {'rep': [1,2,3,1,2,3],
'trat': [1,1,1,2,2,2],
'prod': [3.4,2.3,5.6,5.7,6.3,7.1]}
arqd = [Link](dic)
arqd
Podemos acrescentar uma nova coluna em nosso DataFrame, ou eliminarmos uma já existente, criando um
novo DataFrame para receber o resultado, como mostrado a seguir, com a opção columns, para qualquer
ordem das chaves (nomes das colunas).
24 CHAPTER 1. INTRODUÇÃO AO PYTHON
arqd
arqd1 = [Link](arqd,columns=['trat','prod'])
arqd1 # selecionando 2 variáveis no DataFrame
arqd['alt'] = [1.5,1.3,1.4,1.2,1.1,1.6] #criando variável altura
arqd
trat prod
0 1 3.4
1 1 2.3
2 1 5.6
3 2 5.7
4 2 6.3
5 2 7.1
Para acessarmos as chaves (colunas), os índices e os valores do DataFrame podemos usar os seguintes
códigos ilustrativos. Observamos que a instrução [Link] para acessar a coluna prod, resulta em erro,
pois prod é uma palavra reservada (produto). Devemos usar a alternativa anterior para esta chave.
[Link]
[Link]
[Link]
arqd['prod'][0] # produção do índice 0
arqd['prod'] # toda a coluna de produção
[Link] # cuidado, pois prod é palavra reservada: erro
[Link]
3.4
0 3.4
1 2.3
2 5.6
3 5.7
4 6.3
5 7.1
Name: prod, dtype: float64
0 1
1 2
2 3
3 1
4 2
5 3
Name: rep, dtype: int64
Para extrairmos uma linha inteira usamos o método loc do DataFrame associado ao índice da linha
(registro), como ilustrado no script a seguir.
L = [Link][1] # segunda linha do DataFrame
print(L)
L2 = [Link][[0,5]] # as linhas 1 e 6 de arqd
L2 # com os índices 0 e 5
rep 2.0
trat 1.0
prod 2.3
alt 1.3
Name: 1, dtype: float64
Para selecionar um bloco de registros (linhas) indo de inicio ao fim (excluindo o limite final) usamos
arqd[inicio:fim]. Como ilustrado a seguir, onde extraímos do índice 0 até o índice 3, ou seja, as três
primeiras linhas com os índices 0, 1 e 2 do arqd.
26 CHAPTER 1. INTRODUÇÃO AO PYTHON
arqd[0:3]
Valores perdidos podem fazer parte do DataFrame e neste caso, eles assumem o valor NaN, do inglês not
a number. Podemos deletar uma coluna com o comando del, da seguinte forma.
del arqd['alt']
arqd
Podemos filtrar impondo condições específicas ao DataFrame. Por exemplo, se estivermos interessado no
DataFrame resultante dos elementos em que a produção é maior ou igual a 5, teremos o seguinte resultado.
O resultado filtrado mostra os registros nas mesmas posições originais e o seu DataFrame original arqd
permanece inalterado.
arqd[arqd['prod'] >= 5.0]
arqd
Para gravarmos um DataFrame podemos escolher o diretório (pasta) e o nome do arquivo e gravarmos em um
arquivo csv (arquivo separado por vírgula) com o comando arqd.to_csv('[Link],index=False,header=True)
para não salvar o índice e salvar o cabeçalho. Vamos recuperar em nosso código, o path original com o
1.7. ESTRUTURAS DE CONTROLE DE PROGRAMAÇÃO 27
comando [Link](apath), em que apath foi obtido quando iniciamos o assunto sobre DataFrame e
refere-se ao diretório deste projeto. Escolhemos o nome [Link] para o arquivo. Podemos ler o arquivo
novamente, conforme mostramos nos primeiros passos da abordagem dos DataFrames. Depois de gravado,
repetimos sua leitura e o colocamos no objeto dic, em que devemos atentar para o separador de colunas,
que neste caso é a vírgula.
[Link](apath)
arqd.to_csv('[Link]',index=False,header=True)
dic = pd.read_csv('[Link]',sep=r',')
dic
Em futuras edições, mostraremos mais detalhes dos DataFrames. Nos capítulos posteriores, caso venhamos
a precisar de um DataFrame e de algumas de suas propriedades, então iremos adicionar os conteúdos
necessários nesta ocasião. Falaremos agora das estruturas condicionais e as estruturas de repetições.
A instrução elif significa else if e é opcional. Se a condição1 for verificada, é executado o bloco
denominado instruções1, que estão indentados em relação ao if. Caso a condição seja falsa, é testada a
condição2 e se ela for verdadeira, são executados as instruções denotadas por instruções2, que pode ser
uma simples instrução ou várias instruções, indentadas em relação ao elif. Finalmente, se a condição2
for falsa são executadas as instruções3. Cada linha das instruções if, elif ou else são seguidas por
dois pontos. Veja o exemplo simples a seguir.
x = 5
if x > 6:
print('Recebe mais que 6 salários.')
print('Você está entre os 20% mais ricos!')
else:
print('Você recebe 6 salários ou menos.')
print('Você representa 80% da população!')
0 se x < 0
FX (x) = x2 se 0 ≤ x ≤ 1
1 se x > 1.
Para este modelo, temos o seguinte script, no qual decidimos qual parte do programa rodar, conforme os
valores de x são atribuídos.
x = 0.8
if x < 0:
F = 0
elif 0 <= x <= 1:
F = x**2
else:
F = 1
F
0.6400000000000001
As estruturas de repetição do Python são o for e o while/else. Existe ainda um terceiro tipo de
procedimento em Python para realizarmos iterações. A estrutura geral de um código Python para o for é
apresentada no script a seguir.
x = [1,2,3,4,5,6]
for i in x:
instruçoes1
else: # opcional instrução else
instruções2
Os objetos (x) do comando for são os objetos iteradores ou iteráveis, que são aqueles que contêm um
número contáveis de valores. As listas, tuplas, dicionários e conjuntos ão todos objetos iteráveis. Vamos
ilustrar com um simples exemplo a seguir, para calcularmos a soma, o produtório e a média de uma lista
de valores.
1.7. ESTRUTURAS DE CONTROLE DE PROGRAMAÇÃO 29
A soma é: 14.899999999999999
O produtório é: 152.90744999999995
A média é: 2.9799999999999995
Para realizarmos iterações em um dicionário, temos o seguinte exemplo:
D = {1: 1.3, 2: 3.1, 3: 1.7}
for i in D: # iterar nas chaves
print(i, ' e ', D[i])
print('Agora iterando nas chaves e valores:')
for (i, valor) in [Link]():
print(i, ' e ', valor) # iterar em chave e valor
1 e 1.3
2 e 3.1
3 e 1.7
Agora iterando nas chaves e valores:
1 e 1.3
2 e 3.1
3 e 1.7
Finalmente, um exemplo em um conjunto:
A = {1.1,2.2,3.4,4.7,5.3}
for i in A:
print('elemento: ',i)
elemento: 1.1
elemento: 2.2
elemento: 3.4
elemento: 4.7
elemento: 5.3
Podemos criar um sequência de valores com o comando range(n) que vai de 0 a n-1. Assim, também
podemos usar o for nessa sequência, como ilustrado no exemplo a seguir.
x = [2.3, 4.1, 1.5, 2.3, 4.7]
soma = 0
n = len(x)
for i in range(n):
soma = soma + x[i]
media = soma / n
print('A soma é: ',soma)
30 CHAPTER 1. INTRODUÇÃO AO PYTHON
A soma é: 14.899999999999999
A média é: 2.9799999999999995
O while é uma outra estrutura de repetição, em que o bloco de comandos indentados irão ser executados
até que uma condição seja satisfeita. A estrutura geral é dada a seguir.
while condição:
instruçoes1
else: # opcional instrução else
instruções2
A soma é: 14.899999999999999
A média é: 2.9799999999999995
Podemos usar os comandos break e continue dentro do while (ou do for). O break é usado após uma
segunda condição ser verificada no bloco de comandos do interior da estrutura de repetição e pula a
execução do programa para a primeira linha de instrução após o bloco do loop, ou seja, encerra o loop.
O continue executa a primeira linha testando a condição primária do while ou tomando o próximo valor
do iterador no for, ou seja, vai para o início do loop. Veja o exemplo a seguir.
y = 35 # experimente outro número inteiro > 1
x = y // 2
while x > 1:
if y % x == 0:
print(y, 'tem fator ', x)
break
x = x - 1
else:
print(y,' é primo')
35 tem fator 7
Podemos utilizar a função filter, como já ilustramos anteriormente neste material, para realizarmos
iterações em nosso código. Não daremos mais detalhes disso, por enquanto.
1.8 Funções
As funções em todas as linguagens são uma poderosa ferramenta de programação, que nos permite quebrar
um grande problema em pequenas tarefas (as funções), facilitando assim a resolução do problema como
um todo. Dizemos que é a estratégia de dividir para conquistar. As funções em geral recebem um objeto
1.8. FUNÇÕES 31
e o processa de acordo com as regras definidas em seu bloco de comando. Desta forma a linguagem
ganha grande poder, conveniência e elegância. O aprendizado em escrever funções úteis é uma das muitas
maneiras de fazer com que o uso do Python seja confortável e produtivo. A sintaxe geral de uma função é
dada por:
def nome(arg1, arg2,...,argn):
instruções
As instruções significam um bloco de comandos (indentados) e podem ou não ter o comando return
objeto, que pode acontecer em qualquer parte do bloco de comandos. A função pode não ter este
comando de return, se ela modificar um arquivo apenas gravando um novo resultado ou se imprimir uma
mensagem quando chamada. Os argumentos ou parâmetros são passado para a função e a sua chamada
deve obedecer estritamente a ordem em que eles aparecem, a menos que a chamada seja com chave, ou
seja, do tipo arg1 = 2.3, por exemplo. Neste caso, os argumentos podem ser colocados em qualquer
ordem. Os argumentos de uma função podem conter valores default, ou seja, na declaração do nome
da função podemos ter algo do tipo: def nome(x, theta = 0.5). O argumento theta=0.5 pode ser
omitido na chamada da função, que será atribuído seu valor 0,5 padrão.
Vamos apresentar uma função simples para testar a hipótese H0 : µ = µ0 a partir de uma amostra simples
de uma distribuição normal. Dois argumentos serão utilizados: o vetor (lista) de dados x de tamanho n e
o valor real hipotético µ0 . A função calculará o valor da estatística tc do teste por:
X̄ − µ0
tc = . (1.2)
√S
n
A função resultante, em Python, é apresentada a seguir. Neste exemplo uma amostra de tamanho n = 8 foi
utilizada para obter o valor da estatística para testar a hipótese H0 : µ = 3,0 (se a amostra era proveniente
do povo na’vu). Podemos observar que o resultado final da função é igual ao do último comando executado,
ou seja o valor da estatística e do valor-p, por meio de um objeto do tipo dicionário. Esta função utiliza no
seu escopo três funções do Python (funções básicas do numpy), ainda não apresentadas. As duas primeiras,
var() e mean() retornam a variância e a média do vetor utilizado como argumento, respectivamente, e a
terceira, pt(), retorna a probabilidade acumulada da distribuição t de Student para o primeiro argumento
da função com ν graus de liberdade, que é o seu segundo argumento.
import scipy as sp # para calcular probabilidade da t
def t_test(x, mu0):
n = len(x)
s2 = [Link](x,ddof=1) # ddof=1, divisor n-1 para s2
xb = [Link](x)
t = {'tc':0,'[Link]':0}
t['tc'] = (xb-mu0) / (s2 / n)**0.5
t['[Link]'] = 2*([Link](abs(t['tc']),n-1))
return t
y = [1.76,1.81,1.74,1.71,1.79,1.75]
t = t_test(y, 3.0) # altura de avatares
print('tc = ',t['tc'])
print('[Link] = ',t['[Link]'])
tc = -84.89699641330068
[Link] = 4.297311617662558e-09
Podemos reescrever esta função para colocarmos um valor default para o argumento mu0. Se escolhêssemos
o valor 0 e em alguma chamada da função, esse argumento fosse omitido, seria feito o teste da hipótese
H0 : µ = 0 por padrão.
32 CHAPTER 1. INTRODUÇÃO AO PYTHON
tc = 120.49896265113645
[Link] = 7.465636997494585e-10
[('S2', 0.001280000000000002), ('xbar', 1.76)]
Vamos realizar um teste para a correlação em populações normais bivariadas. Assim, dado o par de
variáveis vetoriais x e y tomados em n indivíduos, temos a hipótese nula
H0 : ρ = 0
√
r n−2
tc = √ ,
1 − r2
tc = 58.469836352468235
[Link] = 0.0002923787083537466
r = 0.9997076212916461
Finalmente, vamos obter uma função para obtermos potências reais de matrizes quadradas simétricas
positivas definidas. Para isso vamos escrever uma função para obter potências reais de uma matriz A
simétrica e positiva definida:
A =PΛP⊤ ,
Aα =PΛα P⊤ .
Deve receber A e retornar Aα . O script a seguir ilustra uma função para obtermos estas potências
matriciais. Observe que não é uma potência elemento a elemento. Também devemos observar que não
é importante que os autovalores estejam ordenados. Mas é óbvio que os autovetores associados a cada
autovalor deve estar corretamente associado e preservado e isso é feito pela biblioteca numpy por meio da
função eig().
def mat_power(A, alpha = 0.5):
e_val, e_vec = [Link](A)
if any(e_val) < 0:
print('Matriz não é positiva definida!')
return
Ap = e_vec.dot([Link](e_val**alpha)).dot([Link](e_vec))
return Ap
A = [[4,1],[1,2]]
print('A = ',A)
mat_power(A) # raiz quadrada
A3 = mat_power(A, 1/3) # raiz cúbica
A3
[Link](A3).dot(A3) # verificando
deste material em uma seção básica e em outra aplicada. As técnicas computacionais são denominadas
de estatística computacional se forem usadas para realizarmos inferências, para gerarmos realizações de
variáveis aleatórias ou para compararmos métodos e técnicas estatísticas.
Vamos explorar métodos de geração de realizações de variáveis aleatórias de diversos modelos probabilísticos,
para manipularmos matrizes, para obtermos quadraturas de funções de distribuição de diversos modelos
probabilísticos e de funções especiais na estatística e finalmente vamos apresentar os métodos de computação
intensiva para realizarmos inferências em diferentes situações reais. Temos a intenção de criar algoritmos
em linguagem Python e posteriormente, quando existirem, apresentar os comandos para acessarmos os
mesmos algoritmos já implementados.
Simular é a arte de construir modelos segundo Naylor et al. [1971], com o objetivo de imitar o funcionamento
de um sistema real, para averiguarmos o que aconteceria se fossem feitas alterações no seu funcionamento
(Dachs [1988]). Este tipo de procedimento pode ter um custo baixo, evitar prejuízos por não utilizarmos
procedimentos inadequados e otimizar a decisão e o funcionamento do sistema real.
Precauções contra erros devem ser tomadas quando realizamos algum tipo de simulação. Podemos
enumerar:
3. erros de implementação.
Devemos fazer o sistema simulado operar nas condições do sistema real e verificar por meio de alguns
testes se os resultados estão de acordo com o que se observa no sistema real. A este processo denominamos
de validação. A simulação é uma técnica que usamos para a solução de problemas. Se a solução alcançada
for mais rápida, com eficiência igual ou superior, de menor custo e de fácil interpretação em relação a
outro método qualquer, o uso de simulação é justificável.
1.10 Exercícios
1. Criar no Python os vetores a⊤ = [4, 2, 1, 5] e b⊤ = [6, 3, 8, 9] e concatená-los formando um único
vetor. Obter o vetor c = 2a − b e o vetor d = b⊤ a. Criar uma sequência cujo valor inicial é igual
a 2 e o valor final é 30 e cujo passo é igual a 2. Replicar cada valor da sequência 4 vezes de duas
formas diferentes (valores replicados ficam agregados e a sequência toda se replica sem que os valores
iguais fiquem agregados).
3. Criar a matriz
1.10. EXERCÍCIOS 35
10 1
A=
1 2
Neste capítulo vamos considerar a geração de números aleatórios para o modelo probabilístico uniforme.
A partir do modelo uniforme podemos gerar realizações de variáveis aleatórias de qualquer outro modelo
probabilístico. A geração de realizações de uma distribuição uniforme não pode ser realizado por máquinas.
Qualquer sequência produzida por uma máquina é uma sequência previsível de números pseudo-aleatórios.
Dois geradores de números aleatórios devem produzir os mesmos resultados nas suas aplicações. Se isso
não ocorrer, um deles não pode ser considerado um bom gerador de números aleatórios [Press et al., 1992].
Os conceitos de números uniformes e números aleatórios são muitas vezes confundidos. Números uniformes
são aqueles que variam aleatoriamente em um intervalo real de valores com probabilidade constante. No
entanto, devemos diferenciar números aleatórios uniformes de outros tipos de números aleatórios, como,
por exemplo, números aleatórios normais ou gaussianos. Estes outros tipos são geralmente provenientes
de transformações realizadas nos números aleatórios uniformes. Então, uma fonte confiável para gerar
números aleatórios uniformes determina o sucesso de métodos estocásticos de inferência e de processos de
simulação Monte Carlo.
37
38 CHAPTER 2. VARIÁVEIS ALEATÓRIAS UNIFORMES
U1 = (7 × 7 + 7) mod 10 = 56 mod 10 = 6
U2 = (7 × 6 + 7) mod 10 = 49 mod 10 = 9
e assim sucessivamente. Obtemos a sequência de números aleatórios:
{7, 6, 9, 0, 7, 6, 9, 0, 7, 6, 9, · · · }
é tão bom quanto o método congruencial com c = ̸ 0, se o módulo m e o multiplicador a forem escolhidos
com cuidado [Press et al., 1992]. Park and Miller [1988] propuseram um gerador “padrão” mínimo baseado
nas escolhas:
a = 75 = 16.807 m = 231 − 1 = [Link] (2.3)
Este gerador de números aleatórios não é perfeito, mas passou por todos os testes a qual foi submetido
e tem sido usado como padrão para comparar e julgar outros geradores. Um problema que surge e que
devemos contornar é que não é possível implementarmos diretamente em uma linguagem de alto-nível a
equação (Equation 2.2) com as constantes de (Equation 2.3), pois o produto de a e Ui excede, em geral, o
limite máximo de 32 bits para inteiros. Podemos usar um truque, devido a Schrage [1979], para multiplicar
inteiros de 32 bits e aplicar o operador de módulo, garantindo portabilidade para implementação em
praticamente todas as linguagens e todas as máquinas. O algoritmo de Schrage baseia-se na fatoração de
m dada por:
em que ⌊z⌋ denota a parte inteira do número z utilizado como argumento. Para um número Ui entre 1 e
m − 1 e para r pequeno, especificamente para r < q, Schrage [1979] mostrou que ambos a(Ui mod q) e
r⌊Ui /q⌋ pertencem ao intervalo 0 · · · m − 1 e que
se verifica. No Python pode-se optar por usar o operador % (mod), que retorna o resto da operação entre
dois inteiros e o operador // (div), que retorna o resultado do dividendo para operações com inteiros.
A quantidade Ui mod q = Ui − (Ui div q) × q pode ser obtida em Python simplesmente por Ui %q.
Atribuímos o resultado a uma variável qualquer definida como inteiro. Para aplicarmos o algoritmo de
Schrage às constantes de (Equation 2.3) devemos usar os seguintes valores: q = 127.773 e r = 2.836.
A seguir apresentamos o algoritmo do gerador padrão mínimo de números aleatórios:
# gerador padrão mínimo de números aleatórios adaptado de Park and
# Miller. Retorna desvios aleatórios uniformes entre 0 e 1. Fazer
# "sem" igual a qualquer valor inteiro para iniciar a sequência;
# "sem" não pode ser alterado entre sucessivas chamadas da sequência
# se "sem" for zero ou negativo, um valor dependente do valor do relógio
# do sistema no momento da chamada é usado como semente. Constantes
# usadas a = 7ˆ5 = 16.807; m = 2ˆ31 - 1 = [Link]
# e c = 0
n = 100000
t1 = [Link]()
x = gnap(n)
t2 = [Link]()
t = t2-t1
print('tempo médio (micros): ',[Link] / n)
print('tempo total (micros): ',[Link])
Existem muitas versões implementadas deste algoritmo, inclusive em Fortran e C e que estão disponíveis
na internet. Felizmente, o Python já possui este algoritmo implementado. Por se tratar de um tópico
mais avançado, que vai além do que pretendemos apresentar nestas notas de aulas, não descreveremos
este tipo de procedimento para incorporações de funções escritas em outras linguagens.
2000
1500
1000
500
0
0.0 0.2 0.4 0.6 0.8 1.0
2.3 Exercícios
1. Utilizar o gerador gnap para gerar n realizações de uma distribuição exponencial f (x) = λe−λx .
Sabemos do teorema da transformação de probabilidades, que se U tem distribuição uniforme,Rx X=
F −1 (U ) tem distribuição de probabilidade com densidade f (x) = F ′ (x); em que F (x) = −∞ f (t)dt
é a função de distribuição de X e F −1 (y) é a sua função
R x inversa para o valor y. Para a exponencial
a função de distribuição de probabilidade é: F (x) = 0 λe−λt dt = 1 − e−λx . Para obtermos a função
inversa temos que igualar u a F (x) e resolver para x. Assim, u = 1 − e−λx e resolvendo para x temos:
x = − ln (1 − u)/λ. Devido à simetria da distribuição uniforme 1 − u pode ser trocado por u. O
resultado final é: x = − ln (u)/λ. Para gerar números da exponencial basta gerar números uniformes
e aplicar a relação x = − ln (u)/λ. Fazer isso para construir uma função que gera n realizações
exponenciais. Aplicar a função para obter amostras aleatórias da exponencial de tamanho n = 100 e
obter o histograma da amostra simulada. Calcule a média e a variância e confronte com os valores
teóricos da distribuição exponencial.
2. Para gerar números de uma distribuição normal, cuja densidade é dada por f (x) =
2.3. EXERCÍCIOS 43
√
1/( 2πσ 2 ) exp{−(x − µ)2 /(2σ 2 )}, qual seria a dificuldade para podermos utilizar o teorema
anunciado no exercício proposto anterior?
3. Como poderíamos adaptar o algoritmo apresentados nesse capítulo para gerar números aleatórios
uniformes utilizando os valores propostos por Park e Miller, ou seja, a = 48.271 e m = 231 − 1?
Implementar o algoritmo, tomando cuidado em relação aos novos multiplicador q e resto r da
fatoração de m?
4. Como você poderia propor um teste estatístico simples para avaliar a aleatoriedade da sequência de
números uniformes gerados por esses algoritmos apresentados no capítulo? Implementar sua ideia.
44 CHAPTER 2. VARIÁVEIS ALEATÓRIAS UNIFORMES
Chapter 3
Neste capítulo vamos apresentar alguns métodos gerais para gerarmos realizações de variáveis aleatórias
de outras distribuições de probabilidade, como, por exemplo, dos modelos exponencial, normal e binomial.
Implementaremos algumas funções em Python e finalizaremos com a apresentação das rotinas otimizadas
e já implementadas.
3.1 Introdução
Vamos estudar a partir deste instante um dos principais métodos, determinado pela lei fundamental
de transformação de probabilidades, para gerarmos dados de distribuições de probabilidades contínuas
ou discretas. Para alguns casos específicos, vamos ilustrar com procedimentos alternativos, que sejam
eficientes e computacionalmente mais simples. Esta transformação tem como modelo fundamental a
distribuição uniforme (0, 1). Por essa razão a geração de números uniformes contínuos é tão importante.
Veremos posteriormente nestas notas de aulas algoritmos para obtermos numericamente a função de
distribuição F (x) e a sua função inversa x = F −1 (p), em que p pertence ao intervalo que vai de 0 a 1.
Este conhecimento é fundamental para a utilização deste principal método.
Neste capítulo limitaremos a apresentar a teoria para alguns poucos modelos probabilísticos, para os quais
podemos facilmente obter a função de distribuição de probabilidade e a sua inversa analiticamente. Para
os modelos mais complexos, embora o método determinado pela lei fundamental de transformação de
probabilidades seja adequado, apresentaremos apenas métodos alternativos, uma vez que, em geral, ele é
pouco eficiente em relação ao tempo gasto para gerarmos cada realização da variável aleatória. Isso se deve
ao fato de termos que obter a função inversa numericamente da função de distribuição de probabilidade
dos modelos probabilísticos mais complexos.
45
46 CHAPTER 3. VARIÁVEIS ALEATÓRIAS NÃO-UNIFORMES
Demonstração: Seja X uma variável aleatória com função de distribuição F e função densidade f . Se
u = F (x), então o jacobiano da transformação é du/dx = F ′ (x) = f (x), em que U é uma variável aleatória
uniforme U (0, 1), com função densidade g(u) = 1, para 0 < u < 1 e g(u) = 0 para outros valores de u.
Assim, a variável aleatória X = F −1 (U ) tem densidade f dada por:
du
fX (x) = g(u) = g (FX (x)) f (x) = f (x).
dx
Em outras palavras a variável aleatória X = F −1 (U ) possui função densidade fX (x), estabelecendo o
resultado almejado e assim, a prova fica completa. ■
Para variáveis aleatórias discretas, devemos modificar o teorema para podermos contemplar funções de
distribuições F em escada, como são as funções de distribuição de probabilidades associadas a essas
variáveis aleatórias.
Na Figura Figure 3.1 representamos como gerar uma realização de uma variável aleatória X com densidade
f e função de distribuição F . Assim, basta gerarmos um número uniforme u0 e invertermos a função de
distribuição F neste ponto. Computacionalmente a dificuldade é obtermos analiticamente uma expressão
para a função F −1 para muitos modelos probabilísticos. Em geral, essas expressões não existem e métodos
numéricos são requeridos para inverter a função de distribuição. Neste capítulo vamos apresentar este
método para a distribuição exponencial.
Figure 3.1: Ilustração do teorema fundamental da transformação de probabilidades para gerar uma variável
aleatória X com densidade f (x) = F ′ (x). A partir de um número aleatório uniforme u0 a função de
distribuição é invertida neste ponto para se obter x0 , com densidade f (x).
Outro método bastante geral que utilizaremos é denominado de método da amostragem por rejeição. Esse
método tem um forte apelo geométrico. Procuraremos, a princípio, descrever esse método de uma forma
bastante geral. Posteriormente, aplicaremos este método para gerarmos variáveis aleatórias de alguns
modelos probabilístico. A grande vantagem deste método contempla o fato de não precisarmos obter a
função de distribuição de probabilidade e nem a sua inversa. Estas estratégias só podem ser aplicadas
em muitos dos modelos probabilísticos existentes, se utilizarmos métodos numéricos iterativos. Seja f (x)
a função densidade de probabilidade para a qual queremos gerar uma amostra aleatória. A área sob a
curva para um intervalo qualquer de x corresponde à probabilidade de gerar um valor x nesse intervalo.
3.2. MÉTODOS GERAIS PARA GERAR REALIZAÇÕES DE VARIÁVEIS ALEATÓRIAS 47
Se pudéssemos gerar um ponto em duas dimensões, digamos (X,Y ), com distribuição uniforme sob a área,
então a coordenada X teria a distribuição desejada.
Para realizarmos de uma forma eficiente a geração de realizações variáveis aleatórias com densidade f (x),
evitando as complicações numéricas mencionadas anteriormente, poderíamos definir uma função qualquer
g(x). Essa função tem que ter algumas propriedades especiais para sua especificação. Deve possuir área
finita e ter para todos os valores x densidade g(x) superior a f (x). Essa função é denominada de função
de comparação. Outra característica importante que g(x) deve ter é possuir função de distribuição G(x)
analiticamente computável e invertível, ou seja, x = G−1 (u). Como a função g(x) não é necessariamente
R∞
uma densidade, vamos denominar a área sob essa curva no intervalo para x de interesse por A = −∞ g(x)dx.
Como G−1 é conhecida, podemos gerar pontos uniformes (x,y) que pertencem à área sob a curva g(x)
facilmente. Para isso basta gerarmos um valor de uma variável aleatória uniforme u1 entre 0 e A e
aplicarmos o teorema (Theorem 3.1). Assim, obtemos o primeiro valor do ponto (x0 ,y0 ) por x0 = G−1 (u1 ).
Para gerarmos a segunda coordenada do ponto não podemos gerar um valor de uma variável aleatória
uniforme no intervalo de 0 a A, sob pena de gerarmos um ponto que não está sob a curva g(x). Assim,
calculamos o valor de g no ponto x0 por g(x0 ). Geramos y0 = u2 , sendo u2 o valor de uma variável
aleatória uniforme entre 0 e g(x0 ). Assim, obtemos um ponto (x0 ,y0 ) uniforme sob a curva g(x). A
dificuldade deste método é justamente estabelecer essa função g(x) com as propriedades exigidas.
Vamos agora traçar as curvas correspondentes a g(x) e f (x) no mesmo gráfico. Se o ponto uniforme (x0 ,y0 )
está na área sob a curva f (x), ou seja se y0 ≤ f (x0 ), então aceitamos x0 como um valor válido de f (x); se por
outro lado o ponto estiver na região entre as densidades f (x) e g(x), ou seja se f (x0 ) < y0 ≤ g(x0 ), então
rejeitamos x0 . Uma forma alternativa de apresentarmos esse critério é tomarmos y0 de uma distribuição
U (0,1) e aceitarmos ou rejeitarmos x0 se y0 ≤ f (x0 )/g(x0 ) ou se y0 > f (x0 )/g(x0 ), respectivamente.
Ilustramos esse método na Figura Figure 3.2, sendo que a A representa a área total sob a curva g(x).
Figure 3.2: Método da rejeição para gerar um valor x0 da variável aleatória X com função densidade f (x)
que é menor do que g(x) para todo x. Nessa ilustração, x0 deve ser aceito.
Vamos ilustrar o primeiro método e salientar que o segundo método é o que ocorre na maioria dos casos
48 CHAPTER 3. VARIÁVEIS ALEATÓRIAS NÃO-UNIFORMES
− ln(1 − u)
x = F −1 (u) = (3.3)
λ
em que u é um número uniforme (0, 1).
Devido à distribuição uniforme ser simétrica, podemos substituir 1 − u na equação (Equation 3.3) por u.
Assim, para gerarmos uma realização de uma variável exponencial X, a partir de uma variável aleatória
uniforme, utilizamos o teorema Theorem 3.1 por intermédio da equação:
− ln(u)
x= . (3.4)
λ
O algoritmo Python para gerarmos realizações de variáveis aleatórias exponenciais é dado por:
# programa demonstrando a geração de n realizações de variáveis
# aleatórias exponenciais com parâmetro lamb, utilizamos
# a função [Link]() para
# gerarmos números aleatórios uniformes
import numpy as np
def rexpon(n, lamb = 1.0):
u = [Link](0.0, 1.0, n) # gera vetor u (U(0,1))
x = -[Link](u) / lamb # gera vetor x com distrib. exp.
return x # retorna o vetor x
# exemplo
rexpon(5, 0.1)
1 (x−µ)2
f (x) = √ e− 2σ 2 (3.5)
2πσ 2
Nenhuma outra função utilizando o teorema Theorem 3.1 será novamente apresentada, uma vez que
podemos facilmente adaptar as funções rexpon se tivermos um eficiente algoritmo de inversão da função
de distribuição do modelo probabilístico alvo. A dificuldade deste método é a necessidade de uma enorme
quantidade de cálculo para a maioria das densidades. Isso pode tornar ineficiente o algoritmo, pois o
tempo de processamento é elevado.
Podemos ainda aproveitar a relação entre algumas funções de distribuições para gerarmos realizações de
variáveis aleatórias de outras distribuições. Por exemplo, se X é normal com densidade (Equation 3.5),
N (µ,σ 2 ), podemos gerar realizações de Y = eX . Sabemos que fazendo tal transformação Y terá distribuição
log-normal, cuja densidade com parâmetros de locação α (µ) e escala β (σ) é:
1
ln(y)−α 2
− 21
f (y) = √ e β
, y > 0. (3.6)
yβ 2π
Um importante método usado para gerar dados da distribuição normal é o de Box-Müller, que é baseado na
generalização do método da transformação de variáveis para mais de uma dimensão. Para apresentarmos
esse método, vamos considerar p variáveis aleatórias X1 , X2 , . . . , Xp com função densidade conjunta
f (x1 ,x2 , . . . , xp ) e p variáveis Y1 , Y2 , . . ., Yp , funções de todos os X’s, então a função densidade conjunta
dos Y ’s é:
∂x1 ∂x1
∂y1 ... ∂yp
.. .. ..
f (y1 , . . . , yp ) = f (x1 , . . . , xp )abs . . . (3.7)
∂xp ∂xp
∂y1 ... ∂yp
√
= −2 ln x1 cos (2πx2 )
y1
√ (3.8)
y2 = −2 ln x1 sin (2πx2 )
= e− 2 (y1 +y2 )
1 2 2
x1
(3.9)
= 1
arctan y2
x
2 2π y1
Sabendo que a seguinte derivada dk arctan(g)/dx é dada por k(dg/dx)/(1 + g 2 ), em que g é uma função
de X, então o Jacobiano da transformação é:
2 2 2 2
∂x1 ∂x1 −y1 e−0,5(y1 +y2 ) −y2 e−0,5(y1 +y2 ) 1 −0,5(y12 +y22 )
∂y1 ∂y2
= − y2 y2 1 y2 =− e (3.10)
∂x2
∂y1
∂x2
∂y2 2πy12 1+ 22 2πy1 1+ 22 2π
y y
1 1
Assim, a função densidade conjunta de Y1 e Y2 é dada por f (x1 ,x2 )|J|, sendo
1 y2 1 − y22
1
f (y1 ,y2 ) = √ e− 2 √ e 2 . (3.11)
2π 2π
Desde que a densidade conjunta de Y1 e Y2 é o produto de duas normais independentes, podemos afirmar
que as duas variáveis geradas são normais padrão independentes, como pode ser visto em Johnson and
Wichern [1998] e Ferreira [2018].
Assim, podemos usar esse resultado para gerarmos variáveis aleatórias normais. A dificuldade, no entanto,
é apenas computacional. A utilização de funções trigonométricas como seno e cosseno pode limitar a
performance do algoritmo gerado tornando-o lento. Um truque apresentado em Press et al. [1992] é
bastante interessante para evitarmos diretamente o uso de funções trigonométricas. Esse truque representa
uma melhoria do algoritmo de Box-Müller e é devido a Marsaglia and Bray [1964].
Ao invés de considerarmos os valores das variáveis aleatórias uniformes x1 e x2 de um quadrado de lado
igual a 1 (quadrado unitário), tomarmos u1 e u2 como coordenadas de um ponto aleatório em um círculo
unitário (de raio igual a 1). A soma de seus quadrados R2 = U12 + U22 é uma variável aleatória uniforme
que pode ser usada como X1 . Já o ângulo que o ponto (u1 , u2 ) determina em relação ao eixo 1 pode
ser usado como um ângulo aleatório dado por Θ = 2πX2 . Podemos apontar que a vantagem da não
utilização direta da expressão (Equation √3.8) refere-se ao fato do√cosseno e do seno poderem ser obtidos
alternativamente por: cos (2πx2 ) = u1 / r2 e sin (2πx2 ) = u2 / r2 . Evitamos assim as chamadas de
funções trigonométricas. Na Figura Figure 3.3 ilustramos os conceitos apresentados e denominamos o
ângulo que o ponto (u1 ,u2 ) determina em relação ao eixo u1 por θ.
Agora podemos apresentar o função boxmuller para gerar dados de uma normal utilizando o algoritmo de
Box-Müller. Essa função utiliza a função polar para gerar dois valores aleatórios, de variáveis aleatórias
independentes normais padrão Y1 e Y2 . A função boxmuller é:
# Função boxmuller retorna uma amostra de tamanho n de
# uma distribuição normal com média mu e variância sigmaˆ2
# utilizando o método Polar Box-Müller
Figure 3.3: Círculo unitário mostrando um ponto aleatório (u1 , u2 ) com r2 = u21 + u22 representando x1 e
θ o ângulo que o ponto (u1 , u2 ) determina em relação ao eixo 1. No exemplo, o ponto está situado no
círculo unitário, conforme é exigido.
52 CHAPTER 3. VARIÁVEIS ALEATÓRIAS NÃO-UNIFORMES
n = 30000
x = boxmuller(n, 10, 2)
graf = [Link](x, bins=20, color='#14e8f3',rwidth=0.95,alpha=0.7)
[Link]()
4000
3000
2000
1000
0
2 4 6 8 10 12 14 16 18
Algumas aproximações são apresentadas em Atkinson and Pearce [1976] e serão apenas descritas na
sequência. Uma das aproximações faz uso da densidade Tukey-lambda e aproxima a normal igualando os
quatro primeiros momentos. Esse algoritmo, além de ser uma aproximação, tem a desvantagem de utilizar
a exponenciação que é uma operação lenta. Utilizando essa aproximação podemos obter uma variável
3.3. VARIÁVEIS ALEATÓRIAS DE ALGUMAS DISTRIBUIÇÕES IMPORTANTES 53
Outro método é baseado na soma de 12 ou mais variáveis uniformes (Ui ∼ U (0,1)) independentes. Assim,
P12
a variável X = i Ui − 6 tem distribuição aproximadamente normal com média 0 e variância 1. Isso
ocorre em decorrência do teorema do limite central e em razão de cada uma das 12 variáveis uniformes
possuírem média 1/2 e variância 1/12.
Estas duas aproximações tem valor apenas didático e não devem ser recomendadas como uma forma de
gerar variáveis normais. Muitas vezes esse fato é ignorado em problemas que requerem elevada precisão e
confiabilidade dos resultados obtidos. Quando isso acontece conclusões incorretas ou no mínimo imprecisas
podem ser obtidas.
De uma forma geral temos
k
X k
Ui −
i=1
2
X= r ∼ N (0,1),
k
12
pois E(Ui ) = (b − a)/2 e V (Ui ) = (b − a)2 /12, em que k é o número de uniformes que devem ser somadas
para cada realização da variável normal. A seguir, implementamos essas duas aproximações para fins de
treinamento em criação de funções com o Python.
A primeira é a Tukey-Lambda:
# Aproximação baseada na Tukey-Lambda
def tukeylambda(n, mu = 0, sigma = 1):
x = [Link](0.0,1.0,n)
x = (x**0.135 -(1-x)**0.135) / 0.1975
x = x * sigma + mu
return x
n = 1000000
x = tukeylambda(n, 10,10)
[Link](x)
[Link](x)
np.float64(10.00977105495767)
np.float64(9.992626739127573)
A segunda é a aproximação da soma de k uniformes, usando o teorema do limite central:
# Aproximação baseada na soma de k
# uniformes
def soma_un(u, k):
if k == 12:
rn = [Link](u) - 6
else:
rn = ([Link](u) - k / 2) / (k / 12)**0.5
return rn
u = [Link](0.0,1.0,k)
x. append(soma_un(u, k))
x = [Link](x) * sigma + mu
return x
n = 100000
k = 12
x = norm_tlc(n, 10, 10, k)
[Link](x)
[Link](x)
np.float64(10.010404136458673)
np.float64(10.02230976843949)
em que U1 , U2 ,· · ·, Un são variáveis uniformes (0,1) i.i.d. e I(•) é uma função indicadora. Então, X tem
distribuição binomial (n, p).
Lemma 3.1 (Soma de binomiais). Se X1 , X2 , · · ·, Xk são variáveis aleatórias binomiais independentes
Pk Pk
com (n1 , p), · · ·, (nk , p), então i=1 Xi tem distribuição binomial, com parâmetros ( i=1 ni , p).
3.4. DISTRIBUIÇÃO BINOMIAL 55
Lemma 3.2 (Tempo de Espera - Propriedade 1). Sejam G1 , G2 , · · · variáveis aleatórias geométricas
independentes e, X, o menor inteiro tal que
X+1
X
Gi > n.
i=1
A geométrica é a distribuição do tempo de espera até a ocorrência do primeiro sucesso no g-ésimo evento,
numa sequência de ensaios Bernoulli independentes. Assim, supõe-se que venham a ocorrer g − 1 fracassos,
cada um com probabilidade de ocorrência constante 1 − p, antes da ocorrência de um sucesso no g-ésimo
ensaio, com probabilidade p. Finalmente, o segundo lema do tempo de espera pode ser anunciado por:
Lemma 3.3 (Tempo de Espera - Propriedade 2). Sejam E1 , E2 , · · · variáveis aleatórias exponenciais
i.i.d., e X o menor inteiro tal que
X+1
X Ei
> −ln(1 − p).
i=1
n−i+1
A maneira mais simples de obtermos uma solução no caso discreto é realizarmos uma busca sequencial
a partir da origem. Para o caso binomial, este algoritmo da inversão, denominado de BINV, pode ser
implementado se utilizarmos a fórmula recursiva:
P (X = 0) = (1 − p)n
(3.14)
P (X = x) = P (X = x − 1) n−x+1 p
x 1−p
u = [Link](0.0,1.0,1)
k = k + 1
if u <= p:
x += 1
return(x)
x = []
for i in range(n):
[Link](func(size, prob))
return x
# Exemplo
prob = 0.5
size = 3
bu(size, prob)
bg(size, prob)
binv(size, prob)
n = 10000 # sample size
x = rbinom(n, size, prob, bg) # pode trocar bg por binv ou bu
graf = [Link](x,bins=size+1,color='#14e8f3',rwidth=0.95,alpha=0.7)
0
1
2
3500
3000
2500
2000
1500
1000
500
0
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Procedimentos de geração de números aleatórios poderiam ser apresentados para muitas outras distribuições
de probabilidades. Felizmente no python não temos este tipo de preocupação, pois estas rotinas já existem
e estão implementadas em linguagem não interpretada. Inclusive para os modelos considerados temos
rotinas prontas. Veremos uma boa parte delas na próxima seção.
beta = 2.0
x = [Link](alpha, beta, n)
graf = [Link](x, bins='auto', color='#14e8f3',rwidth=0.95,alpha=0.7)
[Link](x)
alpha/(alpha+beta) # verdadeira
[Link](x)
alpha*beta/((alpha+beta)**2*(alpha+beta+1)) # verdadeira
np.float64(0.3349115920603343)
0.3333333333333333
np.float64(0.05559456122798455)
0.05555555555555555
2000
1750
1500
1250
1000
750
500
250
0
0.0 0.2 0.4 0.6 0.8 1.0
3.6 Exercícios
1. Seja f (x) = 3x2 uma função densidade de uma variável aleatória contínua X com domínio definido
no intervalo [0; 1]. Aplicar o método da inversão e descrever um algoritmo para gerar realizações de
variáveis aleatórias dessa densidade. Implementar em R e gerar uma amostra de tamanho n = 1.000.
Estimar os quantis 1%, 5%, 10%, 50%, 90%, 95% e 99%. Confrontar com os quantis teóricos.
2. Os dados a seguir referem-se ao tempo de vida, em dias, de n = 40 insetos. Considerando que a
distribuição do tempo de vida é a exponencial e quePo parâmetro λ pode ser estimado pelo estimador
n
de máxima verossimilhança λ̂ = 1/X̄, em que X̄ = i=1 Xi /n, obter o intervalo de 95% de confiança
utilizando o seguinte procedimento: i) gerar uma amostra da exponencial de tamanho n = 40,
utilizando o algoritmo rexpon, considerando o parâmetro igual a estimativa obtida; ii) determinar
a estimativa da média µ = 1/λ por X̄ nesta amostra simulada de tamanho n = 40; iii) repetir
1.000 vezes os passos (i) e (ii) e armazenar os valores obtidos; iv) ordenar as estimativas e tomar os
quantis 2,5% e 97,5%. Os valores obtidos são o intervalo de confiança almejado, considerando como
verdadeira a densidade exponencial para modelar o tempo de vida dos insetos. Este procedimento é
denominado de bootstrap paramétrico. Os dados em dias do tempo de vida dos insetos são:
3.6. EXERCÍCIOS 61
Repetir esse processo, gerando 100.000 amostras de tamanho n = 40. Compare os resultados e
verifique se o custo adicional de ter aumentado o número de simulações compensou a possível maior
precisão obtida.
3. Gerar uma amostra de n = 5.000 realizações de variáveis normais padrão utilizando as aproximações:
X = [U 0,135 −(1 − U )0,135 ]/ 0,1975 e da soma de 12 ou mais variáveis uniformes (Ui ∼ U (0,1))
P12
independentes, dada por X = i Ui − 6. Confrontar os quantis 1%, 5%, 10%, 50%, 90%, 95% e
99% esperados da distribuição normal com os estimados dessa distribuição. Gerar também uma
amostra de mesmo tamanho utilizando o algoritmo Polar-Box-Müller. Estimar os mesmos quantis
anteriores nesta amostra e comparar com os resultados anteriores.
4. Se os dados do exercício 2 pudessem ser atribuídos a uma amostra aleatória da distribuição log-
normal, então estimar os parâmetros da log-normal e utilizar o mesmo procedimento descrito naquele
exercício, substituindo apenas a distribuição exponencial pela log-normal para estimar por intervalo
a média populacional. Para estimar os parâmetros da log-normal utilizar o seguinte procedimento:
a) transformar os dados originais, utilizando Xi∗ = ln(Xi ); b) determinar a média e o desvio padrão
amostral dos dados transformados - estas estimativas são as estimativas de µ e σ. Utilizar estas
estimativas para gerar amostras log-lognormais. Realizar os mesmos procedimentos descritos para
exponencial, confrontar os resultados e discutir a respeito da dificuldade de se tomar uma decisão
da escolha da distribuição populacional no processo de inferência. Como em situações reais nunca se
sabe de qual distribuição os dados são provenientes com precisão, então você teria alguma ideia de
como fazer para determinar qual a distribuição que melhor modela os dados do tempo de vida dos
insetos? Justificar sua resposta adequadamente com os procedimentos numéricos escolhidos.
5. Fazer reamostragens com reposição a partir da amostra do exercício 2 e estimar o intervalo de 95%
para a média populacional, seguindo os passos descritos a seguir e utilizar um gerador de números
uniformes para determinar quais elementos amostrais devem ser selecionados: i) reamostrar com
reposição os n = 40 elementos da amostraPnoriginal e compor uma nova amostra por: Xi∗ ; ii) calcular
a média desta nova amostra por X̄ = i=1 Xi /n; iii) armazenar este valor e repetir os passos (i)
∗ ∗
e (ii) B − 1 vezes; iv) agrupar os valores com o valor da amostra original; e v) ordenar os valores
obtidos e determinar os quantis 2,5% e 97,5% desse conjunto de B valores. Escolher B = 1.000
e B = 100.000 e confrontar os resultados obtidos com os obtidos nos exercícios anteriores para a
distribuição exponencial e log-normal. Os resultados que estiverem mais próximos deste resultado
devem fornecer um indicativo da escolha da distribuição mais apropriada para modelar o tempo de
vida de insetos. Este procedimento sugerido neste exercício é o bootstrap não-paramétrico. A sua
grande vantagem é não precisar fazer suposição a respeito da distribuição dos dados amostrais.
6. Uma importante relação para obtermos intervalos de confiança para a média de uma distribuição
exponencial, f (x) = λe−λx , refere-se ao fato de que a soma de n variáveis exponenciais com parâmetro
y
α−1 − β
λ é igual a uma gama f (y) = λΓ(α)1
(y/β) e com parâmetros α = n e β = 1/λ. Assim, assumir
que os dados do exercício 2 têm distribuição exponencial com parâmetro λ estimado pelo recíproco
da média amostral, ou seja, β̂ = 1/X̄. Considerando que a variável X tem distribuição gama padrão
com parâmetro α = n, então obtenha Y = β̂X = X/λ̂. Neste caso Y |β̂ tem distribuição da soma
de n variáveis exponenciais, ou seja, distribuição gama com parâmetros α = n e β = β̂. Como
queremos a distribuição da média, devemos obter a transformação Ȳ = Y /n. Gerar amostras de
tamanho n = 1.000 e n = 100.000 e estimar os quantis 2,5% e 97,5% da distribuição de Ȳ , em cada
62 CHAPTER 3. VARIÁVEIS ALEATÓRIAS NÃO-UNIFORMES
uma delas. Confrontar os intervalos de confiança, dessa forma obtidos, com os do exercício 2. Quais
são as suas conclusões? Qual é a vantagem de utilizar a distribuição gama?
7. Duas amostras binomiais foram realizadas em duas (1 e 2) diferentes populações. Os resultados
do número de sucesso foram y1 = 2 e y2 = 3 em amostras de tamanho n1 = 12 e n2 = 14,
respectivamente, de ambas as populações. Estimar os parâmetros p1 e p2 das duas populações por:
pˆ1 = y1 /n1 e pˆ2 = y2 /n2 . Para testarmos a hipótese H0 : p1 = p2 ou H0 : p1 − p2 = 0, podemos
utilizar o seguinte algoritmo bootstrap paramétrico: a) utilizar pˆ1 e pˆ2 para gerarmos amostras
de tamanho n1 e n2 de ambas as populações; b) estimar pˆ1j = y1j /n1 e pˆ2j = y2j /n2 na j-ésima
repetição desse processo; c) calcular dj = pˆ1j − pˆ2j ; d) repetir os passos de (a) a (c) B − 1 vezes;
e) unir com o valor da amostra original; f) ordenar os valores e obter os quantis 2,5% e 97,5% da
distribuição bootstrap de dj ; e g) se o valor hipotético 0 estiver contido nesse intervalo, não rejeitar
H0 , caso contrário, rejeitar a hipótese de igualdade das proporções binomiais das duas populações.
Chapter 4
Os modelos multivariados ganharam grande aceitação no meio científico em função das facilidades
computacionais e do desenvolvimento de programas especializados nesta área. Os fenômenos naturais
são em geral multivariados. Um tratamento aplicado em um ser, a um solo ou a um sistema não afeta
isoladamente apenas uma variável, e sim todas as variáveis. Ademais, as variáveis possuem relações entre si
e qualquer mudança em uma ou algumas delas, afeta as outras. Assim, a geração de realizações de vetores
ou matrizes aleatórias é um assunto que não pode ser ignorado. Vamos neste capítulo disponibilizar ao
leitor mecanismos para gerar realizações de variáveis aleatórias multidimensionais.
4.1 Introdução
Os processos para gerarmos variáveis aleatórias multidimensionais são muitas vezes considerados difíceis
pela maioria dos pesquisadores. Uma boa parte deles, no entanto, pode ser realizada no Python com
apenas uma linha de comando. Embora tenhamos estas facilidades, nestas notas vamos apresentar detalhes
de alguns processos para gerarmos dados dos principais modelos probabilísticos multivariados como, por
exemplo, a normal multivariada, a Wishart e a Wishart invertida, a t de Student multivariada e algumas
outras distribuições.
Uma das principais características das variáveis multidimensionais é a correlação entre seus componentes.
A importância destes modelos é praticamente indescritível, mas podemos destacar a inferência paramétrica,
a inferência bayesiana, a estimação de regiões de confiança, entre outras. Vamos abordar nas próximos
seções formas de gerarmos realizações de variáveis aleatórias multidimensionais para determinados modelos
utilizando o Python para implementarmos as rotinas ou para utilizarmos as rotinas pré-existentes. Nossa
aparente perda de tempo, descrevendo funções menos eficientes do que as pré-existentes no Python,
tem como razão fundamental permitir ao leitor ir além de simplesmente utilizar rotinas previamente
programadas por terceiros. Se o leitor ganhar, ao final da leitura deste material, a capacidade de produzir
suas próprias rotinas e entender como as rotinas pré-existentes funcionam, nosso objetivo terá sido
alcançado.
63
64CHAPTER 4. GERAÇÃO DE AMOSTRAS ALEATÓRIAS DE VARIÁVEIS MULTIDIMENSIONAIS
1
−p − 12
fX (x) = (2π) 2 |Σ| exp − (x − µ)⊤ Σ−1 (x − µ) (4.1)
2
X = C−1 Y.
1 −1
− 21
fY (y) =(2π) −p/2
|Σ| exp − (C y − µ) Σ (C y − µ) |C|−1
⊤ −1 −1
2
−1
=(2π)−p/2 |C|−1/2 |Σ| 2 |C|−1/2
1
⊤ −1 −1 −1
× exp − (y − Cµ) C ⊤
Σ C (y − Cµ)
2
1
1
⊤ −2
−1
=(2π) −p/2
CΣC exp − (y − Cµ)⊤ CΣC⊤ (y − Cµ)
2
O teorema Theorem 4.1 nos fornece o principal resultado para gerarmos dados de uma normal multivariada.
Assim, como nosso objetivo é gerar dados de uma amostra normal multivariada com vetor de médias µ
e matriz de covariâncias Σ pré-estabelecidos, devemos seguir os seguintes procedimentos. Inicialmente
devemos obter a matriz raiz quadrada de Σ, representada por Σ1/2 . Para isso, vamos considerar a
decomposição espectral da matriz de covariâncias dada por Σ = PΛP⊤ . Logo, podemos definir a matriz
raiz quadrada de Σ por Σ1/2 = PΛ1/2 P⊤ , em que Λ é a matriz diagonal dos autovalores, Λ1/2 é a matriz
diagonal contendo a raiz quadrada destes elementos e P é a matriz de autovetores, cada um destes vetor
disposto em uma de suas colunas. Desta decomposição facilmente podemos observar que Σ = Σ1/2 Σ1/2 .
Assim, podemos utilizar o seguinte método para gerarmos uma realização p-variada de uma normal
multivariada. Inicialmente devemos gerar um vetor aleatório Z = [Z1 , Z2 , · · ·, Zp ]⊤ de p variáveis normais
padrão independentes utilizando, por exemplo, o algoritmo de Box-Müller. Isto que dizer que µZ = 0 e
que Cov(Z) = I. Este vetor deve sofrer a seguinte transformação linear:
Y = Σ1/2 Z + µ. (4.2)
4.2. DISTRIBUIÇÃO NORMAL MULTIVARIADA 65
De acordo com o teorema Theorem 4.1, o vetor Y possui distribuição normal multivariada com média
µY = Σ1/2 µZ +µ = µ e matriz de covariâncias Σ1/2 IΣ1/2 = Σ. Este será o método que usaremos para obter
a amostra p-dimensional de tamanho n de uma normal multivariada com média µ e covariância Σ. Para
obtermos a matriz raiz quadrada no Python, podemos utilizar o comando para a obtenção da decomposição
espectral [Link]() ou alternativamente o comando [Link](a,/,*,upper=False), que
retorna o fator de Cholesky de uma matriz positiva definida, que na verdade é um tipo de raiz quadrada.
A decomposição do valor singular neste caso se especializa na decomposição espectral, pois a matriz Σ
é simétrica. Geraremos um vetor de variáveis aleatórias normal padrão independentes Z utilizando o
comando [Link](). Em seguida a transformação Equation 4.2 é realizada.
Vamos ilustrar e apresentar o programa de geração de variáveis normais multivariada para um caso
particular bivariado (p = 2) e com vetor de médias µ = [10, 50]⊤ e matriz de covariâncias dada por:
4 1
Σ= .
1 1
# Exemplo de uso
n = 10000
sigma = [Link]([[4,1.9],[1.9,2]])
mu = [Link](2, 0) # cria o vetor com 2 valores 0
x = rnormmv(n, mu, sigma)
[Link](x, axis = 0)
[Link](x,rowvar = False)
a = x[:,0]
b = x[:,1]
[Link](a, b, s = 3, c = '#070808', alpha = 0.5)
[Link]('x1')
[Link]('x2')
[Link]()
array([-0.00108366, 0.00969471])
array([[4.01124863, 1.91681719],
[1.91681719, 2.01578024]])
66CHAPTER 4. GERAÇÃO DE AMOSTRAS ALEATÓRIAS DE VARIÁVEIS MULTIDIMENSIONAIS
Text(0.5, 0, 'x1')
6
4
2
x2
0
2
4
6
7.5 5.0 2.5 0.0 2.5 5.0 7.5
x1
No python temos o gerador do numpy, random.multivariate_normal(mean, cov, size=None,
check_valid='warn', tol=1e-8) para gerarmos dados da distribuição normal multivariada. No script
a seguir aplicamos a nossa função e a do numpy para gerarmos um número grande de observações e
comparamos o desempenho em termos de tempo de processamento. A função do numpy protege o processo,
com, por exemplo verificando se a matriz de covariâncias é positiva definida. Em nossa função não
utilizamos nenhuma proteção, embora seja possível fazer isso.
# Comparativo de desempenho dos dois geradores de
# normais multivariadas
from datetime import datetime
n = 60000
sigma = [Link]([[4,1.9],[1.9,2]])
mu = [Link](2, 0) # cria o vetor com 2 valores 0
t1 = [Link]()
x = rnormmv(n, mu, sigma)
t2 = [Link]()
trnm = t2-t1
print('tempo médio da rnormmv: ',[Link] / n)
t1 = [Link]()
x = [Link].multivariate_normal(mu, sigma, size=n)
t2 = [Link]()
tnp = t2-t1
print('tempo médio da numpy: ',[Link] / n)
# tempo relativo
print('numpy é mais rápido ',[Link] / [Link] ,' vezes')
A rotina do numpy foi muitas vezes superior a nossa implementação. Isso em parte é devido ao fato de
estar compilada e, talvez, também ao possível método utilizado para obter a matriz raiz quadrada. A
troca da matriz raiz quadrada de svd para Cholesky ou pela decomposição espectral, poderia ser feita e
o desempenho médio dos procedimentos avaliados, para definirmos a melhor estratégia de implementação.
Fizemos uma alteração na nossa implementação. Colocamos um argumento com uma das três opções,
para que o usuário escolha qual método utilizar. A obtenção da raiz quadrada foi implementada em uma
função separada.
# Função Python para gerar n vetores aleatórios normais
# multivariados com vetor de médias mu e covariância sigma.
# o resultado é uma matriz n x p, sendo p o número de variáveis
def sigmaroot(sigma, metodo='chol'):
if metodo == 'svd':
u, d, vt = [Link](sigma)
root = u @ [Link](d**0.5) @ vt
elif metodo == 'eig':
d, u = [Link](sigma)
root = u @ [Link](d**0.5) @ [Link](u)
else:
root = [Link](sigma)
return root
t1 = [Link]()
x = [Link].multivariate_normal(mu, sigma, size=n)
t2 = [Link]()
tnp = t2-t1
print('tempo médio da svd: ',[Link] / n)
print('tempo médio da eig: ',[Link] / n)
print('tempo médio da chol: ',[Link] / n)
print('tempo médio da numpy: ',[Link] / n)
# tempo relativo
print('numpy é mais rápido ',tsvd/tnp ,' vezes que svd')
print('numpy é mais rápido ',teig/tnp ,' vezes que eig')
print('numpy é mais rápido ',tchol/tnp,' vezes que chol')
Embora nossa função seja relativamente rápida, levando entre 2 e 14 micro segundos para rodar cada
observação bivariada neste caso, ela ainda foi bem menos eficiente que rotina numpy. A compilação é
fundamental em relação à interpretação. É extremamente simples gerarmos dados de normais multivariadas
utilizando a função numpy, o que nos desobriga de programar os nossos próprios geradores aleatórios.
Reiteramos que fizemos isso, pois queremos que o nosso leitor e nosso estudante consigam desvendar o que
está por trás de cada método deste e também que consiga desenvolver suas habilidades em programação,
implementando rotinas sofisticadas como estas.
n−1
X
W= Xj X⊤
j
j=1
possui distribuição Wishart com n − 1 graus de liberdade e parâmetro Σ (matriz positiva definida).
Da mesma forma, se temos uma amostra aleatória de tamanho n de uma distribuição normal multivariada
com média µ e covariância Σ, a distribuição da matriz aleatória
n
X
W= (Xj − X̄)(Xj − X̄)⊤
j=1
A função densidade Wishart de uma matriz aleatória W de somas de quadrados e produtos e representada
por Wp (ν,Σ) é definida por:
4.3. DISTRIBUIÇÃO WISHART E WISHART INVERTIDA 69
( )
|Σ|−ν/2 |w|(ν−p−1)/2 tr Σ−1 w
fW (w|ν,Σ) = exp − (4.3)
p
+ 1 2
Y ν − i
2νp/2 π p(p−1)/4 Γ
i=1
2
R∞
em que Γ(x) = 0
tx−1 e−x dt é função gama.
Assim, para gerarmos variáveis Wishart com parâmetros n − 1 graus liberdade inteiro e matriz Σ positiva
definida, podemos utilizar um gerador de amostras aleatórias normais multivariadas e obter a matriz de
somas de quadrados e produtos amostrais. Esta matriz será uma realização de uma variável aleatória
Wishart, que é uma matriz de dimensão p × p. A seguinte função pode ser utilizada para obtermos
realizações aleatórias de uma Wishart:
# Exemplificação para gerarmos matrizes de somas de quadrados e produtos
# aleatórias W com distribuição Wishart(nu, Sigma), nu = n - 1
# utiliza o random.multivariate_normal para gerar normais multivariadas.
def rwishart(nu, sigma):
p = [Link][0]
mu = [Link](p, 0)
x = [Link].multivariate_normal(mu, sigma, size=nu + 1)
w = nu * [Link](x, rowvar=False)
return w
# Exemplo de uso
sigma = [Link]([[4, 1], [1, 1]])
nu = 5
w = rwishart(nu, sigma)
print(w)
[[32.07270987 11.02963336]
[11.02963336 7.00507668]]
Outra distribuição relacionada que aparece frequentemente na inferência multivariada é a Wishart invertida.
Considere W uma matriz aleatória Wp (ν, Σ), então a distribuição de S = W−1 , dada pela função densidade
( )
|Σ−1 |ν/2 |s|−(ν+p+1)/2 tr Σ−1 s−1
fS (s|ν,Σ) = exp − (4.4)
p
ν−i+1 2
Y
2νp/2 π p(p−1)/4 Γ
i=1
2
# Exemplo de uso
sigma = [Link]([[4, 1], [1, 1]])
nu = 5
wi = True
res = rw_wi(nu, sigma, wi)
print('Ambas, W e WI', res)
wi = False
res = rw_wi(nu, sigma, wi)
print('Só a Wishart: ', res)
s = [Link](sigma)
else:
t = [Link](df)**0.5
s = sigma**0.5
if tip == int or tip == float:
w = s**2 * t**2
else:
w = s @ t @ [Link](t) @ [Link](s)
if wi == True:
if tip == int or tip == float:
iw = 1 / w
else:
iw = [Link](w)
res = {'W': w, 'WI': iw}
return res
else:
return w
# Exemplo de uso
sigma = [Link]([[4, 1], [1, 1]])
nu = 5
wi = True
res = rwwi_sh(nu, sigma, wi)
print('Ambas, W e WI', res)
wi = False
res = rwwi_sh(nu, sigma, wi)
print('Só a Wishart: ', res)
Vamos chamar a atenção para alguns fatos sobre esta função, pois utilizamos alguns coman-
dos e recursos não mencionados até o presente momento. Inicialmente utilizamos o comando
[Link]([Link](df, p)**0.5) para preencher a diagonal da matriz T com realizações
de variáveis aleatórias qui-quadrado e em seguida do método diag para transformar o vetor em uma
matriz diagonal.
Outro aspecto interessante que merece ser mencionado é o uso do fator de Cholesky no lugar de obter
a matriz raiz quadrada de Σ. O fator de Cholesky utiliza a decomposição Σ = SS⊤ , em que S é uma
matriz triangular inferior. O Python, por meio da função [Link](sigma), retorna a matriz
S. Finalmente, se o número de variáveis é igual a 1, a distribuição Wishart se especializa na qui-quadrado
e a Wishart invertida na qui-quadrado invertida. Assim, quando p = 1 o algoritmo retornará variáveis
σ 2 X e 1/(σ 2 X) com distribuições proporcionais a distribuição qui-quadrado e qui-quadrado invertida,
respectivamente, sendo X uma variável qui-quadrado com ν graus de liberdade.
O processamento para o caso escalar foi controlado com o uso do if, pois operações matriciais não são
aplicáveis se os argumentos forem escalares inteiros ou float (reais).
72CHAPTER 4. GERAÇÃO DE AMOSTRAS ALEATÓRIAS DE VARIÁVEIS MULTIDIMENSIONAIS
Γ ν+p
z −(ν+p)/2
g(z) = 2
1+ .
(πν) Γ(ν/2)
p/2 ν
− ν+p
Γ ν+p 1 ⊤
2
fY (y) = 2
1+ y y . (4.5)
(πν) Γ(ν/2)
p/2 ν
A variável aleatória Y terá vetor de médias nulo e covariâncias νI/(ν − 2) e a densidade terá contornos
esféricos de mesma probabilidade.
Vamos apresentar a forma geral para gerarmos variáveis aleatórias p-dimensionais t multivariada com ν
graus de liberdade e parâmetros µ e Σ. Seja um vetor aleatório Z com distribuição Np (0, I) e a variável
aleatória U com distribuição qui-quadrado com ν graus de liberdade, então o vetor aleatório Y, dado pela
transformação
√ Z
Y = ν√ , (4.6)
U
possui distribuição t multivariada esférica com ν graus de liberdade. O vetor X obtido pela transformação
linear
X = Σ1/2 Y + µ, (4.7)
possui distribuição t multivariada elíptica com ν graus de liberdade e parâmetros µ e Σ.
Assim, devemos aplicar a transformação (Equation 4.7) n vezes a n diferentes vetores aleatórios Y e
variáveis U . Ao final deste processo teremos uma amostra de tamanho n da distribuição t multivariada
almejada com ν graus de liberdade. Assim, para gerarmos dados de uma t multivariada com dimensão p,
graus de liberdade ν (não necessariamente inteiro), vetor de média µ qualquer e matriz positiva definida
Σ podemos utilizar a seguinte função Python, substituindo na expressão (Equation 4.7) a matriz raiz
quadrada pelo fator de Cholesky F de Σ:
# Função para gerarmos variáveis aleatórias t
# multivariadas (n, mu, Sigma, nu).
def rtmult(n, mu=[0,0,0], sigma=[Link](3), df = 1):
4.4. DISTRIBUIÇÃO T DE STUDENT MULTIVARIADA 73
f = [Link](sigma)
p = [Link][0]
x = [Link].multivariate_normal([Link](p),[Link](p),n)
q = ([Link](df,n) / df)**0.5
x = x / q[:, [Link]] @ [Link](f)
x += [Link](mu, n).reshape(n,p)
return x
# Exemplo de uso
n = 3000
nu = 3
mu = [0,0]
sigma = [Link](2)
x = rtmult(n, mu, sigma, nu)
fig, (ax1, ax2) = [Link](1, 2)
[Link]('Distribuições t multivariadas')
[Link](x[:,0], x[:,1], s = 3, c = '#070808', alpha = 0.5)
ax1.set_xlabel('x1')
ax1.set_ylabel('x2')
mu = [10, 5]
sigma = [Link]([[4, 1.9], [1.9,1]])
x = rtmult(n, mu, sigma, nu)
[Link](x[:,0], x[:,1], s = 3, c = '#070808', alpha = 0.5)
ax2.set_xlabel('x1')
ax2.set_ylabel('x2')
[Link](x, rowvar = False)
nu * sigma
Text(0.5, 0, 'x1')
Text(0.5, 0, 'x1')
array([[12.14949354, 5.67204001],
[ 5.67204001, 2.91537628]])
array([[12. , 5.7],
[ 5.7, 3. ]])
74CHAPTER 4. GERAÇÃO DE AMOSTRAS ALEATÓRIAS DE VARIÁVEIS MULTIDIMENSIONAIS
Distribuições t multivariadas
20
15 15
10 10
5 5
0 0
x2
x2
5
5
10
10
15
15 20
10 0 10 20 50 25 0 25
x1 x1
Graus de liberdade reais positivos podem ser utilizados como argumento da função criada. Foram geradas,
no exemplo anterior, duas amostras aleatórias para ilustrar. Foram definidos para os parâmetros µ, Σ e ν,
os valores [0,, 0, 0]⊤ , I3 e 3, respectivamente, como default. Podemos também utilizar a implementação
determinada pela função multivariate_t.rvs(mu, sigma, df n), da biblioteca scipy, para gerarmos
dados da distribuição t de Student multivariada.
from [Link] import multivariate_t
n = 5
X = multivariate_t.rvs([1.0, -0.5], [[2.1, 0.3], [0.3, 1.5]], df=3, size = n)
print(X)
[[ 0.31880772 0.00700063]
[ 1.66875877 -0.62235628]
[ 2.57363526 0.56635033]
[ 0.40409292 -0.90737645]
[ 2.59648256 -0.59989829]]
# Exemplo de uso
n = 7
print(rlgnormalmv(n)) # lgnormalmv padrão
(x − µ1 )⊤ Σ−1
1 (x − µ1 )
fX (x) =δ(2π) −p/2
|Σ1 |−1/2
exp −
2
(x − µ2 )⊤ Σ−1
2 (x − µ2 )
+ (1 − δ)(2π)−p/2 |Σ2 |−1/2 exp −
2
# Exemplo
n = 2000
delta = 0.8
mu1 = [Link](2)
sig1 = [Link](2)
mu2 = [10, 5]
sig2 = [Link]([[2, 1.4], [1.4, 1]])
x = rncm(n, mu1, mu2, sig1, sig2, delta)
[Link](x[:,0], x[:,1], s = 3, c = '#070808', alpha = 0.5)
[Link]()
76CHAPTER 4. GERAÇÃO DE AMOSTRAS ALEATÓRIAS DE VARIÁVEIS MULTIDIMENSIONAIS
8
6
4
2
0
2
4
2.5 0.0 2.5 5.0 7.5 10.0 12.5 15.0
4.6 Exercícios
1. Gerar uma amostra de tamanho (n = 10) uma distribuição normal trivariada com vetor de médias
µ = [5,10, 15]⊤ e matriz de covariâncias
5 −1 2
Σ = −1 3 −2 .
2 −2 7
Estimar a média e a covariância amostral. Repetir este processo 1.000 vezes e estimar a média das médias
amostrais e a média das matrizes de covariâncias amostrais. Estes valores correspondem exatamente aos
respectivos valores esperados? Se não, apresentar a(s) principal(is) causa(s).
2. A partir da alteração da função rnormmv que realizamos, comparar o tempo de processamento
médio quando for utilizado svd, eig ou cholesky que corresponde a um possível tipo de matriz
raiz quadrada de Σ. Verificar se houve melhoria no desempenho em relação ao tempo médio de
processamento de cada variável aleatória gerada, vetor p-dimensional. Testar isso usando modificando
n e p, pois só apresentamos o teste para um único valor de p e pequeno, p = 2.
3. Sabemos que variáveis Wishart possuem média νΣ. Apresentar um programa para verificar se o
valor esperado, ignorando o erro de Monte Carlo, é alcançado com o uso das funções apresentadas
para gerar variáveis Wishart. Utilizar 10.000 repetições de Monte Carlo. Isso seria uma forma
simples, embora não conclusiva, de checar se a função está realizando as simulações corretamente.
Serve ao menos para indicar a presença de erro, mas não garante a assertividade.
4. Implementar funções Python para gerarmos variáveis aleatórias log-normal e normal-contaminada
elípticas multivariadas.
Chapter 5
Muitos algoritmos para o cálculo de médias, variâncias e covariâncias são imprecisos, podendo gerar
resultados finais contendo grandes erros. A necessidade de utilizarmos algoritmos eficientes para realizarmos
estas tarefas simples são evidentes e serão descritos neste capítulo.
5.1 Introdução
Felizmente o Python utiliza algoritmos precisos para cálculo da média, da variância e de covariância.
Vamos buscar esclarecer como a utilização de algoritmo ineficientes podem levar a resultados inconsistentes
e imprecisos. Nosso objetivo neste capítulo é apresentar os algoritmos eficientes para estimarmos estes
parâmetros quando as fórmulas convencionais podem falhar se utilizadas nos algoritmos diretamente.
Estes algoritmos eficientes são particularmente úteis quando os dados possuem grande magnitude ou
estão muito próximos de zero. Neste caso particular, algumas planilhas eletrônicas, como o Excel nas
suas versões mais antigas, podiam falhar [McCullough and Wilson, 1999]. O conhecimento de algoritmos
que conduzirão a maiores precisões numéricas pode levar o pesquisador a não cometer os mesmos erros
encontrados em alguns softwares.
n
X
Xi
i=1
X̄. = , (5.1)
n
77
78 CHAPTER 5. ALGORITMOS PARA MÉDIAS, VARIÂNCIAS E COVARIÂNCIAS
e
n
!2
X
n Xi
1 X
i=1
S2 = X2 − . (5.2)
n − 1 i=1 i n
Alguns algoritmos existentes procuraram melhorar os algoritmos dos livros textos que são as fórmulas das
equações (Equation 5.1 e Equation 5.2) procurando fazer adaptações, repassando a amostra duas vezes.
Estes algoritmos podem ser eficientes do ponto de vista da precisão, mas não são rápidos, justamente por
repassarem duas vezes os dados. West [1979] propôs utilizar um algoritmo que faz uma única passada,
atualizando a média e a variância em cada nova observação. Podemos facilmente mostrar para uma
amostra de tamanho n, que a média é igual a X1 se n = 1, (X1 + X2 )/2 se n = 2 e assim por diante. No
(k − 1)-ésimo passo podemos especificar o estimador da média por:
k−1
X
Xi
i=1
X̄k−1 = .
k−1
k
X
Xi
i=1
X̄k = . (5.3)
k
A pergunta que fazemos é “podemos expressar a média do k-ésimo passo em função da média do (k − 1)-
ésimo passo?” A resposta a esta pergunta é sim e o resultado nos fornece o algoritmo desejado. A partir
da equação (Equation 5.3) obtemos:
k−1
X
Xi + Xk
i=1
X̄k =
k
k−1
X
(k − 1) Xi
i=1 Xk
= +
(k − 1)k k
(k − 1)X̄k−1 Xk
= +
k k
X̄k−1 Xk
=X̄k−1 − +
k k
resultando na equação recursiva final
Xk − X̄k−1
X̄k = X̄k−1 + , (5.4)
k
1 < k ≤ n por:
k
!2
X
k
Xi
X i=1
W 2k = Xi2 −
i=1
k
veremos que a variância correspondente é dada por Sk2 = W 2k /(k − 1). Se expandirmos esta expressão
isolando o k-ésimo termo e simplificarmos a expressão resultante teremos:
k−1
!2
X
k−1
Xi + Xk
(5.5)
X i=1
W 2k = Xi2 + Xk2 −
i=1
k
2
= W 2k−1 + k Xk − X̄k /(k − 1)
A expressão que desenvolvemos (Equation 5.5) é equivalente a apresentada por West [1979]. Isso pode ser
demonstrado facilmente se substituirmos X̄k obtida na equação (Equation 5.4) na equação (Equation 5.5),
de onde obtivemos:
2
W 2k = W 2k−1 + (k − 1) Xk − X̄k−1 /k (5.6)
para 2 ≤ k ≤ n, sendo que W 21 = 0. A variância é obtida por S = W 2n /(n − 1).2
i=1
k
k−1
!2 k−1
X X
k−1
X i + 2X k Xi + Xk2
X
= Xi2 + Xk2 − i=1 i=1
i=1
k
k−1
!2 k−1
X X
k−1
(k − 1) X i + 2(k − 1)X k Xi + (k − 1)Xk2
X i=1 i=1
= Xi2 + Xk2 −
i=1
k(k − 1)
k−1
!2 k−1
!2 k−1
X X X
k−1
X i X i 2(k − 1)X k Xi
X i=1 i=1 i=1 X2
= 2
Xi − + Xk +2
− − k
i=1
k−1 k(k − 1) k(k − 1) k
k−1
! 2 k−1
X X
Xi 2(k − 1)Xk Xi
i=1 i=1 (k − 1)Xk2
=W 2k−1 + − +
k(k − 1) k(k − 1) k
(k − 1)X̄k−1
2
2(k − 1)Xk X̄k−1 (k − 1)Xk2
=W 2k−1 + − + ,
k k k
resultando em
2
W 2k = W 2k−1 + (k − 1) Xk − X̄k−1 /k.
80 CHAPTER 5. ALGORITMOS PARA MÉDIAS, VARIÂNCIAS E COVARIÂNCIAS
Podemos generalizar essa expressão para computarmos a covariância S(x,y) entre uma variável X e outra
Y . A expressão para a soma de produtos é dada por:
media : 4.285714285714286
variância : 6.571428571428572
SQ2 : 39.42857142857143
W3 : 22.04081632653064
W4 : 391.4518950437319
Xk − X̄k−1
X̄k = X̄k−1 + , (5.10)
k
para 2 ≤ k ≤ n, sendo que X̄1 = X1 .
Da mesma forma, a adaptação para p dimensões da expressão (Equation 5.6) é direta e o resultado obtido
é:
⊤
Wk = Wk−1 + (k − 1) Xk − X̄k−1 (5.11)
Xk − X̄k−1 /k
para 2 ≤ k ≤ n, sendo W1 = 0, uma matriz de zeros de dimensões (p × p). O estimador da matriz de
covariâncias é obtido por S = Wn /(n − 1).
Implementamos a função medcov() apresentada a seguir para obtermos o vetor de médias, a matriz de
somas de quadrados e produtos e a matriz de covariâncias. O argumento desta função deve ser uma matriz
de dados multivariados com n linhas (observações) e p colunas (variáveis). O programa resultante e um
exemplo são apresentados na sequência. Escolhemos um exemplo onde geramos uma matriz de dados de
uma normal multivariada.
# função para retornar o vetor de médias, a matriz de
# somas de # quadrados e produtos e a matriz de covariâncias
import numpy as np
def medcov(x):
n = [Link][0]
p = [Link][1]
if (n <= 1):
print('Matriz deve ter mais de 1 linha!')
return
xb = x[0,:]
w = [Link]((p, p))
for j in range(1, n):
82 CHAPTER 5. ALGORITMOS PARA MÉDIAS, VARIÂNCIAS E COVARIÂNCIAS
d = x[j,:] - xb
i = j + 1.0
w = w + (i - 1) * [Link](d, d) / i
xb = xb + d / i
s = w / (n - 1)
res = {'media': xb, 'covariância': s, 'SQP': w}
return res
# Exemplo
n = 1000
p = 5
x = [Link].multivariate_normal([Link](p),[Link](p),n)
medcov(x)
# comparar com resultado da nossa função
[Link](x, rowvar = False)
[Link](x, axis = 0)
5.4 Exercícios
1. Mostrar a equivalência existente entre as expressões (Equation 5.5) e (Equation 5.6).
2. Implementar em Python uma função para obtermos a soma de produtos, equação (Equation 5.7), e
covariância entre as n observações de duas variáveis e cujos valores estão dispostos em dois vetores X
5.4. EXERCÍCIOS 83
e Y. Criar uma matriz com n linhas e p = 2 colunas de algum exemplo e utilizar a função medcov()
para comparar os resultados obtidos.
3. Os coeficientes de assimetria e curtose amostrais são funções
√ das somas de potências dos desvios em
relação a média. O coeficiente de assimetria é dado por b1 = (W 3n /n)/(W 2n /n)3/2 e o coeficiente
de curtose é b2 = (W 4n /n)/(W 2n /n)2 . Implementar uma função Python, utilizando a função
medsq() para estimar o coeficiente de assimetria e curtose univariados.
4. Utilizar uma amostra em uma área de sua especialidade e determinar a média, variância, soma de
quadrados, soma de desvios ao cubo e na quarta potência. Determinar também os coeficientes de
assimetria e curtose.
84 CHAPTER 5. ALGORITMOS PARA MÉDIAS, VARIÂNCIAS E COVARIÂNCIAS
Chapter 6
Aproximação de Distribuições
Algoritmos para a obtenção de probabilidades foram alvos de pesquisas de muitas décadas e ainda
continuam sendo. A importância de se ter algoritmos que sejam rápidos e precisos é indescritível. A maior
dificuldade é que a maioria dos modelos probabilísticos não possui função de distribuição explicitamente
conhecida. Assim, métodos numéricos sofisticados são exigidos para calcularmos probabilidades. Um outro
aspecto é a necessidade de invertermos as funções de distribuição para obter quantis da variável aleatória
para uma probabilidade acumulada conhecida. Nos testes de hipóteses e nos processos de estimação por
intervalo e por região quase sempre utilizamos estes algoritmos indiretamente sem nos darmos conta disso.
Neste capítulo vamos introduzir estes conceitos e apresentar algumas ideias básicas de métodos gerais
para realizar as quadraturas necessárias. Métodos numéricos particulares serão abordados para alguns
poucos modelos. Também abordaremos separadamente os casos discreto e contínuo. Para finalizarmos
apresentaremos as principais funções pré-existentes do Python para uma série de modelos probabilísticos.
6.1 Introdução
Não temos muitas preocupações com a obtenção de probabilidades e quantis, quando utilizamos o Python,
pois a maioria dos modelos probabilísticos e das funções especiais já está implementada. Nosso objetivo
está além deste fato, pois nossa intenção é buscar os conhecimentos necessários para ir adiante e para
entendermos como determinada função opera e quais são suas potencialidades e limitações.
Vamos iniciar nossa discussão com a distribuição exponencial para chamarmos a atenção para a principal
dificuldade existente neste processo. Assim, escolhemos este modelo justamente por ele não apresentar
tais dificuldades. Considerando uma variável aleatória X com distribuição exponencial com parâmetro λ,
temos a função densidade
f (x) = λe−λx , x > 0 (6.1)
e a função de distribuição
F (x) = 1 − e−λx . (6.2)
Para este modelo probabilístico podemos calcular probabilidades utilizando a função de distribuição
(Equation 6.2) e obter quantis com a função de distribuição inversa que é dada por:
− ln(1 − p)
q = F −1 (p) = . (6.3)
λ
85
86 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
# Exemplo
import numpy as np
lamb = 0.1
x = [Link]([29.95732])
p = [Link]([0.95])
dexp(x, lamb)
pexp(x, lamb)
qexp(p, lamb)
array([0.005])
array([0.94999999])
array([-29.95732274])
Estas três funções foram facilmente implementadas, pois conseguimos explicitamente obter a função de
distribuição e sua inversa. Se por outro lado tivéssemos o modelo normal
1 (x − µ)2
f (x) = √ exp − (6.4)
2πσ 2 2σ 2
não poderíamos obter explicitamente a função de distribuição e muito menos a função inversa da função
de distribuição de probabilidade. Como para a grande maioria dos modelos probabilísticos encontramos os
6.2. MODELOS PROBABILÍSTICOS DISCRETOS 87
mesmos problemas, necessitamos de métodos numéricos para obter estas funções. Por essa razão iremos
apresentar alguns detalhes neste capítulo sobre alguns métodos gerais e específicos de alguns modelos. No
caso discreto podemos utilizar algoritmos sequenciais, porém como boa parte dos modelos probabilísticos
possuem relação exata com modelos contínuos, encontramos os mesmos problemas.
x
X n t
F (x) = p (1 − p)n−t . (6.6)
t=0
t
Vimos no capítulo 3, equação (Equation 3.14) que podemos obter as probabilidades acumuladas de forma
recursiva. Sendo P (X = 0) = (1 − p)n , obtemos as probabilidades para os demais valores de X, ou seja,
para x = 1, · · ·, n de forma recursiva utilizando a relação P (X = x) = P (X = x−1)[(n−x+1)/x][p/(1−p)].
Este mesmo algoritmo apropriadamente modificado é utilizado para obtermos as probabilidades acumuladas
e a inversa da função de distribuição. Se os valores de n e de p forem grandes, este algoritmo pode ser
ineficiente. Para o caso do parâmetro p ser grande (próximo de um) podemos utilizar a propriedade da
binomial dada por: se X ∼ Bin(n, p), então Y = n − X ∼ Bin(n, 1 − p). Assim, podemos, por exemplo,
obter P (X = x) de forma equivalente por P (Y = n − x) e P (X ≤ x) = P (Y ≥ n − x), que pode ser
reescrito por FX (x) = 1 − FY (n − x − 1), exceto para x = n, em que FX (x) = 1. Desta forma trocamos
de variável para realizar o cálculo e retornamos o valor correspondente a do evento original.
# função de probabilidade e distribuição
# da binomial(n, p)
def dpbinom(x, n, p = 0.5):
if p > 0.5:
pp = 1 - p
x = n - x
else:
pp = p
qq = 1 - pp
if (x < 0):
f = 0
else:
f = qq**n
r = pp / qq
g = r * (n + 1)
cdf = f
if x > 0:
u = 0
88 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
# Exemplo de uso
p = 0.5713131
x = 5 # retire a tolerância, qbinom falha
n = 20
res = dpbinom(x, n, p)
print(res)
prob = res['cdf']
qbinom(prob, n, p)
r = pp / qq
g = r * (n + 1)
x = 0
cdf = f
while ((probq - cdf) >= 1.e-11):
x += 1
f *= (g / x - r)
cdf += f
if p > 0.5:
if (probq - cdf) <= 1.e-11:
x = n - x - 1
else:
x = n - x
return x
# Exemplo de uso
p = 0.5713131
x = 5 # retire a tolerância, qbinom falha
n = 20
res = dpbinom(x, n, p)
print(res)
prob = res['cdf']
qqbinom(prob, n, p)
Desta forma podemos antever a importância de conhecermos algoritmos de integração numérica das
distribuições contínuas. Isso fica mais evidente quando percebemos que para grandes valores de n os
algoritmos recursivos são ineficientes e podem se tornar imprecisos. Na seção Section 6.3 apresentaremos
alguns métodos gerais de integração para funções contínuas. No script seguinte utilizamos a relação
(Equation 6.7}) para a obtenção da função de distribuição da binomial.
# função de probabilidade e distribuição da
# binomial(n, p) a partir da relação com a
# função beta incompleta - cdf da beta
from [Link] import beta
def pbinombeta(x, n, p = 0.5):
if p <= 0 or p >= 1:
print('p deve estar no intervalo (0, 1)')
return
if (x < n):
a = x + 1
b = n - x
90 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
cdf = 1 - [Link](p, a, b)
else:
cdf = 1
return cdf
# Exemplo
p = 0.5713131
n = 20
x = 3
pbinombeta(x, n, p)
import [Link] as sps
[Link](x, n, p)
np.float64(0.00013459360572398715)
np.float64(0.00013459360572400553)
A distribuição Poisson é a segunda que consideraremos. Existe relação da função de distribuição da
Poisson com a gama incompleta. Se uma variável aleatória discreta X com valores x = 0, 1, 2, 3, · · ·
possui distribuição Poisson com parâmetro λ, então podemos definir a função de probabilidade por:
λx e−λ
P (X = x) = . (6.8)
x!
Podemos obter probabilidades desta distribuição utilizando a relação da gama incompleta com a função
de distribuição Poisson
F (x|λ) = 1 − Iλ (α = x + 1), x = 0, 1, 2, · · · . (6.9)
sendo
1 x
Z
Ix (α) = e−t tα−1 dt,
Γ(α) 0
# Exemplo
lamb = 0.85
x = 3.0
pgama(x, lamb)
6.3. MODELOS PROBABILÍSTICOS CONTÍNUOS 91
np.float64(0.9888689674521238)
np.float64(0.9888689674521238)
1 (x − µ)2
f (x) = √ exp − . (6.10)
2πσ 2 2σ 2
A função de distribuição de probabilidade não pode ser obtida explicitamente e é definida por:
x
1 (t − µ)2
Z
F (x) = √ exp − dt. (6.11)
−∞ 2πσ 2 2σ 2
Em geral utilizamos a normal padrão para aplicarmos as quadraturas. Neste caso a média é µ = 0 e a
variância é σ 2 = 1. Assim, representamos frequentemente a densidade por ϕ(z), a função de distribuição
por Φ(z) e a sua inversa por Φ−1 (p), em que 0 < p < 1. A variável Z é obtida de uma transformação
linear de uma variável X normal por Z = (X − µ)/σ.
Vamos apresentar de forma bastante resumida a regra trapezoidal estendida para realizarmos quadraturas
de funções. Seja fi o valor da função alvo no ponto xi , ou seja, fi = f (xi ), então a regra trapezoidal é
dada por:
Z xn
h
f (x)dx = (f1 + fn ) + O(h3 f ′′ ), (6.12)
x1 2
em que h = xn − x1 e o termo de erro O(h3 f ′′ ) significa que a verdadeira resposta difere da estimada
por uma quantidade que é o produto de h3 pelo valor da segunda derivada da função avaliada em algum
ponto do intervalo de integração determinado por x1 e xn , sendo x1 < xn .
Esta equação retorna valores muito imprecisos para as quadraturas da maioria das funções de interesse
na estatística. Mas se utilizarmos esta função n − 1 vezes para fazer a integração nos intervalos (x1 , x2 ),
(x2 , x3 ), · · ·, (xn−1 , xn ) e somarmos os resultados, obteremos a fórmula composta ou a regra trapezoidal
estendida por:
xn
(xn − x1 )3 f ′′
Z
f1 fn
f (x)dx = h + f2 + f3 + · · · + fn−1 + +O . (6.13)
x1 2 2 n2
92 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
(b − a)
S= (f (a) + f (b)) .
2
return s
# função para executar quadraturas de
# funções definidas em func()
# até que uma determinada
# precisão tenha sido alcançada
def qtrap(func, a, b):
eps = 1.0e-11
nmax = 25
olds = -1.e30 # impossível valor para quadratura inicial
n = 1
fim = 0
s = -1.0e15 # valor arbitrário para iniciar o loop
while abs(s-olds) >= eps * abs(olds) and fim == 0:
olds = s
s = trapzd(func, a, b, s, n)
# evitar convergência prematura espúria
if n > 5:
if s < 0.0 and [Link](olds, 0.0):
fim = 1
n = n + 1
if n > nmax:
fim = 1
if n > nmax:
print('Limite de passos ultrapassado!')
return
return s
Devemos chamar a função qtrap() especificando a função de interesse func() e os limites de integração.
Assim, para a normal padrão devemos utilizar as seguintes funções:
# fdp da normal padrão
def dnorm(x):
fx = (1.0 / (2.0 * [Link])**0.5) * [Link](-x**2 / 2)
return fx
def pnorm(x):
p = qtrap(dnorm, 0.0, abs(x))
if x > 0:
p += 0.5
else:
p = 0.5 - p
return p
# exemplo
z = 1.96
print(pnorm(z))
print([Link](z))
0.9750021048512256
0.9750021048517795
É evidente que temos métodos numéricos gerais mais eficientes do que o apresentado. As quadraturas
gaussianas representam uma classe de métodos que normalmente são mais eficientes que este apresentado.
94 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
Eventualmente, existem métodos específicos para obtermos as integrais dos principais modelos probabilísti-
cos que são mais eficientes. A maior eficiência destes métodos específicos se dá por dois aspectos: velocidade
de processamento e precisão. Se exisitirem, no Python estas rotinas específicas já estão implementadas.
Como ilustração, podemos substituir a função que implementamos pnorm() pela pré-existente no Python
[Link]() do [Link]. Vamos ilustrar um destes algoritmos especializados para obtermos a
função de distribuição da normal padrão. O algoritmo de Hasting possui erro máximo de 1 × 10−6 e é
dado por:
se x ≤ 0
G
Φ(x) = (6.15)
1−G se x > 0
sendo G dado por
G = (a1 η + a2 η 2 + a3 η 3 + a4 η 4 + a5 η 5 )ϕ(x)
em que
1
η=
1 + 0,2316418|x|
e a1 = 0,319381530, a2 = −0,356563782, a3 = 1,781477937, a4 = −1,821255978 e a5 = 1,330274429.
O resultado da implementação deste procedimento é:
# CDF da normal padrão:
# aproximação de Hasting
def phnorm(x):
eta = 1 / (1 + abs(x) * 0.2316418)
a1 = 0.319381530; a2 = -0.356563782
a3 = 1.781477937; a4 = -1.821255978
a5 = 1.330274429
phi = 1 / (2 * [Link])**0.5 * [Link](-x * x / 2)
g = (a1*eta+a2*eta**2+a3*eta**3+a4*eta**4+a5*eta**5)*phi
if (x <= 0):
cdf = g
else:
cdf = 1 - g
return cdf
# exemplo
z = 1.96
p = phnorm(z)
p
0.9750021668514388
Com o uso desta função ganhamos em precisão, principalmente para grandes valores em módulo do
limite de integração superior e principalmente ganhamos em tempo de processamento. Podemos ainda
abordar os métodos de Monte Carlo, que são especialmente úteis para integrarmos funções complexas e
multidimensionais. Vamos apresentar apenas uma versão bastante rudimentar deste método. A ideia é
determinar um retângulo que engloba a função que desejamos integrar e bombardearmos a região com
pontos aleatórios (u1 , u2 ) provenientes da distribuição uniforme.
Contamos o número de pontos sob a função e determinamos a área correspondente, a partir da propor-
cionalidade entre este número de pontos e o total de pontos simulados em relação a área sob a função
6.3. MODELOS PROBABILÍSTICOS CONTÍNUOS 95
na região de interesse em relação à área total do retângulo. Se conhecemos o máximo da função fmax ,
podemos determinar este retângulo completamente. Assim, o retângulo de interesse fica definido pela base
(valor entre 0 e z1 em módulo, sendo z1 fornecido pelo usuário) e pela altura (valor da densidade no ponto
de máximo). Assim, a área
√ deste retângulo é A = |z1 |fmax . No√ caso da normal padrão, o máximo obtido
para z = 0 é fmax = 1/ 2π e a área do retângulo A = |z1 |/ 2π. Se a área sob a curva, que desejamos
estimar, for definida por A1 , podemos gerar números uniformes u1 entre 0 e |z1 | e números uniformes u2
entre 0 e fmax . Para cada valor u1 gerado calculamos a densidade f1 = f (u1 ). Assim, a razão entre as
áreas A1 /A é proporcional a razão n/N , em que n representa o número de pontos (u1 , u2 ) para os quais
u2 ≤ f1 e N o número total de pontos gerados. Logo, a integral é obtida por A1 = |z1 | × fmax × n/N , em
que z1 é o valor da normal padrão para o qual desejamos calcular a área que está compreendida entre 0 e
|z1 |, para assim obtermos a função de distribuição no ponto z1 , ou seja, para obtermos Φ(z1 ). Assim, se
z1 ≤ 0, então Φ(z1 ) = 0,5 − A1 e se z1 > 0, então Φ(z1 ) = 0,5 + A1 . Veja a seguinte figura ilustrativa:
# exemplo
96 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
z = 1.96
N = 1500000
mcpnorm(z, N)
print('Erro de MC: ',1 / N**0.5)
0.9756810408364969
Erro de MC: 0.0008164965809277261
Outra forma de obtermos uma aproximação da integral
abs(z)
1 abs(z)
Z Z
2
√ e−t /2 dt = ϕ(t)dt
0 2π 0
por Monte Carlo é gerarmos m números uniformes entre 0 e abs(z), digamos z1 , z2 , . . ., zm e obter
" m
#
abs(z)
1 X
Z
ϕ(t)dt ≈ ∆z ϕ(zi ) ,
0 m i=1
√
em que ϕ(t) = 1/ 2π × exp{−t2 /2} é a função densidade normal padrão avaliada no ponto t e ∆z =
abs(z) − 0. A ordem de erro desse processo é dada por O(m−1/2 ). O programa Python para obter o valor
da função de distribuição normal padrão Φ(z) utilizando essas ideias é apresentado a seguir. Por meio de
uma comparação dessa alternativa Monte Carlo com a primeira podemos verificar que houve uma grande
diminuição do erro de Monte Carlo no cálculo dessa integral, nessa nova abordagem. Muitas variantes e
melhorias nesse processo podem ser implementadas, mas nós não iremos discuti-las aqui.
# Quadratura da normal padrão via simulação
# Monte Carlo. Segunda forma de obter a
# integral: forma clássica
import numpy as np
import math
import [Link] as sps
import matplotlib
import [Link] as plt
def pmcnorm(z, m = 2000):
x = [Link](0, abs(z), m)
p = (1/(2 * [Link])**0.5)*[Link](-(x**2)/2)
p = abs(z) * [Link](p)
if z < 0:
p = 0.5 - p
else:
p += 0.5
return p
# Exemplo
m = 15000 # número de pontos muito inferior ao caso anterior
z = 1.96
pmcnorm(z, m) # Estimativa de Monte Carlo
p = [Link](z) # valor real
p
6.4. QUADRATURAS GAUSSIANAS 97
np.float64(0.9772224159439482)
np.float64(0.9750021048517795)
0.00816496580927726
Text(0.5, 0, 'Iterações')
1.05
1.00
P(Z<=z)
0.95
0.90
0.85
0 2000 4000 6000 8000 10000 12000 14000
Iterações
6.4 Quadraturas Gaussianas
As quadraturas gaussianas desempenham um papel preponderante na estatística, pois os principais
algoritmos de obtenção de probabilidades e quantis utilizam tais métodos implicitamente. Nesta seção,
apresentaremos, de forma bastante simplificada e sem aprofundar nos aspectos matemáticos mais técnicos,
o principal método de quadratura gaussiana. Os interessados em uma maior pormenorização desses
98 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
aspectos podem consultar inúmeros livros específicos. Recomendamos por exemplo a leitura de Quarteroni
et al. [2000].
A base dessas quadraturas são as interpolações baseadas em polinômios ortogonais. Seu uso extrapola o
tema das quadraturas numéricas, podendo ser usados para aproximar soluções de quadrados mínimos e
diferenciação numérica. Uma quadratura gaussiana de n pontos é aquela que fornece resultados exatos
para a integração de um polinômio de grau igual ou inferior a 2n − 1. O domínio das quadraturas é, sem
perda de generalidade e por convenção, assumido como sendo [−1, 1]. A regra geral para a quadratura
gaussiana de n pontos para uma função f (x), no domínio [−1, 1], é dada por
Z 1 n
X
f (x)dx ≈ wi f (xi ), (6.16)
−1 i=1
em que g(x) é tal que f (x) = w(x)g(x) e wi nesse caso são pesos alternativos.
√
Por exemplo,
√ para o caso específico da função f (x) = e2x , usando pesos w1 = w2 = 1 e nós x1 = 3/3 e
x2 = − 3/3, no domínio de −1 a 1 resulta em
Z 1 √ √
e2x dx ≈ f ( 3/3) + f (− 3/3) = 3,4882.
−1
É surpreendente que a soma de f (x1 ) + f (x2 ) resulte em valores exatos para polinômios de grau até
3, 2n − 1, pois n = 2. O que temos que fazer é encontrar mecanismos para obter os pesos e os nós de
integração. Para um intervalo de integração diferente de [−1, 1], ou seja, [a, b], b > a, temos que a seguinte
transformação não altera a precisão da integração
a+b
Rb b−a R1 b−a
f (x)dx = f x + dx
a 2 −1 2 2
a+b
b − a Pn b−a
≈ w
i=1 i f x i + .
2 2 2
Vale a pena ressaltar que os nós utilizados nas quadraturas são as raízes dos polinômios ortogonais de
alguma das famílias gaussianas, normalmente empregadas. Entre elas, podemos citar os polinômios de
Legendre, Laguerre, Hermite, Chebyshev e Jacobi. Descreveremos apenas a quadratura Gauss-Legendre
na sequência.
A quadratura denominada Gauss-Legendre é utilizada para intervalos de integração definidos por [−1, 1],
podendo ser extrapolado para intervalos mais gerais, se utilizarmos a expressão (Equation 6.18). Assim,
6.4. QUADRATURAS GAUSSIANAS 99
para o domínio [−1, 1], a quadratura Gauss-Legendre pode ser aplicada por
n
X
wi f (xi ) para a = −1 e b = 1
Z b
f (x)dx ≈ i=1
Xn (6.18)
a
wi g(yi ) para a e b reais quaisquer,
i=1
Z b
µ0 = w(x)dx
a
j
bj = p , j = 1, 2, . . . , n − 1.
4j 2 − 1
Definidas as quantidades necessárias a, b e µ0 , devemos aplicar a segunda parte do algoritmo para obter
os nós xi ’s e os pesos wi ’s, i = 1, 2, . . ., n. Utilizamos, para isso, o algoritmo de Golub–Welsch, pelo qual
determinamos a matriz J a partir dessas quantidades e determinamos seus autovalores e autovetores. A
matriz J é tridiagonal simétrica definida por
a1 b1 0
... ... ...
b1 a2 b2 0 ... ...
0 0
b2 a3 b3 ...
J = 0
.
... ... ... ... 0
... ... 0
bn−2 an−1 bn−1
... ... ... 0 bn−1 an
100 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
respectivamente.
Modificamos a notação usual para representar os elementos de uma matriz. Assim, verificamos que na
matriz P, o primeiro índice do elemento pij refere-se ao índice do i-ésimo autovetor e o segundo índice ao
j-ésimo elemento desse autovetor. Por isso, a notação não usual dos índices nessa matriz. Portanto, cada
coluna de P, corresponde a um dos n autovetores associados a cada um dos autovalores x∗i ordenados em
ordem decrescente.
Em seguida, podemos obter os nós (pontos), tomando simplesmente os autovalores em ordem crescente,
ou seja, o vetor dos nós x, com componentes xi , é definido por
x∗n
x1
x2 x∗n−1
x= .. = .. , (6.20)
. .
xn x∗1
em que x∗i é o i-ésimo autovalor da matriz tridiagonal J.
Podemos constatar que os nós são os autovalores dessa matriz apresentados em ordem reversa, ou seja,
em ordem crescente ao invés da habitual forma de apresentá-los, que é a ordem decrescente. Finalmente,
podemos obter os pesos como sendo uma função de µ0 e dos primeiros elementos de cada autovetor de
J, também considerado em ordem reversa, considerando as ordenações clássicas do maior para o menor
autovalor. Essa função busca realizar uma normalização adequada para garantir a validade da quadratura.
Assim, temos que os pesos das quadraturas são obtidos por
2
w1 pn,1
w2 p2n−1,1
2
w3 pn−2,1
w= . = µ0 . , (6.21)
.. ..
wn−1 p22,1
wn p21,1
if n < 0:
print('É necessário um número não negativo de nós!')
return
if n == 0:
return {'xi':[], 'wi': []}
i = [Link](1, n + 1)
i1 = [Link](i,n - 1)
mu0 = 2
b = i1 / [Link](4 * i1**2 - 1)
J = [Link]((n*n))
J[(n + 1) * (i1-1) + 1] = b
J[(n + 1) * i1 - 1] = b
J = [Link](n, n)
va, ve = [Link](J)
w = ve[0,:]
w = mu0 * w**2
return {'xi': va, 'wi': w}
# função polinomial
def p3(x):
return x**3-0.4*x**2+2.3*x+5
# Exemplos
n = 4
gausslegendquad(n)
n = 8
quadgl(e2x, n)
print('Valor exato: ',(e2x(1)-e2x(-1))/2)
quadgl(p3, 2)
print('Valor exato wolfram: ',9.73333)
# X ~ N(mu, sig)
import math
import scipy as sp
def pnormal(x, mu = 0, sigma = 1, n = 16):
if [Link](x, mu):
return 0.5
else:
if [Link](mu, 0) and [Link](sigma, 1):
z = x
else:
z = (x - mu) / sigma
res = quadglab([Link], n, 0, abs(z))
if (x < mu):
res = 0.5 - res
else:
res += 0.5
return res
np.float64(198.01986869689384)
Valor exato: 198.01986869690222
np.float64(35.73333333333333)
6.5. NEWTON-RAPHSON 103
np.float64(0.9750021048517794)
Nas quadraturas anteriores, podemos substituir nossa função que calcula os nós e os pesos, por outra, que
utiliza um método de obtenção de autovalores e autovetores mais eficiente. Esse método é aplicável a
matrizes simétricas tridiagonais. O script a seguir apresenta essa implementação. Precisamos apenas
fornecer um vetor com a diagonal da matriz e outro com os elementos das duas diagonais secundárias, que
são idênticos. A função utilizada foi a [Link].eigh_tridiagonal() da biblioteca scipy.
# Função para obter nós e pesos
# da quadratura Gauss-Legendre
# Essa versão dispensa a criação
# da matriz J, por aplicar um
# processo de obtenção de
# autovalores e autovetores em
# matrizes tridiagonais - + rápido
import numpy as np
import scipy as sp
def gausslegendquad2(n):
n = int(n)
if n < 0:
print('É necessário um número não negativo de nós!')
return
if n == 0:
return {'xi':[], 'wi': []}
i = [Link](1, n + 1)
i1 = [Link](i,n - 1)
mu0 = 2
b = i1 / [Link](4 * i1**2 - 1)
a = [Link](n)
va, ve = [Link].eigh_tridiagonal(a, b)
w = ve[0,:]
w = mu0 * w**2
return {'xi': va, 'wi': w}
n = 4
gausslegendquad2(n)
6.5 Newton-Raphson
Vamos apresentar o método numérico de Newton-Raphson para obtermos a solução da equação
z =Φ−1 (p),
em que 0 < p < 1 é o valor da função de distribuição da normal padrão, Φ−1 (p) é a função inversa da
função de distribuição normal padrão no ponto p e z o quantil correspondente, que queremos encontrar
dado um valor de p. Podemos apresentar esse problema por meio da seguinte equação
104 CHAPTER 6. APROXIMAÇÃO DE DISTRIBUIÇÕES
Φ(z) − p =0,
em que Φ(z) é a função de distribuição normal padrão avaliada em z. Assim, nosso objetivo é encontrar
os zeros da função f (z) = Φ(z) − p. Em geral, podemos resolver essa equação numericamente utilizando o
método de Newton-Raphson, que é um processo iterativo. Assim, devemos ter um valor inicial para o
quantil para iniciarmos o processo e no (n + 1)-ésimo passo do processo iterativo podemos atualizar o
valor do quantil por
f (zn )
zn+1 = zn − , (6.22)
f ′ (zn )
em que f ′ (zn ) é derivada de primeira ordem da função para a qual queremos obter as raízes avaliada no
ponto zn . Para esse caso particular temos que
[Φ(zn ) − p]
zn+1 = zn − , (6.23)
ϕ(zn )
sendo ϕ(zn ) a função densidade normal padrão. Como valores iniciais usaremos uma aproximação grosseira,
pois nosso objetivo é somente demonstrar o método. Assim, se p for inferior a 0,5 utilizaremos z0 = −0,1,
se p > 0,5, utilizaremos z0 = 0,1. Obviamente se p = 0, a função deve retornar z = 0. A função qPhi a
seguir retorna os quantis da normal, dados os valores de p entre 0 e 1, da média real µ e da variância
real positiva σ 2 , utilizando para isso o método de Newton-Raphson e a função phnorm de Hasting para
obtermos o valor da função de distribuição normal padrão.
# função auxiliar para retornar o valor da
# função densidade normal padrão
def phi(z):
return (1 / (2 * [Link])**0.5 * [Link](-z * z / 2))
# Exemplo
p = 0.975
mu = 0
sig = 1
qphi(p, mu, sig, phnorm)
[Link](p, mu, sig) #para fins de comparação
np.float64(1.959963984540054)
Poderíamos, ainda, ter utilizado o método da secante, uma vez que ele não necessita da derivada de
primeira ordem, mas precisa de dois valores iniciais para iniciar o processo e tem convergência mais lenta.
O leitor é incentivado a consultar Press et al. [1992] para obter mais detalhes. Também poderia ter sido
usada a função pnorm, que utiliza o método trapezoidal. Nesse caso a precisão seria maior, mas o tempo
de processamento também é maior. Para isso bastaria usar na chamada como valor do argumento func
a função pnorm. Também seria possível chamar mcpnorm ou pmcnorm ou qualquer outra implementação
eficiente que tivermos para a quadratura.
Felizmente o Python também possui rotinas pré-programadas para este e para muitos outros modelos prob-
abilísticos, que nos alivia da necessidade de programar rotinas para obtenção das funções de distribuições
e inversas das funções de distribuições dos mais variados modelos probabilísticos existentes.
6.7 Exercícios
1. Comparar a precisão dos três algoritmos de obtenção da função de distribuição normal apresentados
neste capítulo. Utilizar a função [Link]() como referência.
2. Utilizar diferentes números de simulações Monte Carlo para integrar a função de distribuição normal
e estimar o erro de Monte Carlo relativo e absoluto máximo cometidos em 30 repetições para cada
tamanho. Utilize os quantis 1,00, 1,645 e 1,96 e a função [Link]() como referência.
1
f (x) =
π(1 + x2 )
Utilizar o método trapezoidal estendido para implementar a obtenção dos valores da distribuição
Cauchy. Podemos obter analiticamente também a função de distribuição e sua inversa. Obter tais
funções e implementá-las no Python. Utilize as funções [Link] pré-existentes para checar os
resultados obtidos.
4. Utilizar o método Monte Carlo descrito nesse capítulo, para obter valores da função de distribuição
Cauchy, apresentada no exercício proposto 3. Utilizar alguns valores numéricos para ilustrar e
comparar com funções implementadas em Python para esse caso. Qual deve ser o número mínimo
de simulações Monte Carlo requeridas para se ter uma precisão razoável.
Chapter 7
por [Link](n, x). Já á listagem das combinações de n tomados x a x é obtida pelo comando
combinations(range(1,n+1), x) do pacote itertools que deve ser importado. O primeiro argumento
da função combinations é uma lista, que no caso foi de uma lista indo de 1 a n, em que usamos a função
107
108 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
range para obter essa lista. Podemos, em vez disso, usar qualquer outra lista. Podemos usar ainda a
função comb do [Link] para obtermos o número de combinações de forma alternativa.
O comando combinations(n, x) do pacote itertools retorna um objeto que deve ser convertido para
uma lista com nx elementos, em que cada elemento da lista corresponde a uma tupla de x elementos (a
combinação).
O script a seguir ilustra um caso particular destes comandos:
from itertools import combinations
from scipy import special as sps
import math
n = 4
x = 2
[Link](n, x)
[Link](n, x) # resultado é um float
comb = list(combinations(range(1,n+1), x))
print(comb)
# listando a combinação 0 da lista
print(comb[0])
len(comb) # número de combinações alternativo
# retorna uma combinação de letras
list(combinations(['a','b','c'], 2))
6
np.float64(6.0)
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]
(1, 2)
6
[('a', 'b'), ('a', 'c'), ('b', 'c')]
perm = list(permutations(range(1,n+1)))
print(perm)
# arranjos
x = 2 # x < n
[Link](n, x)
arran = list(permutations(range(1,n+1), x))
print(arran)
# Imprime as permutações
for i in list(perm):
print(i)
type(perm[0])
[(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)]
6
[(1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2)]
(1, 2, 3)
(1, 3, 2)
(2, 1, 3)
(2, 3, 1)
(3, 1, 2)
(3, 2, 1)
tuple
Para sortearmos permutações ou arranjos de forma aleatória, o comando [Link](range(1,n+1),n)
retorna uma permutação de n tomados n a n. Se quisermos um sorteio de n tomados x a x
[Link](range(1,n+1),x), ou seja, um arranjo de um resultado com x elementos sem repetição
(amostra sem reposição). Não é referente ao caso em questão (permutações e arranjos), mas se
desejarmos amostragens com reposição, como requerido nos métodos bootstrap, usamos a função
[Link](range(1,n+1),x), em que x é um valor entre 1 e n. Ambos os comando requerem que
importemos a biblioteca random do Python. O script a seguir ilustra estes comandos:
# Obter amostras aleatórias de
# permutações ou arranjos
# (sem reposição), ou amostras com
# reposição
import random
n = 5
# amostras de permutações
amostra = [Link](range(1,n+1), n)
print(amostra)
x = 3 # x < n
# amostras de arranjos
amostra = [Link](range(1,n+1), x)
print(amostra)
# amostras com reposição
x = 3
amos = [Link](range(1,n+1), k=x)
print(amos)
[2, 4, 3, 5, 1]
[5, 3, 2]
[4, 3, 1]
110 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
7.3 Contagem
O princípio fundamental da contagem, conhecido por princípio multiplicativo, é utilizado para encontrar
o número de possibilidades para um evento constituído de k etapas. As etapas devem ser sucessivas e
independentes. Se um evento tem duas (k = 2) etapas, sendo que a primeira possui n1 possibilidades e a
segunda, n2 possibilidades, então existem n1 × n2 possibilidades.
Assim, o princípio fundamental da contagem é a multiplicação das opções de cada etapa para determinar
o total de possibilidades. Esse conceito é importante para a análise combinatória, área da matemática que
reúne os métodos para resolução de problemas incluindo a contagem entre eles e, por isso, é muito útil na
investigação para determinar a probabilidade de fenômenos naturais.
Para obtermos todas as possibilidades em k etapas, cada uma com a mesma quantidade n de possibilidades
criamos uma classe Python com algumas funções. A função contag(n, k) retorna todas as contagens
a partir de n possibilidades em k etapas, como, por exemplo, com n = 2 sexos de animais em k = 3
nascimentos, ou n = 2 portas de um prédio em k = 2 possibilidades, correspondentes às entradas e às
saídas do prédio. Também fizemos o mesmo quando temos uma lista (iterador qualquer) de tamanho k
correspondendo às k etapas em que cada elemento refere-se ao número de possibilidades daquela etapa.
A pergunta, do segundo exemplo, de quais maneiras diferentes uma pessoa pode entrar e sair do edifício
com 2 entradas e duas saídas é o que pretendemos listar (enumerar). A função proposta usa um processo
recursivo, em que em cada chamada é atualizada a contagem em particular. A recursividade ocorre na
função permuta ou permutav, para o caso de uma lista de etapas.
# Classe para obter as contagens de x
# possibilidades em cada k etapas: x, k vezes
# ou x = [n1,n2,...,nk](lista com k nós).
# Usa overload e ilustra a criação
# de uma classe Python
from multipledispatch import dispatch
import numpy as np
class Cont:
def permuta(self,z, pos, permn, n, k, prim):
if prim:
z = [Link](k + 1, 1)
prim = False
permun = [Link](permn,[z[1:(k+1)]],axis=0)
else:
permun = permn
if z[pos] < n:
z[pos] = z[pos] + 1
permun = [Link](permun,[z[1:(k+1)]],axis=0)
else:
nachei = True
if pos == k:
i = pos - 1
paratras = True
else:
i = pos + 1
paratras = False
while nachei and i >= 1 and i <= k:
if z[i] < n:
z[i] = z[i] + 1
if paratras:
7.3. CONTAGEM 111
z[(i+1):(k+1)] = 1
permun = [Link](permun,[z[1:(k+1)]],axis=0)
pos = k
nachei = False
else:
if (pos == k):
i = i - 1
else:
i = i + 1
if [Link](z[1:(k+1)]) == n * k:
return permun
else:
return [Link](z, pos, permun, n, k, prim)
def permutav(self,z, pos, permn, x, k, prim):
if prim:
k = len(x)
pos = k
z = [Link](k + 1, 1)
prim = False
permun = [Link](permn,[z[1:(k+1)]],axis=0)
else:
permun = permn
if z[pos] < x[pos-1]:
z[pos] = z[pos] + 1
permun = [Link](permun,[z[1:(k+1)]],axis=0)
else:
nachei = True
if pos == k:
i = pos - 1
paratras = True
else:
i = pos + 1
paratras = False
while nachei and i >= 1 and i <= k:
if z[i] < x[i-1]:
z[i] = z[i] + 1
if paratras:
z[(i+1):(k+1)] = 1
permun = [Link](permun,[z[1:(k+1)]],axis=0)
pos = k
nachei = False
else:
if pos == k:
i = i - 1
else:
i = i + 1
if [Link](z[1:(k+1)]) == sum(x):
return permun
else:
return [Link](z, pos, permun, x, k, prim)
@dispatch(int, int)
112 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
[[1 1 1]
[1 1 2]
[1 1 3]
[1 2 1]
[1 2 2]
[1 2 3]
[1 3 1]
[1 3 2]
[1 3 3]
[2 1 1]
[2 1 2]
[2 1 3]
[2 2 1]
[2 2 2]
[2 2 3]
[2 3 1]
[2 3 2]
[2 3 3]
[3 1 1]
[3 1 2]
[3 1 3]
[3 2 1]
[3 2 2]
[3 2 3]
[3 3 1]
[3 3 2]
7.3. CONTAGEM 113
[3 3 3]]
array([[1, 1, 1],
[1, 1, 2],
[1, 2, 1],
[1, 2, 2],
[1, 3, 1],
[1, 3, 2],
[1, 4, 1],
[1, 4, 2],
[2, 1, 1],
[2, 1, 2],
[2, 2, 1],
[2, 2, 2],
[2, 3, 1],
[2, 3, 2],
[2, 4, 1],
[2, 4, 2],
[3, 1, 1],
[3, 1, 2],
[3, 2, 1],
[3, 2, 2],
[3, 3, 1],
[3, 3, 2],
[3, 4, 1],
[3, 4, 2]])
No exemplo anterior tivemos várias novidades. Criamos uma classe pela primeira vez. Nesta classe
implementamos quatro funções. Entre elas, as duas com a função de realizar as listagens de todas
as contagens foram implementadas com chamadas recursivas. As outras duas, foram duas versões
de uma mesma função com diferentes argumentos, a função contag. Para isso usamos a técnica de
overloading. Assim, usamos o pacote multipledispatch de onde importamos from multipledispatch
import dispatch. O comando @dispatch(int, int) ou @dispatch(list) antes da definição da função
contag indica que ela tem diferentes argumentos. A primeira definição possui dois argumentos inteiros e
a segunda, um argumento, correspondente a uma lista. Para exemplificar, criamos um objeto da classe
Cont, objeto res, por meio do qual chamamos o método contag duas vezes com argumentos diferentes, o
que faz com que ao ser executado, o interpretador escolha a opção apropriada.
Se, por exemplo, quiséssemos obter todas as possibilidades de nascimento quanto ao sexo de famílias de
até 4 filhos, teríamos o seguinte espaço amostral de nosso experimento com 24 = 16 elementos:
# Espaço amostral dos filhos
# quanto ao sexo
# com n = 4 e 0 <= x <= n
filhos = Cont()
n = 4
res = [Link](2, n) - 1
sexo = ['F','M']
type(res[0])
print([Link](sexo)[res])
[Link]
[['F' 'F' 'F' 'F']
['F' 'F' 'F' 'M']
114 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
(1, 2, 3),
(1, 3, 1),
(1, 3, 2),
(1, 3, 3),
(2, 1, 1),
(2, 1, 2),
(2, 1, 3),
(2, 2, 1),
(2, 2, 2),
(2, 2, 3),
(2, 3, 1),
(2, 3, 2),
(2, 3, 3),
(3, 1, 1),
(3, 1, 2),
(3, 1, 3),
(3, 2, 1),
(3, 2, 2),
(3, 2, 3),
(3, 3, 1),
(3, 3, 2),
(3, 3, 3)]
[(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2), (1, 3, 1), (1, 3, 2), (1, 4, 1), (1, 4, 2), (2, 1, 1), (
{1, 2, 3, 4, 5, 6, 7}
{1, 2, 3, 4, 5, 6, 7, 8}
{1, 2, 3, 4, 5, 6, 7}
{1, 2, 3, 4, 5, 6, 7}
1
{2, 3, 4, 5, 6, 7}
O Python inclui algumas operações com conjuntos, sendo algumas delas:
• [Link](x, y): união dos conjuntos x e y, que corresponde ao conjunto de todos os elementos
presentes em x e em y, considerando apenas uma vez aqueles elementos com multiplicidade maior
que 1;
• [Link](x,y): interseção dos conjuntos x e y, que corresponde ao conjunto dos elementos
presentes simultaneamente em x e y;
• [Link](x,y): conjunto da diferença entre os conjuntos x e y, consistindo no conjunto de
todos os elementos de x que não estão em y;
• x == y ou x != y: testa se dois conjuntos x e y são iguais ou se são diferentes, respectivamente;
• c in y: pertencimento, ou seja, testa se c é um elemento do conjunto y;
7.4. CONJUNTOS EM PYTHON 117
{1, 2, 3, 4, 5}
{1, 2, 3, 4, 5}
{3, 4, 5}
{3}
{4, 5}
False
True
(1, 2, 3)
tuple
{1, 2, 4, 5}
Outras opções de operações básicas com conjuntos são as seguintes:
• [Link](x): retorna True se y está contido ou é igual ao conjunto x;
• [Link](y): retorna True se x contém ou é igual ao conjunto y;
• [Link](y): retorna True se x e y não tiverem elementos comuns, ou seja, se a interseção for
um conjunto vazio.
# outras operações com conjuntos
x = {1,2,3,4,5}
y = {1, 3, 5}
# y está contido em x
[Link](x)
# x contém y
[Link](y)
# se x e y são disjuntos
[Link](y)
True
118 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
True
False
A365,n 365!
P (A) = = .
365n (365 − n)!365n
Desejamos a probabilidade do evento complementar, que é dada por P (Ac ) = 1 − P (A). Logo,
A365,n 365!
P (Ac ) =1 − =1− ,
365n (365 − n)!365n
que é a probabilidade de haver pelo menos uma coincidência de aniversários na mesma data.
Implementamos duas formas. Na primeira usamos a fórmula completa. Tomamos o logaritmo (usamos a
função gammaln da [Link]) e ao final recuperamos tomando o exp do resultado. Na segunda alter-
nativa usamos a biblioteca math, em que o total de arranjos A365,n é obtido pelo comando [Link](365,
n).
n = [Link](1, 100+1)
# Gráfico
[Link](n, ca(n))
[Link]()
7.5. ALGUNS PROBLEMAS DE PROBABILIDADE 119
1.0
0.8
0.6
0.4
0.2
0.0
0 20 40 60 80 100
1.0
0.8
0.6
0.4
0.2
0.0
0 20 40 60 80 100
No segundo exemplo, vamos considerar a probabilidade de que uma comissão de tamanho m ao ser retirada
de um grupo de tamanho n aleatoriamente não contenha representantes de um grupo existente específico
de tamanho k entre os n elementos do grupo todo. Essa probabilidade é:
n−k
P (Não ser representado) = n .m
O programa Python foi implementado e ilustramos com o exemplo dos senadores do Brasil, que possui
representantes das diferentes unidades federativas brasileiras. No Brasil temos n = 27 × 3 senadores, ou
seja, 81 senadores no total, sendo 3 de cada estado ou do Distrito Federal. Se vamos formar uma comissão
aleatória de m = 10 senadores, qual é a probabilidade de que uma unidade federativa qualquer contendo
k = 3 senadores não seja representada. Vimos por meio da análise da fórmula anterior que temos m n
[Link]()
# probabilidade de MG e DF
# não serem representados na comissão
k = 6
m = 10
pnr(n,m,k)
np.float64(0.6698898265353962)
1.0
0.8
0.6
0.4
0.2
0.0
0 10 20 30 40 50 60 70 80
np.float64(0.44129815640475195)
Nosso próximo exemplo refere-se às probabilidades em uma mão de pôquer, construindo as possibilidades
e probabilidades. Neste jogo podemos ter uma mão sem nada (sem pares, sem dois pares, sem trincas,
etc.), uma com um par, uma com dois pares, uma com uma trinca, uma com uma sequência (qualquer de
naipe), uma com flush (diferentes valores que não estão em sequência do mesmo naipe), uma com um full
house (trinca e par), uma com o four (quadra), uma com uma sequência do mesmo naipe e uma sequência
real do mesmo naipe (royal flush - do 10 ao Ás).
O número de possibilidades totais de distribuir 52 em um sorteio de 5 cartas (uma mão) é:
52 52!
= = 2.598.960.
5 5!(52 − 5)!
• Para obtermos uma mão sem valor, temos que entender que o baralho é constituído de 52 cartas,
sendo 13 valores das cartas (Ás, 1, 2, · · ·, 10, Valete, Dama e Rei) e 4 naipes para cada valor
(espadas, paus, ouros e copas). O número de possibilidades de uma mão de cinco cartas (sorteio de
5 cartas sem reposição) conter uma mão sem valor, ouseja, sem pares, sem trincas, sem quadras,
sem sequências de naipes diferentes ou sem sequências do mesmo naipe é:
5 5
13 4 4 13
− 10 −4 + 4 × 10 = = 1.302.540,
5 1 1 5
ou seja, é o número de cartas de diferentes valores com qualquer um dos 4 naipes 13 5 4 , subtraído de
5
todas as possibilidades das sequências de naipes diferentes 10 × 45 e das possibilidades de ter cartas de
122 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
diferentes valores com o mesmo naipe 4 135 . Este valor deve ser adicionado das dez possíveis sequências
para cada um dos naipes 4 × 10. Esta adição se dá em razão de termos retirado as possibilidades de de ter
cartas de diferentes valores com o mesmo naipe, que incluí as 10 sequências do mesmo naipe possíveis e
novamente foi retirada quando eliminamos as possibilidades de termos 5 cartas de valores diferentes com
o mesmo naipe, que incluí as 10 sequências do mesmo naipe possíveis. Daí precisamos adiciona-las, uma
vez que elas foram retiradas duas vezes.
• Para a mão de um par simples temos:
13 4 12 3
4 ,
1 2 3
remanescentes (não escolhido para o par) e 43 , escolhe os naipes dos 3 valores diferentes do par.
• Para os dois pares temos:
13 4 12 4 11
4,
1 2 1 2 1
em que 13 1 escolhe o valor do primeiro par, 2 escolhe os naipes para o primeiro par, 1 escolhe o
4 12
segundo par entre os 12 valores remanescentes, 42 escolhe os naipes para o segundo par e 11 1 escolhe os
valores entre os 11 valores remanescentes (não escolhido para os 2 pares) e 4, escolhe os naipes do valor
diferente dos dois pares.
• Para uma tripla temos:
13 4 12 2
4 ,
1 3 2
10 × 45 − 10 × 4,
13
× 4 − 10 × 4,
5
em que 13 5 escolhe os 5 diferentes valores entre os 13 e 4 escolhe um dos naipes para os 5 valores. O
resultado incluí as sequências do mesmo naipe, que devem ser retiradas, que correspondem à 10 × 4
(sequências de mesmo naipe incluindo as reais).
7.5. ALGUNS PROBLEMAS DE PROBABILIDADE 123
13 4 12 4
1 3 1 2
13 12
4
1 1
em que 13 1 escolhe o valor para a quadra que necessariamente terá uma de cada naipe, 1 escolhe um
12
valor para a carta remanescente dos 12 valores restantes e 4 escolhe um dos naipes da carta remanescente,
que não é o da quadra.
• Para as sequências do mesmo naipe (straight flush):
10 × 4 − 4,
pois são 10 sequências para um dos 4 naipes, subtraída das 4 sequências reais entre elas, do 10 ao Ás de
cada um naipe entre os quatro naipes possíveis.
• A sequência real (royal straight flush): tem uma única sequência possível do 10 ao Ás para cada um
dos quatro naipes, totalizando 4 possíveis sequências reais.
O programa a seguir cria um dataframe com todas essas contagens e calcula as probabilidades dividindo-as
pelo número total de combinações:
#Mão de pôquer
# probabilidades
import pandas as pd
import numpy as np
from scipy import special as sps
poker = {
'none': [Link](13,5)*4**5 - 10*4**5 - 4*[Link](13,5) + 10*4,
'pair': 13*[Link](4,2)*[Link](12,3)*4**3,
'[Link]': [Link](13,2)*[Link](4,2)**2*11*4,
'triple': 13*[Link](4,3)*[Link](12,2)*4**2,
'straight': 10*4**5 - 10*4,
'flush': 4*[Link](13,5) - 10*4,
'[Link]': 13*[Link](4,3)*12*[Link](4,2),
'four': 13*[Link](4,4)*12*4,
'[Link]': 10*4 - 4,
'[Link]': 4 }
sum([Link]()) - [Link](52,5) # conferir
poker
mão = list([Link]())
possibilidades = list([Link]())
124 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
np.float64(0.0)
{'none': np.float64(1302540.0),
'pair': np.float64(1098240.0),
'[Link]': np.float64(123552.0),
'triple': np.float64(54912.0),
'straight': 10200,
'flush': np.float64(5108.0),
'[Link]': np.float64(3744.0),
'four': np.float64(624.0),
'[Link]': 36,
'[Link]': 4}
1.0
Nosso próximo exemplo refere-se à probabilidade da máxima diferença de uma lista enumerada ordenada.
Obter o maior gap (salto) (diferença máxima) entre dois números consecutivos no sorteio sem reposição
de m números entre os n primeiros inteiros de 1 a n, para m < n. No Python, o comando [Link](x),
retorna as diferenças entre números consecutivos da lista x.
O programa a seguir, usa a função checkgap para retornar a máxima diferença de um interador qualquer.
A função maxgap, recebe n, m e k e calcula as probabilidades de cada caso. Para isso ela percorre todas as
combinações geradas por combinations retornando o máximo de cada combinação obtida em y e com o
comando mean(y == k) é obtida a probabilidade exata de P (X = k), sendo X a variável que representa a
maior diferença entre números inteiros consecutivos de uma amostra sem reposição de tamanho m obtida
entre n números inteiros, para m de 1 a n e k o valor da máxima diferença (valor da variável aleatória)
para o qual queremos calcular a probabilidade de ocorrência. Se m = 2, então a máxima diferença será
considerada a diferença consecutiva, pois a diferença em cada resultado do espaço amostral é única, em
razão de termos amostras de tamanho 2.
# Max gap probabilidades
from itertools import combinations
import pandas as pd
def checkgap(cmb):
x = [Link](cmb)
return max(x)
7.5. ALGUNS PROBLEMAS DE PROBABILIDADE 125
# exemplo de uso
n = 6
m = 3
k = 2
maxgap(n,m,k)
maxgapsx(n,m)
np.float64(0.4)
x P(X=x)
0 1 0.2
1 2 0.4
2 3 0.3
3 4 0.1
Nosso último exemplo é para a probabilidade das somas das faces no lançamento de n dados e da diferença
em valor absoluto das faces no lançamento de 2 dados. Para este caso, podemos usar uma versão similar,
ou aplicando a função contag anteriormente apresentada ou aplicando a função product da biblioteca
intertools. Optamos por usar a função contag da classe cont.
Com uso de alguns pequenos detalhes adicionais, obtivemos os resultados para os dois casos. O segundo
caso, da diferença, por razões óbvias, são resultantes do lançamento de apenas 2 dados. Nos dois casos
apresentamos também os resultados das probabilidades obtendo todas as possibilidades que corresponde
ao espaço amostral do experimento aleatório. Temos uma função para obter uma probabilidade para um
valor específico da variável aleatória e outra para todos os valores de probabilidade relativos ao suporte da
126 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
variável aleatória, nos dois casos, o da soma de n dados e o da diferença de dois dados. Os resultados
estão apresentados no script a seguir:
# Probabilidades para o lançamento de
# n dados e obtenção da soma e também
# para a diferença absoluta de 2 dados
# somafaces: retorna P(X=x)
# depende da class cont
def somafaces(n, x):
if x<n or x>6*n:
print('x deve estar entre n e 6n!')
return
res = Cont()
perm = [Link](6, n)
y = [Link](0)
for pmb in perm:
y = [Link](y,[Link](pmb))
return [Link](y == x)
def somafacessx(n):
if x<n or x>6*n:
print('x deve estar entre n e 6n!')
return
res = Cont()
perm = [Link](6, n)
y = [Link](0)
for pmb in perm:
y = [Link](y,[Link](pmb))
res = {'x':[],'P(X=x)':[]}
for k in range(n,6*n+1):
res['P(X=x)'].append([Link](y == k))
res['x'].append(k)
return [Link](res)
# Probabilidades da diferença entre
# duas faces de dois dados equilibrados
# em módulo - retorna P(X=x), X=|i-j|
def difabs2faces(x):
if x<0 or x>5:
print('x deve estar entre 0 e 5!')
return
res = Cont()
perm = [Link](6, 2)
y = [Link](0)
for pmb in perm:
y = [Link](y,abs([Link](pmb)))
return [Link](y == x)
# Probabilidades da diferença entre
# duas faces de dois dados equilibrados
# em módulo - retorna P(X=x), para todo
# o suporte de x
def difabs2facessx():
if x<0 or x>5:
print('x deve estar entre 0 e 5!')
return
7.5. ALGUNS PROBLEMAS DE PROBABILIDADE 127
res = Cont()
perm = [Link](6, 2)
y = [Link](0)
for pmb in perm:
y = [Link](y,abs([Link](pmb)))
res = {'x':[],'P(X=x)':[]}
for k in range(0,5+1):
res['P(X=x)'].append([Link](y == k))
res['x'].append(k)
return [Link](res)
# Exemplo
n = 4
x = 4
print('P(X =',x,') = ',somafaces(n,x))
dist = somafacessx(n)
dist
sum(dist['P(X=x)']) # checando
print('P(X =',x,') = ',difabs2faces(x))
difabs2facessx()
P(X = 4 ) = 0.0007716049382716049
x P(X=x)
0 4 0.000772
1 5 0.003086
2 6 0.007716
3 7 0.015432
4 8 0.027006
5 9 0.043210
6 10 0.061728
7 11 0.080247
8 12 0.096451
9 13 0.108025
10 14 0.112654
11 15 0.108025
12 16 0.096451
13 17 0.080247
14 18 0.061728
15 19 0.043210
16 20 0.027006
17 21 0.015432
18 22 0.007716
19 23 0.003086
20 24 0.000772
1.0
P(X = 4 ) = 0.1111111111111111
x P(X=x)
0 0 0.166667
128 CHAPTER 7. CONJUNTOS E ELEMENTOS DE ANÁLISE COMBINATÓRIA EM PYTHON
x P(X=x)
1 1 0.277778
2 2 0.222222
3 3 0.166667
4 4 0.111111
5 5 0.055556
7.6 Exercícios
1. Obter uma função para obter por meio de simulação Monte Carlo as probabilidades dos lançamentos
de n dados e obtenção da soma e também para as diferenças em módulo entre as faces nos lançamentos
de dois dados. Para fazer isso escolha um grande número de repetições e em cada uma delas simule
os resultados dos lançamentos dos dados. Obtenha os resultados das variáveis aleatórias e ao final
calcule as proporções de cada ocorrência, conforme fizemos de forma exata anteriormente. Estas
proporções são estimativas empíricas das probabilidades almejadas.
2. Obtenha também por Monte Carlo as probabilidades empíricas de uma mão de pôquer, simulando
um grande número de mãos. Para isso é preciso criar um baralho, que pode ser feito com um
dicionário e depois amostrar as mãos de cinco cartas e calcular as probabilidades empíricas de cada
possível resultado. Confrontar com as probabilidades exatas obtidas neste capítulo. O que você
espera que ocorra com essa comparação quando aumenta-se o número de simulações?
Chapter 8
Neste capítulo, pretendemos apresentar os conceitos básicos de bootstrap de uma forma bastante simples.
Vamos fazer algumas funções simples em alguns casos particulares.
8.1 Introdução
Um dos principais temas envolvendo os métodos computacionalmente intensivos e, talvez, o maior
responsável pela popularidade desses métodos é a técnica bootstrap. O método bootstrap envolve
reamostragens com reposição da amostra original. Inspirado na teoria da amostragem, o bootstrap é
utilizado para a realização de testes de hipótese, estimação de parâmetros por intervalos e estimação de
erros padrões.
A grande vantagem de se utilizar os métodos de reamostragens, como o bootstrap, para realizar inferência
é quando não conhecemos a distribuição de probabilidade da população e o modelo normal não é adequado
para os dados ou resíduos. O problema surge, quando violamos as condições assumidas para a aplicação
de um teste ou obtenção de um intervalo ou região de confiança.
Neste capítulo, estudaremos os métodos bootstrap para realizar inferência sobre parâmetros de interesse.
129
130 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
A ideia do bootstrap é substituir a distribuição desconhecida da população pela distribuição empírica. Por
essa razão, que denominamos esse método de bootstrap não-paramétrico. Ao obtermos um série finita de
reamostragens de mesmo tamanho da amostra original e construirmos uma amostragem da distribuição
do estimador, estamos realizando um processo de mímica em relação à teoria de amostragem.
Esse processo se baseia no fato de que a amostra obtida da população contém toda a informação disponível
dessa população subjacente. Então, ela passa a representar a “população” para o processo de reamostragem.
Efron (1979) foi o pesquisador que organizou as teorias a respeito desse método e conectou o bootstrap
não-paramétrico com as técnicas estatísticas para estimar erros, aceitas desde os anos de 1930, tais como
jackknife e o método delta. Tanto a técnica bootstrap quanto os métodos de permutação assumem apenas
que as variáveis aleatórias atendam ao princípio de serem permutáveis (exchangeable), que é uma suposição
mais fraca que sejam independentes e identicamente distribuídas. Por permutáveis (intercambiáveis)
devemos entender que a distribuição de probabilidade de quaisquer k observações consecutivas (k = 1, 2,
. . ., n) não se altera quando a ordem das observações é trocada por meio de alguma permutação.
O algoritmo geral das reamostragens bootstrap, para estimar um pa-râ-me-tro θ desconhecido, cuja
distribuição de probabilidade é f (de Y ), considerando, ainda, que o estimador θ̂ é uma função dos valores
amostrais, ou seja, θ̂ = g(Y1 , Y2 , . . ., Yn ), é dado por:
1. atribuir massa de probabilidade 1/n para cada observação amostral Yj , j = 1, 2, . . ., n, criando a
distribuição de probabilidade empírica fˆ;
2. gerar amostras aleatórias com reposição da distribuição de probabilidade empírica fˆ, denominada
de amostra de bootstrap e dada por Ỹ1 , Ỹ2 , . . ., Ỹn ;
3. calcular uma estimativa de θ̂ por θ̃, usando a amostra de bootstrap no lugar da amostra original, da
seguinte forma θ̃ = g(Ỹ1 , Ỹ2 , . . ., Ỹn );
4. repetir os passos 2 e 3 B vezes.
O número de reamostragens foi denotado por B e a escolha de valores adequados conforme a situação
será discutida posteriormente. A inferência que pretendemos realizar depende da distribuição descon-
hecida de θ̂ − θ. Quando aplicamos o algoritmo anterior, obtemos uma amostra Monte Carlo, da
distribuição de bootstrap, da quantidade θ̃ − θ̂. Se n for suficientemente grande, pelas ideias do método
bootstrap, extremamente comprovadas e aceitas atualmente, esperamos que as duas distribuições sejam
aproximadamente idênticas.
8.3 Estimação
Para estimar parâmetros de uma população, a técnica bootstrap é uma excelente alternativa, principalmente
se a distribuição não for a normal, ou outra distribuição, cujos processos de estimação sejam bem
fundamentados teoricamente, como, por exemplo, a estimação de proporções binomiais p. Se o método
bootstrap for aplicado a esses casos, não será um problema, pois os resultados serão equivalentes ao da
teoria clássica.
No entanto, se as suposições de um processo de estimação clássico forem violadas, em geral, os métodos
bootstrap resultarão em processos de estimação com melhores propriedades. Inicialmente, devemos
entender que os métodos bootstrap não melhoram as estimativas pontuais, mas fornecem mecanismos
apropriados para estimar os erros padrões, intervalos de confiança e a distribuição de amostragem do
estimador, por mais complexo que seja esse estimador.
parâmetro θ, que pretendemos estimar, então aplicamos o algoritmo anterior e obtemos uma amostra
bootstrap de tamanho B da distribuição desse estimador.
Se a amostra for θ̃1 , θ̃2 , θ̃3 , . . ., θ̃B , podemos estimar o erro padrão de θ̂ por
v
v u P 2
B B B
u 1 X u 1 X θ̃
u u
θ̃j − θ̃¯ = u
2 j=1 j
Sθ̂ = t θ̃j2 − , (8.1)
B − 1 j=1 B − 1 j=1 B
t
em que
B
X
θ̃j
θ̃¯ =
j=1
.
B
Depreendemos da expressão (Equation 8.1), que o erro padrão do estimador de θ é o desvio padrão da
distribuição bootstrap desse estimador. Se, entretanto, o parâmetro θ for a média µ da distribuição de
probabilidade amostrada, podemos simplificar a expressão do estimador do erro padrão da média amostral
Ȳ , pois, nesse caso, não precisamos da amostra bootstrap do estimador.
Da teoria estatística, sabemos que o erro padrão da média amostral Ȳ é dado por SȲ = (1/n)1/2 S, sendo
S dado por v
u P 2
n
v
n Y
u 1 X u 1 X
u u
j=1 j
S=t (Yj − Ȳ )2 = u Yj2 − . (8.2)
n − 1 j=1 n − 1 j=1 n
t
import numpy as np
import [Link] as plt
import random
# função para retornar a mediana
def est_med(x):
y = [Link](x)
n = len(x)
if n%2 == 0:
est = (y[n//2 - 1] + y[(n+2)//2 - 1]) / 2
else:
est = y[(n+1) // 2 - 1]
return est
return estm
300
250
200
150
100
50
0
0.2 0.1 0.0 0.1 0.2
Como na distribuição de bootstrap de um estimador θ̂, a média representa a esperança em (Equation 8.3)
e o estimador θ̂ na amostra original, faz o mesmo papel de θ na população, podemos estimar o viés por
É comum pensarmos que θ̃¯ seria o estimador corrigido para viés, o que não é verdade. Podemos obter um
estimador ajustado por
θ̂aj. =θ̂ − δ̃θ̂ = θ̂ − θ̃¯ + θ̂
=2θ̂ − θ̃¯;
Na sequência, abordaremos alguns dos diferentes tipos de intervalos de confiança usualmente empregados
com o método denominado bootstrap não-paramétrico.
# função para obter o erro padrão do estimador
# e aplicar a correção de viés
def stderr_vies_est(x, est, B = 2000, *args):
result = dist_boot(x, est, B, *args)
se_boot = [Link](result)
est_boot = est(x, *args)
mean_boot = [Link](result)
vies_boot = mean_boot - est_boot
est_nvies = est_boot - vies_boot
return {"Estimativa": est_boot, "Erro padrão": se_boot, \
"Viés ": vies_boot, "Est. Não Viesada ": est_nvies}
134 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
Estimativa: 6.462984935204323
Erro padrão: 0.4599467403095527
Viés : -0.07137307722057162
Est. Não Viesada : 6.534358012424895
np.float64(6.931471805599453)
350
300
250
200
150
100
50
0
5.0 5.5 6.0 6.5 7.0 7.5 8.0
A principal propriedade requerida de um intervalo de confiança é que ele possua probabilidade de cobertura
igual ao valor nominal de confiança 1 − α adotado, o que, em geral, não é atendida nesse caso. A ideia é
admitir o efeito do teorema do limite central:
˙
θ̃∼N θ̂, Sθ̂2 , (8.5)
em que, onde se encontra o símbolo ∼,˙ deve se ler “se distribui aproximadamente como”, θ̂ é o valor do
estimador na amostra original e Sθ̂ é o estimador do erro padrão de θ̂, obtido na distribuição bootstrap.
Em consequência de (Equation 8.5}), podemos construir o intervalo bootstrap padrão com confiança de
aproximadamente 1 − α por
h i
IC1−α (θ) : θ̂ − Z1−α/2 Sθ̂ ; θ̂ − Zα/2 Sθ̂ , (8.6)
em que Zα/2 e Z1−α/2 são os quantis da distribuição normal padrão. O intervalo de confiança (Equation 8.6),
como já havíamos comentado, possui probabilidade de cobertura, em geral, inferior ao valor nominal de
100(1 − α)%.
Podemos melhorar este intervalo obtendo os limites do intervalo de confiança usando um estimador
corrigido para viés. . O resultado para este caso é:
h i
IC1−α (θ) : θ̂aj. − Z1−α/2 Sθ̂ ; θ̂aj. − Zα/2 Sθ̂ , (8.7)
estimativa: 9.292005812796827
erro padrão: 2.767882006829041
viés: -0.7257103215318583
136 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
estimativa: 0.2450601001212278
erro padrão: 0.0750243886992481
viés: 0.014427725897615207
limite inferior: 0.08358727441095248
limite superior: 0.3776774740362727
np.float64(0.227468211559786)
dos percentis 100(α/2)% e 100(1 − α/2)% obtidos diretamente da distribuição de bootstrap na escala
transformada.
Isso pode ser feito, em virtude do que acabamos de comentar, ou seja, que a distribuição de bootstrap
obtida pela transformação, exceto pelo fato de ser uma amostra finita, é normal e os dois intervalos na
escala transformada se equivalem assintoticamente.
Posto dessa forma, parece ser fácil obter um intervalo como o apresentado em (Equation 8.8). Isso não
é verdade, pois necessitaríamos obter a função ϕ, que levaria a uma distribuição de bootstrap normal.
Encontrar tal função para cada caso real é uma tarefa inexequível. Por outro lado, não precisamos conhecer
tal função, basta pressupor sua existência.
Como os valores do intervalo, obtidos na escala transformada, quando mapeados na escala original são
equivalentes aos percentis 100(α/2)% e 100(1−α/2)% da distribuição de bootstrap de θ̃, então a construção
do intervalo de confiança por esse método é trivial. Esse resultado é extremamente interessante e foi
preponderante para o sucesso dos métodos de estimação bootstrap. Como afirmam Efron e Tibshirani
(1993), podemos entender o método dos percentis para a obtenção do intervalo de confiança, como sendo
um algoritmo automático de incorporação de tais transformações que levam à normalidade. Para que fique
mais claro, podemos dizer que os limites obtidos de (Equation 8.8) correspondem exatamente ao intervalo
baseado em percentis da distribuição de bootstrap na escala original, portanto não precisamos conhecer a
transformação que irá conduzir a distribuição de bootstrap à normalidade, mas apenas pressupor sua
existência!
Para aplicar o intervalo de confiança baseado em percentis bootstrap, devemos aplicar o algoritmo
anteriormente apresentado. Obtemos: θ̃1 , θ̃2 , θ̃3 , . . ., θ̃B .
Devemos ordenar esses valores, obtendo as estatísticas de ordem por: θ̃(1) , θ̃(2) , θ̃(3) , . . ., θ̃(B) . Em seguida,
podemos obter o intervalo de confiança a partir dessa função de distribuição empírica, obtendo os percentis
100(α/2)% e 100(1 − α/2)%. Assim, construímos o intervalo de confiança baseado em percentis bootstrap
por
IC1−α (θ) : θ̃(k1 ) ; θ̃(k2 ) , (8.9)
em que k1 = ⌊(B + 1)(α/2)⌋ e k2 = ⌊(B + 1)(1 − α/2)⌋, ou seja, são os maiores inteiros que não são
maiores que (B + 1)(α/2) e (B + 1)(1 − α/2), respectivamente; θ̃(k1 ) é o percentil 100(α/2)% da função
de distribuição bootstrap empírica; e θ̃(k2 ) , o percentil 100(1 − α/2)% da função de distribuição bootstrap
empírica.
Esse intervalo possui propriedade melhores que o anterior e funciona bem na maioria dos casos mais
simples. Entretanto, existem alternativas melhores. Veremos, a seguir, alguns outros tipos de intervalos
que são melhores que os dois primeiros, embora em algumas situações o método percentil é superior a
alguns deles.
Uma grande vantagem do intervalo percentílico é que ele possui a propriedade de respeitar a transformação,
ou seja, os limites desse intervalo obedecem ao domínio do parâmetro que está sendo estimado, o que, em
alguns métodos, não ocorre.
# função para obter o IC Percentis
def ic_perc_boot(x, est, B = 2000, alpha = 0.05, *args):
res = dist_boot(x, est, B, *args)
est_boot = est(x, *args)
se_boot = [Link](res)
mean_boot = [Link](res)
vies_boot = mean_boot - est_boot
li, ls = [Link](res, [alpha/2 * 100, (1-alpha/2) * 100])
return {'estimativa': est_boot, 'erro padrão': se_boot, \
'viés': vies_boot, 'limite inferior': li, 'limite superior': ls}
138 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
estimativa: 0.2450601001212278
erro padrão: 0.0743344694305294
viés: 0.015100460649265346
limite inferior: 0.12254957974793343
limite superior: 0.39111047902372975
np.float64(0.227468211559786)
Como
d
θ̃ − θ̂ −
→ θ̂ − θ,
então, substituindo na expressão anterior, temos
h i
P θ̃(α/2) − θ̂ ≤ θ̂ − θ ≤ θ̃(1−α/2) − θ̂ = 1 − α.
Se isolarmos θ teremos h i
P 2θ̂ − θ̃(1−α/2) ≤ θ ≤ 2θ̂ − θ̃(α/2) = 1 − α. (8.10)
Os limites da afirmativa probabilística em (Equation 8.10) são os limites de uma intervalo de 100(1 −
α)% de confiança de bootstrap básico, dado por
h i
IC1−α (θ) : 2θ̂ − θ̃(1−α/2) ; 2θ̂ − θ̃(α/2) . (8.11)
Esse intervalo pode também ser visto como o intervalo de confiança t de bootstrap, que será o próximo a ser
estudado, quando o erro padrão em cada reamostragem for considerado igual a unidade. As quantidades
θ̃(1−α/2) e θ̃(α/2) do intervalo (Equation 8.11) são, respectivamente, os percentis θ̃(k1 ) e θ̃(k2 ) do intervalo
(Equation 8.9) baseado em percentis. Observamos que o presente método é apenas uma variação do
método percentil.
8.3. ESTIMAÇÃO 139
estimativa: 0.2450601001212278
erro padrão: 0.07622856459718486
viés: 0.010441124683899128
limite inferior: 0.09900972121872587
limite superior: 0.36865770046085516
np.float64(0.227468211559786)
em que k1 = ⌊(B + 1)p1 ⌋ e k2 = ⌊(B + 1)p2 ⌋, ou seja, são os maiores inteiros que não são maiores que
(B + 1)p1 e (B + 1)p2 , respectivamente; θ̃(k1 ) é o percentil 100(p1 )% da função de distribuição bootstrap
empírica; e θ̃(k2 ) , o percentil 100p2 % da função de distribuição bootstrap empírica.
Para valores conhecidos de â e ẑ0 , os valores de p1 e p2 , usados para determinar as ordens k1 e k2 do valor
de θ̃ na distribuição empírica de bootstrap, conforme apresentado em (Equation 8.12), são determinados
140 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
por ! !
ẑ0 + zα/2 ẑ0 + z1−α/2
p1 = Φ ẑ0 + e p2 = Φ ẑ0 + , (8.13)
1 − â ẑ0 + zα/2 1 − â ẑ0 + z1−α/2
em que Φ(x) é a função de distribuição da normal padrão avaliada no valor x e zα/2 e z1−α/2 são os
quantis 100(α/2)% e 100(1 − α/2)% dessa mesma distribuição, respectivamente.
Se tanto â quanto ẑ0 forem nulos, os valores de p1 e p2 serão, respectivamente, α/2 e 1 − α/2, como pode
ser facilmente deduzido de (Equation 8.13). Nesse caso, os intervalos BCa , expressão (Equation 8.12) e
bootstrap percentil, expressão (Equation 8.9), são equivalentes. Valores não nulos dessas quantidades,
modificam os limites do intervalo (Equation 8.12) e reduzem algumas deficiências do intervalo percentílico
(Equation 8.9), conforme afirmam Efron e Tibshirani (1993).
O valor de ẑ0 é definido como o quantil da normal padrão cuja probabilidade acumulada é dada pela
proporção de estimativas na distribuição de bootstrap que é menor que a estimativa de θ na amostra
original. Esse valor reflete a discrepância que existe entre a mediana da distribuição de θ̃ e a mediana da
distribuição de θ̂, expressa em unidades normais padrão.
Assim, definimos ẑ0 por
XB
I θ̃ i < θ̂
−1 i=1
ẑ0 = Φ (8.14)
,
B
em que Φ−1 (p) é a inversa da função de distribuição da normal padrão avaliada em p entre 0 e 1 e I(·) é a
função indicadora que retorna 1 se o valor de seu argumento for verdadeiro e 0, se for falso.
O valor ẑ0 será zero somente se a proporção de valores da distribuição de bootstrap que são inferiores
a θ̂ for igual a 50%. O valor de â pode ser calculado utilizando estimativas jackknife do parâmetro de
interesse.
Admitamos que o estimador θ̂ seja obtido por uma função do vetor de observações amostrais Y = [Y1 ,
Y2 , . . ., Yn ]⊤ por θ̂ = g(Y). Se eliminarmos a i-ésima observação Yi teremos o vetor Y(i) = [Y1 , Y2 , . . .,
Yi−1 , Yi+1 , . . ., Yn ]⊤ , que corresponde ao vetor original sem a i-ésima observação. O estimador jackknife
é obtido aplicando-se essa mesma função ao vetor resultante sem a i-ésima observação, sendo dado θ̂(i) =
g(Y(i) ).
Se calcularmos a média das n estimativas jackknife, obtidas com a eliminação de cada uma das n
observações originais, por
Xn
θ̂(i)
¯ i=1
θ̂(.) = ,
n
a aceleração â pode ser obtida por
n
¯
X 3
θ̂(.) − θ̂(i)
i=1
â = # . (8.15)
2 3/2
" n
X ¯
6 θ̂(.) − θ̂(i)
i=1
estimativa: 0.2450601001212278
erro padrão: 0.07390301412647568
viés: 0.011660164735816492
limite inferior: 0.12158628771838918
limite superior: 0.38694808161865485
np.float64(0.227468211559786)
142 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
sendo ẑ0 obtido da mesma forma que foi descrito anteriormente, por intermédio da equação (Equation 8.14).
Esse método é particularmente interessante em algumas situações em que os valores de â são difíceis de
ser obtidos. Esses casos são discutidos em Efron e Tibshirani (1993). Entretanto, o método BCa é, em
geral, mais eficiente e deve ser preferido.
# função para obter o IC com correção de viés
def ic_cv_boot(x, est, B = 2000, alpha = 0.05, *args):
res = dist_boot(x, est, B, *args)
est_boot = est(x, *args)
p = [Link](res < est_boot) / B
z0 = [Link](p)
z1 = [Link](alpha / 2)
z2 = [Link](1 - alpha / 2)
aux = 2 * z0 + z1
p1 = [Link](aux)
aux = 2 * z0 + z2
p2 = [Link](aux)
se_boot = [Link](res)
mean_boot = [Link](res)
vies_boot = mean_boot - est_boot
li, ls = [Link](res, [p1 * 100, p2 * 100])
return {'estimativa': est_boot, 'erro padrão': se_boot, \
'viés': vies_boot, 'limite inferior': li, 'limite superior': ls}
estimativa: 0.2450601001212278
erro padrão: 0.07638326075105437
viés: 0.01574564116660082
limite inferior: 0.10992136987364592
limite superior: 0.38694808161865485
np.float64(0.227468211559786)
# Demonstrando o [Link]
# Exemplo para a mediana da gama
B = 2000
alpha = 0.05
x = (x,) # dados devem estar numa tupla de array
res_med = [Link](x, est_med, n_resamples=B, \
confidence_level=1-alpha, method='percentile')
print(res_med.confidence_interval)
res_med = [Link](x, est_med, n_resamples=B, \
method='basic', confidence_level=1-alpha)
print(res_med.confidence_interval)
res_med = [Link](x, est_med, n_resamples=B, \
method='BCa', confidence_level=1-alpha)
print(res_med.confidence_interval)
# mediana teórica
[Link](0.5, a)
ConfidenceInterval(low=np.float64(0.12254957974793343), high=np.float64(0.39111047902372975))
ConfidenceInterval(low=np.float64(0.09900972121872587), high=np.float64(0.36757062049452216))
ConfidenceInterval(low=np.float64(0.12146249978160047), high=np.float64(0.3952728764288047))
np.float64(0.227468211559786)
144 CHAPTER 8. MÉTODOS BOOTSTRAP EM PYTHON
References
A. C. Atkinson and M. C. Pearce. the computer generation of beta, gamma and normal random variable.
journal of the royal statistical society, series a, 139(4):431–461, 1976.
J. N. W. Dachs. Estatìstica computacional: uma introdução em turbo Pascal. LTC, Rio de Janeiro, 1988.
L. Devroy. generating the maximum of independent identically distributed random variables. Computers
and mathematics with applications, 6:305–315, 1980.
L. Devroy. Non-uniform random variate generation. springer-verlag, new york, 1986.
D. F. Ferreira. Estatística multivariada. Editora UFLA, Lavras, 3 edition, 2018.
R. A. Johnson and D. W. Wichern. Applied multivariate statistical analysis. Prentice Hall, New Jersey, 4
edition, 1998.
V. Kachitvichyanukul and B. W. Schmeiser. binomial random variate generation. Communications of the
ACM, 31(2):216–222, 1988.
Donald E. Knuth. Literate programming. Comput. J., 27(2):97––111, May 1984.
G. Marsaglia and T. A. Bray. A convenient method for generating normal variables. Siam Review, 6(3):
260–264, 1964.
M. Matsumoto and T. Nishimura. Mersenne Twister: A 623-dimensionally equidistributed uniform
pseudo-random number generator. ACM Transactions On Modeling and Computer Simulation, 8(1):
3–30, 1998.
B. D. McCullough and B. Wilson. on the accuracy of statistical procedures in Microsoft Excel 97.
Computational Statistics and Data Analysis, 31:27–37, 1999.
T. H. J. Naylor, J. L. Balintfy, D. S. Burdick, and K. Chu. Técnicas de simulação em computadores.
Vozes, Petrópolis, 1971.
S. K. Park and K. W. Miller. Random number generators: good ones are hard to find. Communications
of the ACM, 31(10):1192–1201, 1988.
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical recipes in Fortran: the art
of scientific computing. Cambridge University Press, Cambridge, 1992.
A. Quarteroni, R. Sacco, and F. Saleri. Numerical mathematics. springer, Berlin, 2000.
L. Schrage. A more portable FORTRAN random number generator. ACM transactions on mathematical
software, 5(2):132–138, 1979.
W. B. Smith and R. R. Hocking. Algorithm AS 53: Wishart variate generator. Applied Statistics - Journal
of the Royal Statistical Society - Series C, 21(3):341–345, 1972.
D. H. D. West. Updating mean and variance estimates: an improved method. ACM transactions on
mathematical software, 22(9):532–535, 1979.
145
146 References
Index
algoritmo Wishart, 68
Hasting, 94 Wishart invertida, 69
amostragem funções
por rejeição, 46 trigonométricas, 50
beta gerador
incompleta, 89 padrão
mínimo, 39
classe geração
Python, 110 de variáveis
t multivariada, 72
densidade
exponencial, 48 Jacobiano
log-normal, 49 da transformação, 50
normal, 49
Tukey-lambda, 52 lema
distribuição soma de binomiais, 54
binomial, 54 tempo de espera
de combinações da exponencial, 55
lineares, 64 da geométrica, 55
t multivariada linguagem
esférica, 72 de alto-nível, 38
matriz
equação
covariâncias, 81
recursiva, 78
soma de quadrados e produtos, 81
estatísticas
Mersenne
de ordem, 137
Twister, 40
função método
de distribuição Box-Müller, 49
binomial, 87 congruencial, 37
exponencial, 85 da inversão
normal, 91 discreto, 56
de distribuição inversa métodos
binomial, 56 listas, 12
exponencial, 48, 85
números
de probabilidade
aleatórios, 37, 39
binomial, 54, 87
pseudo-aleatórios, 37
geométrica, 55 uniformes, 37
Poisson, 90
densidade precisão
exponencial, 85 dupla, 40
normal, 86
normal multivariada, 63 quadraturas
147
148 INDEX
gaussianas, 93
Monte Carlo, 94
quantis
exponencial, 85
random, 37
regra
trapezoidal
estendida, 91
simulação
matrizes aleatórias
Wishart, 70
Wishart invertidas, 70
soma
de produtos, 80
teorema
da transformação
de probabilidades, 45
tuplas, 14
Wishart
invertida, 69