Os Probabilistic Generative Models (PGM) são modelos de classificação para estimar a probabilidade de que uma determinada observação pertença a uma classe especifica. O modelo se utiliza do teorema de Bayes de modo que a partir de \( p(x|C)\) se obter \(p(C|x)\), ou seja, a partir da probabilidade de que se observe os valores \(x\) dado que a observação pertence a certa classe se obter a probabilidade de que uma outra observação pertença à class \(C\) a partir dos valores de suas variáveis. Trabalharemos aqui com apenas duas classes, de modo que o teorema de Bayes nos dá:
\[p(C_1|x) = \frac{p(x|C_1)p(C_1)}{p(x|C_1)p(C_1)+p(x|C_2)p(C_2)}\]
Desenvolveremos aqui o modelo seguindo de perto a seção 4.2 do livro "Pattern Recognition And Machine Learning" de Christopher Bishop.
Utilizaremos então a função Sigmoide de modo que:
\[ p(C_1|x) = \frac{1}{1+exp(-a)} = \sigma (a) \]
onde,
\[ a=ln\frac{p(x|C_1)p(C_1)}{p(x|C_2)p(C_2)} \]
Obteremos então \(a\) a partir da função logit, de modo que \(a=w^Tx+w_0\). Os coeficientes por sua vez, advém de:
\[ w = \Sigma^{-1}(\mu_1 - \mu_2) \]
\[w_0 = \frac{1}{2}(\mu_2\Sigma^{-1}\mu_2 - \mu_1\Sigma^{-1}\mu_1) + ln\frac{p(C_1)}{p(C_2)}\]
Onde, finalmente, \( \mu_1\) e \( \mu_2\) são as médias aritméticas das variáveis \( x\) pertencentes às classes 1 e 2, respectivamente, e \(\Sigma\) é a matriz de covariância de \(x\).
Passemos então ao código da implementação do PGM em Python.
Como parte do estudo envolve a aplicação do modelo à dados de Monte-Carlo, implementamos a função geradora dos dados.
import matplotlib.pyplot as plt
import numpy as np
def datagen(n,b0,b1,b2):
# Processo de geracao de dados
# semelhante ao dos slides da aula 23.
x1=np.random.normal(0,2,n)
x2=np.random.normal(10,4,n)
z=np.empty(n)
w=np.empty(n)
y=np.empty(n)
yr=np.empty(n)
for i in range(n):
z[i] = b0+b1*(x1[i])+b2*x2[i]
w[i] = np.random.normal(z[i],3)
for i in range(n):
if (w[i]>b0+10*b2) :
yr[i]=1
else:
yr[i]=0
if (z[i]>b0+10*b2):
y[i]=1
else:
y[i]=0
yr = np.array(yr)
return yr, y, x1, x2 # x1 e x2 sao as variaveis,
# y e o target e
# yr o target randomizado
Em seguida desenvolvemos as funções de geração da matriz de covariância e sigmoide.
def cov(a,b):
n1 = len(a)
n2 = len(b)
A = np.zeros((n1,n2))
for i in range(n1):
for j in range(n2):
A[i][j] = a[i]*b[j]
return A
def sgm(a): # Funcao sigmoide logistica de (a)
return 1/(1+np.exp(-a))
E finalmente à implementação do PGM.
def PGM(X,y): # Modelo gerador de probabilidades
# classificador binario.
# X[j][i] = variavel j observacao [i].
# Elaborado seguindo a secao 4.2
# do livro do Bishop.
# Na medida do possivel nomeei as variaveis
# da mesma forma que o autor.
n = len(X[0]) # Numero de observacoes
k = len(X) # Numero de variaveis
mu = np.zeros((2,k))
n1 = sum(y)
n0 = n-n1
pi = n1/n # p(C_1)
y0 = np.ones(n) - y
for j in range(k):
mu[1][j] = (np.dot(y, X[j]))/n1
mu[0][j] = (np.dot((y0), X[j]))/n0
S=np.zeros((k,k)) # Matriz de Covariancia (Sigma)
for i in range(n):
Si = X[:,i] - y[i]*mu[1] - y0[i]*mu[0]
S = S + cov(Si,Si)
S = S/n
Sv = np.linalg.inv(S) # Inversa de Sigma
# formas quadraticas de mu 0 e 1
mSvm0 = np.matmul(np.matmul(mu[0].T, Sv),mu[0])
mSvm1 = np.matmul(np.matmul(mu[1].T, Sv),mu[1])
w = np.matmul(Sv,mu[1] - mu[0])
w0 = (mSvm0 - mSvm1)/2 + np.log(pi/(1-pi))
ye = np.zeros(n) # Target estimado
yp = np.empty(n) # Probabilidade do target
a = np.empty(n)
for i in range(n):
a[i] = np.dot(X[:,i], w) + w0
yp[i] = sgm(a[i]) # p( y_i = 1 | x_i )
if yp[i] > 0.5 :
ye[i] = 1
erros = y - ye
R2 = 1 - sum(erros**2)/n
# Porcentagem de acertos, chamado aqui
# de R2 por conveniencia.
return w, w0, yp, ye , erros, R2
Em seguida precisamos estabelecer a função para testar o modelo encontrado sobre os dados de teste.
def Test(w,w0,X,y):
# Funcao para testar os valores encontrados
# no PGM no conjunto de teste.
n = len(X[0])
a = np.empty(n)
yp = np.zeros(n)
ye = np.zeros(n)
for i in range(n):
a[i] = np.dot(X[:,i], w) + w0
yp[i] = sgm(a[i]) # p( y_i = 1 | x_i )
if yp[i] > 0.5 :
ye[i] = 1
erros = y - ye
R2 = 1 - sum(erros**2)/n # Porcentagem de acertos
return yp, ye, erros, R2
Passemos agora para a aplicação do modelo aos dados de Monte-Carlo. Primeiro geraremos os dados usando a função já exposta. Em seguida treinaremos o modelo e depois geraremos novos dados para testar o algoritmo.
if __name__ == '__main__':
# Gerando dados por Monte Carlo para treinar o algoritmo
[yr,y,x1,x2] = datagen(300,10,-5,1)
X = np.array([x1,x2])
[w, w0, yp, ye, erros, R2] = PGM(X,yr)
print 'Montecarlo, treino, % de acertos: ', R2
# Gerando novos dados por Monte Carlo para testar o algoritmo,
# usando os mesmos parametros
[yr,y,x1,x2] = datagen(300,10,-5,1)
X = np.array([x1,x2])
[yp, ye, erros, R2] = Test(w,w0,X,yr)
print 'Montecarlo, teste, % de acertos: ', R2
# Plotando os resultados obtidos nos testes.
# Em azul e vermelho os valores atribuidos corretamente,
# em amarelo os falsos positivos,
# em verde 'cyan' os falsos negativos.
for i in range(300):
if(yr[i]==1 and ye[i] == 1):
plt.plot(x1[i],x2[i],u'ro',markersize=8, markeredgewidth=0)
if(yr[i]==1 and ye[i] == 0):
plt.plot(x1[i],x2[i],u'yo',markersize=10, markeredgewidth=0)
if(yr[i]==0 and ye[i] == 1):
plt.plot(x1[i],x2[i],u'co',markersize=10, markeredgewidth=0)
if(yr[i]==0 and ye[i] == 0):
plt.plot(x1[i],x2[i],u'bo',markersize=8, markeredgewidth=0)
plt.show
Como resultado obtemos: Montecarlo, treino, % de acertos: 0.906666666667, Montecarlo, teste, % de acertos: 0.9. E plotamos o gráfico com as classes estimadas e os erros de estimação separados em falsos positivos e falsos negativos.
Aplicamos então o modelo aos dados indicados no exercício 8.8 do livro "Modern Multivariate Statistical Techniques" de Alan Izenman. Os dados, chamados de The Insurance Company Benchmark foram obtidos no endereço: kdd.ics.uci.edu/databases/tic. A base possui 86 variáveis mais o target, ter o não comprado seguro. E está divida entre a amostra de treino, com 5.822 observações, e a amostra de teste, com 4.000 observações.
# Passando agora para o treino e teste com os dados
# reais obtidos em kdd.ics.uci.edu/databases/tic,
# ticdata2000.txt sao os dados para treino, ja com
# o target incluido. ticeval2000.txt sao os dados para teste
# sem target e tictgts sao os targets.
# Dados para treinar o algoritmo
ticdata = np.genfromtxt('ticdata2000.txt')
print np.shape(ticdata)
Xd = np.array(ticdata[:,0:85].T)
yd = np.array(ticdata[:,85])
[w, w0, yp, ye , erros, R2t] = PGM(Xd,yd)
print 'Dados reais, treino, % de acertos: ', R2t
# Dados para teste
ticdatatest = np.genfromtxt('ticeval2000.txt')
Xt = np.array(ticdatatest.T)
# Targets para teste
ticdatatarget = np.genfromtxt('tictgts2000.txt')
yt = np.array(ticdatatarget)
print 'shape Xt:', np.shape(Xt), 'shape yt:', np.shape(yt)
print 'Dados reais, teste, % de acertos', Test(w,w0,Xt,yt)[3]
Tendo a base uma quantidade muito grande de variáveis, não nos preocupamos em elaborar um gráfico, uma vez que não haveriam muitas informações que pudessem ser extraídas de uma visualização bidimensional. Como resultado obtivemos: "Dados reais, treino, % de acertos: 0.937306767434" e "Dados reais, teste, % de acertos 0.935".
A elevada taxa de acertos demonstra a robustez do PGM e a reduzida queda entre a taxa de acertos do treino e do teste mostra que o modelo é eficaz em absorver apenas o efeito real, sem incorporar o ruído dos dados em seu potencial explicativo.