In [19]:
# ignore esse código inicial, é apenas para preparar a página
from IPython.display import YouTubeVideo, HTML, Image, display
css_file = './modelo.css'
HTML(open(css_file, "r").read())
Out[19]:
Na matemática, uma Matriz é uma tabela de valores com $m$ linhas e $n$ colunas.
$$ \left( \begin{array}{cccc} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \end{array} \right) $$Um Vetor é uma Matriz com apenas $1$ coluna.
$$ \left( \begin{array}{c} 1 \\ 2 \\ 3 \\ \end{array} \right) $$Essas estruturas matemáticas possuem seu conjunto próprio de operações:
Com nosso conhecimento atual poderíamos representar vetores e matrizes como listas, porém tais operações não funcionariam.
In [1]:
u = [1,2,3]
v = [4,5,6]
print(u + v)
print(u + 1)
Alternativamente poderíamos implementar funções que realizassem tais operações:
In [2]:
def SomaVetores(u,v):
z = []
if len(u)!=len(v):
return z
for i in range(len(u)):
z.append( u[i] + v[i] )
return z
print (SomaVetores(u,v))
Mas a biblioteca numpy contém o tipo array que permite trabalhar com vetores e matrizes de maneira simplificada.
In [4]:
import numpy as np
u = np.array([1,2,3])
v = np.array([4,5,6])
print (u + v)
print (u + 1)
print (u*v )# produto elemento a elemento
print (np.dot(u,v)) # produto interno
print (np.outer(u,v) )# produto externo
Analogamente com matrizes:
In [4]:
A = np.array( [ [1,2], [3,4] ] )
B = np.array( [ [5,6], [7,8] ] )
print (A+B, '\n')
print (A*B, '\n' ) # multiplicação elemento a elemento
print (np.dot(A,B), '\n' )# multiplicação de matrizes
print (A.T, '\n') # transposta
print (np.linalg.inv(A), '\n' ) # inversa
Para acessar os elementos utilizamos pares de colchetes exatamente como nas listas.
Para matrizes colocamos os dois índices dentro dos colchetes.
In [5]:
print (u[0])
print (v[0:1])
print (A[0,0])
print (A[0,:])
Adicionalmente podemos passar uma lista de índices que queremos acessar:
In [6]:
u = np.array([i*i for i in range(10)])
print (u[ [1,3,5] ])
O tamanho de um vetor ou matriz pode ser obtido com o comando .shape.
In [7]:
print (A.shape)
Podemos criar uma Matriz com todos elementos iguais a $0$ ou $1$ ou até mesmo uma constante $c$.
In [8]:
m,n = 5,3
A = np.zeros((m,n)) # repare no parentes duplo
B = np.ones((m,n))
C = np.ones((m,n))*3
print (A, '\n')
print (B, '\n')
print (C, '\n')
A biblioteca numpy conta com todas as funções da biblioteca math e outras diversas funções úteis para cálculo numérico e álgebra.
In [9]:
print (np.cos(C), '\n')
print (np.arange(0,1,0.1), '\n')
In [10]:
D = np.random.random( (m,n) ) # matriz de números aleatórios de tamanho m x n
print (D)
In [11]:
print (D[ D<0.5]) # somente os elementos menores que 0.5
In [12]:
idx = np.where(D<0.5) # índices dos elementos < 0.5
print (idx, '\n')
print (D[idx])
In [13]:
print (D.max(), D.argmax()) # máximo valor e índice do máximo
print (D.min(), D.argmin()) # mínimo valor e índice do mínimo
In [14]:
print( D.mean(), D.std(), D.var()) # média e desvio-padrão e variância.
In [16]:
print (D.mean( axis=0 )) # média das colunas
print (D.mean(axis=1) ) # média das linhas
print (D.sum(axis=0)) # soma das colunas
In [17]:
# matriz aleatória gaussiana de média 1.7 e variância 0.5
D = np.random.normal( loc=1.7, scale=0.5, size=(m,n) )
print (D)
In [18]:
# random.choice() escolhe um elemento da lista dada a dist. de probabilidade p
print (np.random.choice( ['cara', 'coroa'], p=[0.5,0.5] )) # moeda honesta
print (np.random.choice( ['cara', 'coroa'], p=[0.8,0.2] )) # moeda desonesta
Vamos estimar a entropia de uma moeda honesta e de uma moeda desonesta utilizando Método de Monte Carlo:
In [19]:
def MonteCarlo( p, maxit ):
"""
Retorna a probabilidade de cara e coroa utilizando a função random.choice com probabilidade p
"""
caras = 0
for it in range(maxit):
moeda = np.random.choice( ['cara','coroa'], p=p )
if moeda == 'cara':
caras += 1
return np.array( [caras/float(maxit), 1.0 - caras/float(maxit)] )
def Entropia( p ):
"""
Calcula a entropia do vetor de probabilidades p
"""
return -(p*np.log2(p)).sum()
honesta = np.array([0.5, 0.5])
desonesta = np.array([0.8,0.2])
print( Entropia( MonteCarlo( honesta, 10000 ) ), Entropia( honesta ))
print (Entropia( MonteCarlo( desonesta, 10000 ) ) , Entropia( desonesta))
Simulação de Incêndio em Floresta
Dada uma floresta representada por uma matriz $m \times n$ em que os elementos são:
Utilizando a função random.choice() da biblioteca numpy, inicialize a matriz com $0$'s e $1$'s de tal forma que uma árvore é criada com probabilidade $p$ para criação de árvores.
Utilize a opção size para determinar o tamanho da matriz.
In [2]:
def CriaFloresta( m,n,p ):
"""
Cria uma floresta com 0s e 1s
"""
return np.random.choice( [0,1], size=(m,n), p=[1.0-p,p])
Crie uma função que faz o seguinte:
Podemos imaginar o código da seguinte maneira:
In [9]:
def SimulaIncendio( F ):
"""
Simula o incêndio em uma floresta
"""
FocoIncendio(F)
while temFogo(F):
AtualizaFloresta(F)
return F
Agora vamos fazer cada uma dessas funções em separado. Começando pela FocoIncendio(F).
In [11]:
def FocoIncendio(F):
# vamos pegar todos os pares (i,j) de índices contendo o valor 1
i,j = np.where(F==1)
# escolhe um deles para ser o foco de incêndio
quantos = len( i )
foco = np.random.choice( range(quantos) )
F[ i[foco], j[foco] ] = 2
Agora vamos implementar a função que verifica se ainda tem fogo na floresta:
In [12]:
def temFogo(F):
return len(np.where(F==2)[0]) > 0
Finalmente implementamos a função AtualizaFloresta(F):
In [13]:
def AtualizaFloresta(F):
# verifica os pontos de incêndio
I, J = np.where(F==2)
# para cada foco de incêndio
for i,j in zip(I, J):
EspalhaFogo(F, i, j)
ApagaFogo(F, i, j)
In [14]:
def EspalhaFogo(F, i, j):
# queremos verificar esses vizinhos de cada ponto de incêndio
vizinhos = [ (-1,-1), (-1,0), (-1,1), (0,-1), (0,1), (1,-1), (1,0), (1,1) ]
for k,l in vizinhos:
# se está dentro do limite e for árvore, bota fogo e avisa que botou fogo
if DentroFloresta(F, i+k, j+l) and temArvore(F, i+k, j+l):
F[i+k][j+l] = 2
In [15]:
def temArvore(F, i, j):
return F[i][j] == 1
def DentroFloresta(F, i, j):
return 0 <= i < F.shape[0] and 0 <= j < F.shape[1]
In [16]:
def ApagaFogo(F, i, j):
F[i,j] = 0 # apaga o fogo
In [17]:
F = CriaFloresta( 10,10, 0.3 )
F = SimulaIncendio(F)
print (F)
In [18]:
# Vamos plotar a floresta
%pylab inline
import matplotlib.pyplot as plt
plt.imshow(F, cmap="hot") # o comando imshow plota uma matriz de valores entre 0 e 1 com tonalidades de cores
plt.show();
Finalmente realizaremos um método de Monte Carlo para estimar o percentual de árvores restantes em um incêndio em uma floresta de densidade $p$.
Para cada densidade na lista [0.1, 0.2, ..., 0.9] e para cada repetição da simulação, crie uma floresta com essa densidade e simule o incêndio.
Calcule o percentual de árvores sobrevientes (média da matriz F) e calcule a média desse percentual em relação a todas as simulações.
In [19]:
def MonteCarlo( m, n, N ):
"""
Simulação de Monte Carlo para verificar a maior densidade em que um incêndio não causa muito estrago em uma floresta
"""
# queremos verificar as probabilidades de 0.1 até 0.9
Probabilidades = np.arange(0.1,1.0,0.1)
PF = []
# para cada probabilidade repete o processo de Monte Carlo
for p in Probabilidades:
media = Simula(p, m, n, N)
PF.append( (p, media) )
return PF
In [22]:
def Simula(p, m, n, N):
pf = 0
for it in range(N):
F = CriaFloresta( m,n,p )
F = SimulaIncendio(F)
pf += F.mean()
return pf/N
In [23]:
PF = MonteCarlo( 50,50,500 )
print ('Inicial \t Final')
for pf in PF:
print( '%.2f \t\t %.2f' % (pf[0], pf[1]))
Esse experimento mostra que em uma floresta com até 40% da área ocupada com árvores, um incêndio danifica apenas uma pequena região.
Iniciando em 50%, um incêndio elimina praticamente toda a floresta. Isso é conhecido como ponto crítico.
Para isso, vamos utilizar o arquivo dolar.csv contendo as informações do câmbio. Vejamos o formato desse arquivo.
In [36]:
!head dolar.csv
Cada linha representa uma das datas analisadas, a primeira coluna indica qual é data (sendo 0 a primeira data), a segunda coluna é o valor do Euro em relação ao Dólar e a terceira coluna o valor do Real em função do Dólar.
Notem que as colunas estão separadas por ';'
A leitura de um arquivo nesse formato para uma Matriz do numpy é feita com o comando genfromtxt.
Essa função pode receber diversos parâmetros, porém apenas dois são de nosso interesse:
genfromtxt( nome_do_arquivo, delimiter = 'delimitardor')
Nosso delimitador é o ';'
In [37]:
import numpy as np
data = np.genfromtxt('dolar.csv', delimiter=';')
print (data[:10,:]) # imprime os 10 primeiros registros
print (data.shape[0], "registros")
In [38]:
# As médias do US$/EURO e US$/R$ são:
print ("Média: ", data[:,1].mean(), data[:,2].mean())
# E o desvio-padrão:
print ("Desvio-padrão: ", data[:,1].std(), data[:,2].std())
In [44]:
# Vamos plotar a tendência de cada um
# O eixo x é a primeira coluna [0,1,2...] e o y serão a segunda e terceira colunas
plt.figure(figsize=(10,10))
plt.style.use('fivethirtyeight')
plt.title("Valorização do Dólar")
plt.plot(data[:,0], data[:,1], color='green', label="Euro")
plt.plot(data[:,0], data[:,2], color='blue', label="Real")
plt.legend()
plt.show()
In [45]:
polyEUR = np.polyfit( data[:,0], data[:,1], 5 )
polyR = np.polyfit( data[:,0], data[:,2], 5 )
print (polyEUR)
print( polyR)
In [57]:
polyEURf = np.poly1d(polyEUR)
polyRf = np.poly1d(polyR)
plt.figure(figsize=(10,10))
plt.style.use('fivethirtyeight')
plt.title("Valorização do Dólar")
plt.plot(data[:,0], polyEURf(data[:,0]), color='green', label="Euro fit", linewidth=5)
plt.plot(data[:,0], polyRf(data[:,0]), color='blue', label="Real fit", linewidth=5)
plt.plot(data[:,0], data[:,1], 'o', color='green', label="Euro", markersize=10, fillstyle='none')
plt.plot(data[:,0], data[:,2], 'o', color='blue', label="Real", markersize=10, fillstyle='none')
plt.legend()
plt.show()
In [58]:
# Vamos ver o que nos aguarda:
print ("Próximos dez dias do Real: ", polyRf(range(180,190)))
print()
print( "Próximos dez dias do Euro: ", polyEURf(range(180,190)))
Vamos salvar a projeção para os próximos 60 dias!
Para isso temos que criar uma matriz chamada proj de tamanho 60x3.
Na primeira coluna vamos inserir os valores de 180 até 240; na segunda coluna a projeção para o Euro e na terceira coluna a projeção para o Real.
Para salvar um arquivo csv usando o numpy utilizamos o comando savetxt:
savetxt( nome_arquivo, matriz, delimiter='delimitador', fmt=[ lista_formato ])
A lista de formatos é referente ao tipo de dado de cada coluna da matriz, no nosso caso a primeira coluna é inteiro ( %d) e as duas colunas seguintes são floats padronizados com 4 casas decimais ( %.4f).
In [59]:
# Vamos gravar as projeções:
proj = np.zeros( (60,3) )
proj[:,0] = range(180,180+60)
proj[:,1] = polyEURf(proj[:,0])
proj[:,2] = polyRf(proj[:,0])
np.savetxt('dolarProj60.csv', proj, delimiter=';', fmt=['%d', '%.4f', '%.4f'])
In [60]:
!head dolarProj60.csv