Primeiramente, vamos criar uma série temporal sintética. É um vetor de quantidade de chuva (milímetros) em cada dia do ano. Para modelar, vamos supor que a quantidade de chuvas é uma onda senoidal ao longo do ano.
Uma onda senoidal pode ser definida como $y(t) = H+A~sin(2\pi f~t + \phi)$, com H sendo a altura, ou o valor do eixo y onde o seno é zero (múltiplos de $\pi$), A é a amplitude da onda, ou a altura máxima que a onda chega (em relação à altura H), f é a frequência, ou o número de ondas por ciclo, e por fim $\phi$ é a fase, ou o deslocamento da onda.
In [42]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pandas as pd
from sklearn.model_selection import cross_val_score
from sklearn.naive_bayes import GaussianNB
from sklearn import tree
from sklearn import linear_model
from scipy.optimize import leastsq
In [28]:
def onda_senoidal(altura, amplitude, frequencia, fase, passos):
t = np.arange(passos)
return altura + amplitude * np.sin(2 * np.pi * frequencia * t + fase)
x = np.arange(100)
y = onda_senoidal(altura=0, amplitude=1, frequencia=0.03, fase=0, passos=100)
plt.plot(x, y);
Para modelar nosso exemplo da chuva, pensemos que a onda pode variar de 0 a 300 milímetros de chuva. Portanto, utilizaremos uma amplitude de 150 e uma altura de 150. A frequência é de uma onda a cada ano, portanto $f = \frac{1}{365}$. Modelaremos ainda o mês de janeiro, ou mais especificamente o dia 1 de janeiro, como o pico das chuvas, e o mês de junho como o dia mais baixo de chuva; para isso, a defasagem deve ser de um quarto do tamanho da onda, multiplicado pela frequência angular, logo, $\frac{365}{4}\omega = \frac{365}{4}\frac{1}{365}2\pi = \frac{\pi}{2}$.
O gráfico abaixo mostra a quantidade de chuva medida a cada dia, ao longo de 5 anos, com $t=0$ indicando o dia 1 de janeiro do ano 0, e $t=1825$ (o último dia) sendo 31 de dezembro do ano 4.
In [29]:
t = 365 * 5 # 5 anos
x = np.arange(t)
y = onda_senoidal(altura=150, amplitude=150, frequencia=1.0/365, fase=np.pi/2, passos=t)
plt.plot(x, y);
Para analisar algoritmos de classificação para este problema, devemos adicionar um pouco de ruído para gerar algumas dificuldades no aprendizado, facilitando a análise do viés indutivo. Adicionaremos uma variável aleatória $R \sim N(\mu=20, \sigma^2=30)$ na nossa série temporal. Precisaremos também truncar a série para ser positiva, não existe valores negativos de chuva!
In [30]:
n = len(y)
r = 30.0 * np.random.randn(n) + 20
y = y + r
y[y < 0] = 0
plt.plot(x, y);
Para nos retermos na classificação e não na tarefa de regressão, faremos a discretização da série acima. Para valores entre 0 e 80, consideraremos que houve "pouca ou nenhuma chuva". Para valores entre 80 e 200, consideraremos que houve "chuva leve", e para valores maiores que 200 consideraremos que houve "chuva forte".
In [31]:
def converter(v):
if 0 < v < 80:
return 0
elif 80 <= v < 200:
return 1
else:
return 2
c = [converter(i) for i in y]
print(c[70:90])
In [32]:
data = []
for i in range(101, len(c)):
data.append(c[i - 101 : i])
colunas = list(reversed(["1 dia antes"] + ["{0} dias antes".format(i) for i in range(2, 101)])) + ["Classe"]
df = pd.DataFrame(data, columns=colunas)
df.head()
Out[32]:
In [33]:
index = []
for i in range(len(df)):
index.append("Y{0}D{1}".format(i // 365, i % 365))
df["dia"] = index
df = df.set_index("dia")
df.head()
Out[33]:
In [34]:
nb = GaussianNB()
dt = tree.DecisionTreeClassifier()
ols = linear_model.LinearRegression()
lr = linear_model.LogisticRegression()
In [35]:
X = df.drop(['Classe'], axis=1)
Y = df['Classe']
In [36]:
for classificador, nome in [(nb, "Naive-Bayes"), (dt, "DecisionTree"), (ols, "Regressão Linear"), (lr, "Regressão Logística")]:
score = cross_val_score(classificador, X, Y, cv=10)
print "Acurácia do classificador {0}: {1}".format(nome, score.mean())
Como nós sabemos que o comportamento da chuva (no nosso exemplo) segue uma forma senoidal, podemos utilizar esse conhecimento de domínio para a regressão linear. Em vez de procurarmos a reta $Y = \alpha_0 + \alpha_1x_1 + ... + \alpha_nx_n$, podemos otimizar a curva $Y = \alpha_0 + \alpha_1~sin(\alpha_2x_1+\alpha_3)$, com os valores de $\alpha$ equivalendo aos parâmetros da onda senoidal.
Vamos também utilizar o conhecimento de domínio da frequência ser $\frac{1}{365}$ para facilitar o processamento do otimizador.
In [52]:
chute_altura = 100
chute_amplitude = 200
chute_fase = 0
alvo = lambda k: k[0] + k[1] * np.sin(2 * np.pi * (1.0/365) * x + k[2]) - y
parametros = leastsq(alvo, [chute_altura, chute_amplitude, chute_fase])[0]
print parametros
In [56]:
curva_ajustada = parametros[0] + parametros[1] * np.sin(2 * np.pi * (1.0/365) * x + parametros[2])
plt.plot(x, y, 'g', label='original')
plt.plot(x, curva_ajustada, 'r', label='regressao')
plt.legend()
plt.show()
Finalmente, utilizando a nova curva (bem) ajustada pela regressão linear, podemos fazer a mesma discretização que fizemos anteriormente no resultado da regressão para podermos comparar o resultado de sua acurácia.
In [58]:
predicao_curva = [converter(i) for i in curva_ajustada]
acertos = sum(1 for predicao, correto in zip(predicao_curva, c) if predicao == correto)
print 1.0 * acertos / len(c)