Escolhendo um valor de K para o Kmeans usando o critério de Calinski

No último post, aqui, a gente viu que quando agrupamos espécies em conjuntos, que tenham algum significado para a gente, existe coisas interessantes que podem ser feitas, como achar espécies que identifiquem aquele conjunto.

Mas isso é um tipo de classificação supervisionada, já que nos classificamos os locais de coleta, a priori, cada ponto está num grupo que nós definimos. Agora suponha que achemos que existem grupos para um conjunto de amostras, mas não sabemos a que grupo cada amostra está associada? E agora?

Existe um algorítimo chamado kmeans, que podemos usar. Mas o foco agora não é explicar como ele funciona, mas como escolher quantos grupos existem nos dados. Vamos tentar pensar num exemplo para ser mais claro. Pense que temos duas espécies só (ou se quiser, pense em dois atributos de uma flor, tamanho de pétalas e sépalas, similar aos dados do iris usado pelo Fisher), x e y, e podemos representá-las em um gráfico 2d.

01

Então temos quatro grupos, que podemos identificar visualmente na figura certo? Cada um dos grupos vem de distribuições normal com média e desvios diferentes, tanto para x quanto para y.

Note também que o preto e o verde estão um pouquinho mais embolados entre si (nos que quisemos assim) em relação ao vermelho e azul. Mas aqui essencialmente temos 120 amostras, e cada uma pertence a um dos 4 grupos que definimos no nosso exemplo, ao simular os dados.

Até aqui tudo bem, mas imagine, que não soubéssemos a cor de cada ponto, não soubéssemos a que grupo cada um pertence, e agora? Veja como é isso abaixo.

02

Aqui é um exemplo trivial, então só no olhômetro da para ver que existem grupos nessas 120 amostras, alias quatro grupos. Agora podemos usar o k-means para atribuir um grupo a cada uma dessas.

Então nos mandamos o kmeans classificar os dados em 4 grupos.

1
ajuste<-kmeans(exemplo,4)

Simples assim, e podemos comparar o resultado com os grupos originais.

03

Veja que a classificação erra um pouco, principalmente nos pontos verdes que estão misturados no pretos (lembrando que uma cor é um grupo, a gente troca as cores, porque o que interessa é pertencer ao grupo certo, não ter a cor certa), mas de modo geral, o kmeans recupera bem o grupo original de cada amostra. Claro que ele tem uma série de pressupostos a serem seguidos, mas isso fica para outro post.

Agora problema é que eu coloquei para separar em 4 grupos, porque eu quis assim. E isso é um problema, como a gente vai saber quantos grupos um conjunto de dados tem? Veja que assim como eu usei 4, eu poderia ter usado outras quantidades de grupos.

04

Lembrado, que quando realmente vamos usar o kmeans, a gente não sabe quantos grupos é a melhor opção. Ai entra o critério de Calinski, é uma medida de o quão bom é o ajuste de um determinado número de grupos, de forma que podemos descobrir quantos grupos usar. Podemos discernir qual dos quadros acima é mais adequado, usando o critério e todos e vendo quando ele é maximizado.

Ele está implementado no pacote vegan do R, na função cascadeKM, e é usado assim:

1
ajuste<-cascadeKM(exemplo,2,10,criterion = 'calinski')

A gente fala a nossa tabela de dados, o número mínimo de grupos, o número máximo e qual critério usar. Dai para cada valor de k(número de grupos), o cascadeKM vai rodar o kmeans, calcular o valor de calinski e dizer a qual grupo cada amostra esta nessa classificação.

Basicamente ele calcula o seguinte:

calinski=\frac{SSB/(k-1)}{SSW/(n-k)}

onde n é o número de amostras e k é o número de classes que estamos selecionando (que vai ser 2,3,4,…10, veja que colocamos 2 e 10 como argumentos). SSW é a soma dos quadrados dentro dos clusters e SSB é a soma dos quadrados entre clusters, sendo análogo a anova (que lembrem que F, quanto mais, mais quer dizer que existem grupo, veja aqui), ou seja, a gente vai procurar o maior valor de F.

Por padrão, o cascadeKM faz o seguinte plot

05

Que é a qual grupo cada ponto foi atribuído e uma figura com os valores de calinski e o número de k. Mas ele está invertido, então podemos também extrair os valores e produzir uma figura com o valor de k no eixo x e o resultado a métrica de calinski no eixo y.

06

Veja que o maior valor está exatamente no 4, que é como geramos os dados. Porém nem tudo é tão simples, esse critério vai se comportar melhor se os dados estiverem em uma região especifica do n-espaço dos dados, de preferencia com variação homogenia, além de outros casos. Mas de qualquer forma, ele pode ser informativo, sendo que podemos avaliar o resultado graficamente para ver se faz sentido, e ele é bastante simples de entender e calcular.

Bem é isso ai, como estava olhando essas coisas, resolvi aproveitar e deixar tudo que estava testando como post, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre, Peter R. Minchin, R. B. O’Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens and Helene Wagner (2016). vegan: Community Ecology Package. R package version 2.3-5. https://CRAN.R-project.org/package=vegan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
##Pacotes utilizados
library(vegan)
 
##Gerando um exemplo
x<-c(1,1,2,2)
y<-c(1,2,1,2)
 
set.seed(123)
exemplo<-data.frame(
    x=c(rnorm(30,x[1],0.3),rnorm(30,x[2],0.2),rnorm(30,x[3],0.3),rnorm(30,x[4],0.2)),
    y=c(rnorm(30,y[1],0.2),rnorm(30,y[2],0.2),rnorm(30,y[3],0.2),rnorm(30,y[4],0.2))
)
grupo<-factor(rep(letters[1:4],each=30))
 
##Exemplo gerado
plot(exemplo$x,exemplo$y,col=as.numeric(grupo),pch=19,frame=F,xlab="eixo x",ylab="eixo y")
legend("topright",legend=c(levels(grupo)),col=unique(as.numeric(grupo)),pch=19,bty="n")
 
##Dados não classificados
plot(exemplo$x,exemplo$y,pch=19,frame=F,xlab="eixo x",ylab="eixo y")
 
ajuste<-kmeans(exemplo,4)
 
##Resultado vc original
par(mfrow=c(2,1))
plot(exemplo$x,exemplo$y,col=as.numeric(grupo),pch=19,xlab="eixo x",ylab="eixo y")
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y")
 
##Comparando varios valores de k numa única figura
par(mfrow=c(2,2),mar=c(4,4,2,2))
ajuste<-kmeans(exemplo,2)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=2")
ajuste<-kmeans(exemplo,3)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=3")
ajuste<-kmeans(exemplo,6)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=6")
ajuste<-kmeans(exemplo,8)
plot(exemplo$x,exemplo$y,col=ajuste$cluster,pch=19,xlab="eixo x",ylab="eixo y",frame=F,main="K=8")
 
##
ajuste<-cascadeKM(exemplo,2,10,criterion = 'calinski')
 
##Plot original
plot(ajuste)
 
##Fazendo nosso próprio plot
plot(2:10,ajuste$results[2,],type="b",pch=19)
abline(v=which.max(ajuste$results[2,])+1,lty=3,lwd=3,col="red")

Regressão logística regularizada de mais graus como classificador

Então, lembram desse post aqui, onde a gente realizou a regressão logística como um classificador.

Então, tudo muito belo, muito lindo, mas vamos supor que os dados sejam assim:

01

Ai ferrou né? Se o mundo não for muito simples, e o cara que tire notas altas não seja aprovado, porque é meio maluco, mas o cara que não tira notas altas também não passa, e quem é médio, tira nota somente em umas das provas também não é aprovado, ai tudo fica pra la de complicado.

Usando a mesma ideia do post passado, não conseguiremos criar um bom classificador.

02

Veja que ele vai falhar nervosamente. Não conseguiremos classificar satisfatoriamente com esse modelo.

O problema é que veja como são os dados, eles tem uma estrutura complexa, parece que a área de aceite é um redondo no meio do plano desses dois preditores. Agora é a hora de usar um modelo complexo. Se você nunca entendeu porque as pessoas fazem aqueles modelos com preditores ao quadrado, ao cubo, interações. Bem são essas coisas que fazem os modelos fazerem curvas.

Como a gente estava fazendo um modelo aditivo, onde temos o primeiro preditor mais o segundo preditor e pronto, ficávamos com 3 parâmetros e um preditor que parecia uma reta, agora quando adicionamos complexidade a esses preditores, “versões” deles ao quadro, ao cubo, etc, esse modelo vai ter curvas, o que vai se adequar melhor a esse caso.

Bem o primeiro passo, é construir uma função, que chamaremos de mapFeature, que faz essas “versões” dos preditores ao quadrado, ao cubo, etc.

  maFeature(x) =\begin{bmatrix}  1 \\[0.3em]  x_1 \\[0.3em]  x_2 \\[0.3em]  x_1^2 \\[0.3em]  x_1 x_2 \\[0.3em]  x_2^2 \\[0.3em]  \vdots \\[0.3em]  x_1 x_2^5 \\[0.3em]  x_2^6  \end{bmatrix}

Então vamos fazer essa função, que fica assim

#Função para mapear a entrada em mais graus
mapFeature<- function(X1,X2,degree = 6) {
    out <-matrix(0,ncol=1,nrow=length(X1))
    out[,1]<-rep(1,length(X1))
     for(i in 1:degree) {
        for(j in 0:i) {
            out = cbind(out,(X1^(i-j))*(X2^j))
        }
    }
    return(out)
}

Veja que a função basicamente vai pegar dois preditores, X1 e X2 (ela não vai funcionar para mais preditores), e o grau que queremos transformar, e vamos pegar e fazer uma matriz, e na primeira coluna colocar ela inteira de 1, com o número de linhas igual o número de amostras, e a partir da segunda coluna, nos vamos basicamente resolver o binômio de newton. Veja que no começo, o j va valer zero e o i vale 1 menos zero, então basicamente é acrescentar o X1 como segunda coluna da matriz (lembre-se que qualquer coisa elevada a zero da 1, e qualquer coisa vezes 1 da ela mesmo), depois o i vale 1 e o j também, i-j da zero, dai acrescentamos X2, depois o i vai valer 2 e o j 0, então acrescentamos X1 ao quadrado, e assim vai, para o grau 6, da um total de 28 preditores, 28 parâmetros para estimar, e esses quadrados, cubos, vão fazer o modelo ficar com curvas.

Agora ao realizar o mesmo procedimento do post passado, obtemos o seguinte.

03

Melhor que uma reta, mais ainda assim, bem ruim. Isso porque para usar tantos parâmetros assim, precisamos penalizar o custo para esse tanto de parametros, precisamos de uma função de custo regularizada. Isso é um negocio duro de acreditar que funciona.

Direto do wikipedia, ele diz:

refers to a process of introducing additional information in order to solve an ill-posed problem or to prevent overfitting. This information is usually of the form of a penalty for complexity, such as restrictions for smoothness or bounds on the vector space norm.

Penalidade pela complexidade, então de alguma forma tem que diminuir o quanto os parâmetros mudam, diminuir o update dos parâmetros.

O duro é que parece algo do outro mundo, mas a implementação é relativamente simples.

Lembra como era o custo?
J(\theta)= \frac{1}{m} \sum\limits_{i=1}^m [-y^{(i)}log(h_\theta(x^{(i)}))-(1-y^{(i)})log(1-h_{\theta}(x^{(i)}))]

Então a gente vai adicionar essa penalidade, que vamos chamar de \lambda (lambda)
J(\theta)= \frac{1}{m} \sum\limits_{i=1}^m [-y^{(i)}log(h_\theta(x^{(i)}))-(1-y^{(i)})log(1-h_{\theta}(x^{(i)}))] + \frac{\lambda}{2m}\sum\limits_{j=1}^n \theta_j^2

E o update de parâmetros, que era assim:

\frac{\partial J(\theta)}{\partial \theta_J} = \frac{1}{m} \sum\limits_{i=1}^m (h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}_{j}

Vai ficar assim:

\frac{\partial J(\theta)}{\partial \theta_J} = \frac{1}{m} \sum\limits_{i=1}^m (h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}_{j} \ for\ j=0
\frac{\partial J(\theta)}{\partial \theta_J} = \frac{1}{m} \sum\limits_{i=1}^m (h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}_{j}+\frac{\lambda}{m}\theta_j \ for\ j\ge1

Agora como a gente so ta somando um componente a mais, um elemento a mais, a gente pode usar a função de custo anterior e so somar essa parte.

#Função de custo regularizado
costFunctionReg <-function(theta, X, y, lambda) {
    m <- length(y)
    J <- 0
    grad <- matrix(0,nrow=length(theta),ncol=1)
    saida<-costFunction(theta, X, y);
    theta<-saida$grad
    J<-saida$J
    J = J + ((lambda / (2*length(y))) * (t(theta)%*%theta))
    grad[1,1] = theta[1,1]
    grad[-1,1] = grad[-1,1] + ((lambda /length(y)) * theta[-1,1]) 
    return(list(J=c(J),grad=grad))
}

A gente pega usa a função anterior aqui

saida<-costFunction(theta, X, y);

Dai a gente tira o valor do custo nessa linha aqui

J<-saida$J

E a gente adiciona o \frac{\lambda}{2m}\sum\limits_{j=1}^n \theta_j^2

J = J + ((lambda / (2*length(y))) * (t(theta)%*%theta))

E basicamente é o mesmo esquema com os parâmetros, pegamos eles do saida

theta<-saida$grad

E fazemos o update, lembrando que aqui, para o primeiro parâmetro, o intercepto, ele não recebe a regularização, não sei explicar exatamente porque, mas é assim que funcionou.

    grad[1,1] = theta[1,1]
    grad[-1,1] = grad[-1,1] + ((lambda /length(y)) * theta[-1,1])

E mais um detalhe importante, agora a gente vai devolver os resultados como uma lista

return(list(J=c(J),grad=grad))

Isso porque, agora para fazer o update dos parâmetros no optim, a gente vai usar a nossa função, diferente do que fizemos no post anterior.
Agora fica assim

lambda <- 1
ajuste<-optim(rep(0,ncol(X)),method="L-BFGS-B",control = list(trace=T),
              fn =function(theta) {costFunctionReg(theta,X,y,lambda)$J},
              gr =function(theta) {costFunctionReg(theta,X,y,lambda)$grad})
ajuste$par

Primeiro a gente definiu o lambda, e no optim, a primeira coisa é que mudamos o algorítimo de otimização, isso porque o nelder-mead não vai permitir fornecermos os próximos valores dos parâmetros, o valores de theta, com esse ai, na verdade se você olhar no help da função, existem três algorítimos que vão permitir a gente usar o argumento gr, que é gr de gradient, ou seja, podemos fornecer uma função para calcular o próximo gradiente.

Eu num exatamente como o nelder-mead funciona, esse L-BFGS-B então não faço a menor ideia, mas é meio lógico, ele otimiza.

Mas agora podemos então usar esses parâmetros e ver como ficou o modelo.

04

Olha ai, com o \lambda=1 o modelo ficou bem bonzinho. Agora estamos pegando bem a área onde a maioria dos pontos aprovados estão, ou seja, nossa predição é bastante satisfatória, melhor que no modelo sem regularização, ou quando não usávamos mais graus nos modelo.

Mas o \lambda, o fator de regularização é complicado, aqui o exemplo é simples, mas ele precisa ter um valor razoável também, se for muito baixo, tipo \lambda=0.01 ficaremos com o seguinte

05

Um modelo com uma capacidade de predição ruim, e se ele for muito alto, por exemplo \lambda=1000

06

Ele pode acabar muito bom para esses dados, mas sofrer de overfitting, e se lembrar que queremos o melhor modelo para o sistema, e não para esse conjunto de dados, vemos que estamos fazendo merda.

Bem e acho que isso acabou esse exercício.
O legal de tentar implementar essas coisas é ver o tanto de coisa que está por traz de um simples modelo que a gente faz.

Por exemplo, o mapFeature é algo bem complicado de fazer, e eu tive que perder um tempo porque não consegui calcular antes de começar o modelo, o tamanho que a matriz tinha que ser, depois de um tempo eu abandonei isso e comecei a simplesmente grudar cada nova coluna com cbind na matriz. Esse trabalho no R é feito pela função formula, que consegue decompor o que a gente escreve depois do ~ e criar uma matriz como a que o mapFeature faz, mas de forma muita mais ampla, ja que aceita uma sintaxe muito maior.

Outra coisa é ver como os algorítimos de optimização atuam de forma diferentes, veja que se a gente mudar os algoritmo do L-BFGS-B pro BFGS por exemplo, parece que independente do lambda que a gente usa, o resultado sai muito mais similar. Eu não sei qual algorítimo na verdade o fminsc do matlab usa, porque meu resultado para lambda grande e pequeno não saiu idêntico, mas parecido.

Bem é isso ai, esse é um exercício do curso de Machine Learning do Coursera, deem uma olhada que é bem legal. Os dados e o script está no repositório do recologia, eu não coloquei aqui, mas da para abrir direto do repositório os dados.. E até o próximo post.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
#Função Sigmoid
sigmoid<- function(z) {
    return(1/(1+exp(-z)))
}
 
#Função de custo
costFunction<- function(theta, X, y) {
    m <- length(y)
    J <- 0
    grad <- matrix(0,nrow=length(theta)+1,ncol=1)
    h <- sigmoid(X %*% theta)
    costPos <- t(-y) %*% log(h)
    costNeg <- (1 - t(y)) %*% log(1 - h)
    J = (1/m) * (costPos - costNeg)
    grad = (1/m) * (t(X) %*% (h - y))
    return(list(J=c(J),grad=grad))
}
 
#Função de custo regularizado
costFunctionReg <-function(theta, X, y, lambda) {
    m <- length(y)
    J <- 0
    grad <- matrix(0,nrow=length(theta),ncol=1)
    saida<-costFunction(theta, X, y);
    theta<-saida$grad
    J<-saida$J
    J = J + ((lambda / (2*length(y))) * (t(theta)%*%theta))
    grad[1,1] = theta[1,1]
    grad[-1,1] = grad[-1,1] + ((lambda /length(y)) * theta[-1,1]) 
    return(list(J=c(J),grad=grad))
}
 
#Função para mapear a entrada em mais graus
mapFeature<- function(X1,X2,degree = 6) {
    out <-matrix(0,ncol=1,nrow=length(X1))
    out[,1]<-rep(1,length(X1))
     for(i in 1:degree) {
        for(j in 0:i) {
            out = cbind(out,(X1^(i-j))*(X2^j))
        }
    }
    return(out)
}
 
####################
exemplo2<-read.table("ex2data2.txt",sep=",",stringsAsFactors =FALSE)
exemplo2[,1]<-as.numeric(exemplo2[,1])
exemplo2[,2]<-as.numeric(exemplo2[,2])
colnames(exemplo2)<-c("QA1","QA2","Status")
 
#
plot(exemplo2[,1],exemplo2[,2],col=ifelse(exemplo2$Status==0,"red","green"),pch=19,xlab="Tesde de microship 1",
     ylab="Teste de microship 2",frame=F)
legend("topright",c("Aceito","Rejeitado"),pch=19,col=c("green","red"),bty="n")
 
#Parametros de entrada
X<-as.matrix(cbind(1,exemplo2[,1:2]))
theta<-matrix(0,nrow=ncol(exemplo2[,1:2])+1,ncol=1)
y<-exemplo2[,3]
 
ajuste<-optim(c(0,0,0),function(theta) {costFunction(theta,X,y)$J},method="Nelder-Mead",control = list(trace=T))
ajuste
 
 
#Figura 2
plot_x<-c(min(X[,2])-2,max(X[,2])+2)
plot_y<-(-1/ajuste$par[3])*(ajuste$par[2]*plot_x + ajuste$par[1])
plot(exemplo2[,1],exemplo2[,2],col=ifelse(exemplo2$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F)
lines(plot_x,plot_y,lty=3,lwd=4,col="black")
legend("topright",c("Adminitido","Reprovado","Margem de decisão"),pch=c(19,19,NA),lty=c(0,0,3),lwd=4,
       col=c("blue","red","black"),bty="n")
 
 
#################################
X<-mapFeature(exemplo2[,1],exemplo2[,2])
theta<-matrix(0,nrow=ncol(X),ncol=1)
y<-exemplo2[,3]
 
ajuste<-optim(rep(0,ncol(X)),function(theta) {costFunction(theta,X,y)$J},method="Nelder-Mead",control = list(trace=T))
ajuste$par
 
#Figura 3
plot(exemplo2[,1],exemplo2[,2],col=ifelse(exemplo2$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F)
u<-seq(-1, 1.5,length= 50);
v<-seq(-1, 1.5,length= 50);
z<-matrix(0,ncol=length(u),nrow=length(v))
for(i in 1:length(u)) {
    for(j in 1:length(v)) {
        z[i,j]<-mapFeature(u[i], v[j])%*%ajuste$par
    }
}
z = t(z)
contour(u,v,z,1,drawlabels=F,add=T,lwd=3,lty=2)
legend("topright",c("Adminitido","Reprovado","Margem de decisão"),pch=c(19,19,NA),lty=c(0,0,3),lwd=4,
       col=c("blue","red","black"),bty="n")
 
 
###############################
lambda <- 1
ajuste<-optim(rep(0,ncol(X)),method="L-BFGS-B",control = list(trace=T),
              fn =function(theta) {costFunctionReg(theta,X,y,lambda)$J},
              gr =function(theta) {costFunctionReg(theta,X,y,lambda)$grad})
ajuste$par
 
#Figura 4
#jpeg("04.jpg")
plot(exemplo2[,1],exemplo2[,2],col=ifelse(exemplo2$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F)
u<-seq(-1, 1.5,length= 50);
v<-seq(-1, 1.5,length= 50);
z<-matrix(0,ncol=length(u),nrow=length(v))
for(i in 1:length(u)) {
    for(j in 1:length(v)) {
        z[i,j]<-mapFeature(u[i], v[j])%*%ajuste$par
    }
}
z = t(z)
contour(u,v,z,1,drawlabels=F,add=T,lwd=3,lty=2)
legend("topright",c("Adminitido","Reprovado","Margem de decisão"),pch=c(19,19,NA),lty=c(0,0,3),lwd=4,
       col=c("blue","red","black"),bty="n")
#dev.off()
###############################
lambda <- 0.01
ajuste<-optim(rep(0,ncol(X)),method="L-BFGS-B",control = list(trace=T),
              fn =function(theta) {costFunctionReg(theta,X,y,lambda)$J},
              gr =function(theta) {costFunctionReg(theta,X,y,lambda)$grad})
ajuste$par
 
#Figura 5
#jpeg("05.jpg")
plot(exemplo2[,1],exemplo2[,2],col=ifelse(exemplo2$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F)
u<-seq(-1, 1.5,length= 50);
v<-seq(-1, 1.5,length= 50);
z<-matrix(0,ncol=length(u),nrow=length(v))
for(i in 1:length(u)) {
    for(j in 1:length(v)) {
        z[i,j]<-mapFeature(u[i], v[j])%*%ajuste$par
    }
}
z = t(z)
contour(u,v,z,1,drawlabels=F,add=T,lwd=3,lty=2)
legend("topright",c("Adminitido","Reprovado","Margem de decisão"),pch=c(19,19,NA),lty=c(0,0,3),lwd=4,
       col=c("blue","red","black"),bty="n")
#dev.off()
 
###############################
 
lambda <- 1000
ajuste<-optim(rep(0,ncol(X)),method="L-BFGS-B",control = list(trace=T),
              fn =function(theta) {costFunctionReg(theta,X,y,lambda)$J},
              gr =function(theta) {costFunctionReg(theta,X,y,lambda)$grad})
ajuste$par
 
#Figura 6
#jpeg("06.jpg")
plot(exemplo2[,1],exemplo2[,2],col=ifelse(exemplo2$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F)
u<-seq(-1, 1.5,length= 100);
v<-seq(-1, 1.5,length= 100);
z<-matrix(0,ncol=length(u),nrow=length(v))
for(i in 1:length(u)) {
    for(j in 1:length(v)) {
        z[i,j]<-mapFeature(u[i], v[j])%*%ajuste$par
    }
}
z = t(z)
contour(u,v,z,1,drawlabels=F,add=T,lwd=3,lty=2)
legend("topright",c("Adminitido","Reprovado","Margem de decisão"),pch=c(19,19,NA),lty=c(0,0,3),lwd=4,
       col=c("blue","red","black"),bty="n")
#dev.off()

Regressão Logística como um Classificador em aprendizado de máquina

Normalmente, sempre a gente estuda algo focada para algum uso.
Eu mesmo, mais do que saber teoria de estatística, eu quero usar para analisar dados em biologia, ver padrões, fazer afirmações. Mas olhando focado as coisas, as vezes a gente perde o todo.

Por exemplo na regressão logística, ou regressão binomial. Agente já falou dela aqui, aqui e um pouquinho aqui.

Pra mim essa era uma análise para ver relação entre atributos e prevalência e pronto. Vendo mais um pouco, vi que dava para pensar coisas em relação a chance de algum evento acontecer, mas como modelo de classificação, só mesmo depois de olhar sobre aprendizado de máquina.

Aqui tem parte de um dos exercícios da disciplina que fiz no coursera, sobre machine learning, muito legal ela, recomendo, o post sobre Gradient Descent também era exercício dessa disciplina, que é feita em matlab ou octave, mas é relativamente fácil adaptar as coisas para o R.

Mas vamos la, basicamente a gente vai construir uma modelo de regressão logística do zero, mais ou menos do zero, que a gente vai usar o optimizador do R (o gradient descent é um optimizador, de uma olhada no post sobre ela para visualizar o que o optim vai ta fazendo mais ou menos), mais especificamente, uma regressão para a chance de ser admitido em uma universidade dados suas notas em exames anteriores.
A gente tem então os dados históricos de candidatos as vagas (leia as notas deles em exames anteriores) e temos o resultado deles na universidade, se ele foram aceitos ou não. Aqui, cada candidato vai ser um exemplo de treinamento, e a nota deles nos exames anteriores vão ser os preditores da chance deles de serem admitidos.

Bom os dados vão estar no repositório do recologia no github, pode baixar de la o arquivo e abrir no R.
Dai vamos dar uma olhadinha nesses dados.

01

Bem, a princípio parece bem obvio, veja que quem tem notas boas nos exames, normalmente é admitido, isso é de certa forma trivial, mas veja que temos os dois cantos, mesmo caras que foram bem no exame 1 mas ruim no exame 2 foram admitidos. Agora a modelagem começa quando começamos a querer transformar bom em um número. Veja que há algo como uma fronteira que separa candidatos admitidos de não admitidos, só que ela não é exata. E isso que vamos buscar com a regressão logística aqui, os parâmetros da equação que normalmente olhamos para ver como é a relação entre as variáveis, aqui vai ser o valores da fronteira que vamos traçar entre os candidatos admitidos e não admitidos.

Bem então primeiro de tudo precisamos da nossa função de ligação (ela ja existe no R, é a função plogis, mas vamos implementar aqui para entender melhor).

Ela é o que vai “enrretar” nossas contas, lembre-se que binomial é probabilidade, que vai de 0 a 1, e a função sigmoide vai mudar isso para ir de menos infinito para mais infinito.

Vamos lembrar que a hipótese da regressão logística é:
h_{\theta}(x) = g(\theta^Tx)

E o g é a função sigmoid que é:

g(z) = \frac{1}{1+e^{-z}}

Veja que a função sigmoid da para passar para o R de forma direta praticamente.

sigmoid<- function(z) {
    return(1/(1+exp(-z)))
}

Um teste rápido é ver o valor para zero que tem que ser 0.5

> sigmoid(0) [1] 0.5

Mas vamos fazer um gráfico e ver que ela é a função que vai converter valores do conjunto do menos infinito até o mais infinito para o conjunto de 0 a 1.

02

Então a regressão, vai ser o de sempre uma hipótese, e a função de custo para os parâmetros dessa hipótese. A hipótese é o h_{\theta} ali de cima, sendo que o \theta vão ser os parâmetros da equação da reta, mas usando a função sigmoid para transformar em probabilidade (valores entre 0 e 1). E o custo de um set de parâmetros vai ser.

J(\theta)= \frac{1}{m} \sum\limits_{i=1}^m [-y^{(i)}log(h_\theta(x^{(i)}))-(1-y^{(i)})log(1-h_{\theta}(x^{(i)}))]

e o gradiente do custo é um vetor de mesmo tamanho de \theta onde o j ésimo elemento (for j=0,1, … , n) é definido como:

\frac{\partial J(\theta)}{\partial \theta_J} = \frac{1}{m} \sum\limits_{i=1}^m (h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}_{j}

A implementação disso no R fica assim.

costFunction<- function(theta, X, y) {
    m <- length(y)
    J <- 0
    grad <- matrix(0,nrow=length(theta)+1,ncol=1)
    h <- sigmoid(X %*% theta)
    costPos <- t(-y) %*% log(h)
    costNeg <- (1 - t(y)) %*% log(1 - h)
    J = (1/m) * (costPos - costNeg)
    grad = (1/m) * (t(X) %*% (h - y))
    return(c(J))
}

Bem alguns comentários, primeiro %*% no R é multiplicação de matriz, enquanto somente o * é multiplicação caso a caso, que nas matrizes vai ser elemento por elemento. A multiplicação de matriz resolve aquelas somatórios, mas tem que prestar atenção ali. O melhor jeito (que eu uso para ir construindo) é ir rodando a função linha por linha, e tentar entender como as operações estão sendo feitas, mas no fim das contas é a formula ali de cima.

No caso do optim do R, a gente so precisa calcular o J. O grad dar a derivada parcial dos parâmetros, isso é preciso no matlab, mas aqui o R cuida disso no optim. Agora no Matlab, e o Octave que é a versão open-source do Matlab, a gente usava a função fminunc para minimizar o custo J dado um conjunto de parâmetros\theta, e la tinha que retornar essa derivada parcial.

Agora para um set de parâmetros qualquer, a gente entra eles na função, mais os dados de entrada e saída, e ele fala custo de usar esse conjunto de parâmetros

> X<-as.matrix(cbind(1,exemplo1[,1:2])) > y<-exemplo1[,3] > theta<-matrix(0,nrow=ncol(exemplo1[,1:2])+1,ncol=1) > theta [,1] [1,] 0 [2,] 0 [3,] 0 > costFunction(theta,X,y) [1] 0.6931472

E o que a gente quer? Escolher o grupo de parâmetros, essa matriz de 3 números que vai gerar o menor número J, o menor custo, aquele que tem o maior likelihood de acontecer.

Certo, agora vamos usar o optim para milagrosamente conseguir chegar a esses valores

ajuste<-optim(c(0,0,0),function(theta) {costFunction(theta,X,y)},method="Nelder-Mead",control = list(trace=T))

Bem, um detalhe, é que veja que no optim tem como entrada o conjunto de parâmetros iniciais e a função que da o custo (verosimilhança), no caso ela precisa acessar os dados também, que estão nos vetores X (notas dos exames) e y (resultado, se o cara foi admitido ou não), mas o valores de \theta tem que ser a entrada da função, assim a gente cria uma função de uma única entrada, ali mesmo como argumento, sem dar um nome pra ela, algo como uma função anônima, que só existe dentro do optim, a ela tem como entrada o \theta e dentro dela ela usa essa entrada, mas chama também o X e o y, que como estão no environment, o central que fica tudo que a gente vai abrindo na sessão, o R olha la também para achar esse X e y.

Agora, esses são os argumentos principais, função e conjunto de parâmetros para começar, mas além disso, a gente definiu qual o método de optimização que vamos usar, o Nelder-Mead, que na verdade é o default, mas colocamos apenas para ver qual estamos usando, e para saber melhor o que esta acontecendo, a gente modificou um parâmetro do control, esse control é uma função, com configurações do algorítimo de optimização, o nelder-mead aqui, que contem um série de configurações, no trace=T estamos falando que queremos ver iteração por iteração o que esta acontecendo, que a função imprima no console qual a situação do custo, os valores dos parâmetros, essas coisas.

A agora é so pimba na gorduchinha, vamos rodar.

> ajuste<-optim(c(0,0,0),function(theta) {costFunction(theta,X,y)},method="Nelder-Mead",control = list(trace=T)) Nelder-Mead direct search function minimizer function value for initial parameters = 0.693147 Scaled convergence tolerance is 1.03287e-08 Stepsize computed as 0.100000 BUILD 4 2.189712 0.684397 LO-REDUCTION 6 2.087842 0.684397 HI-REDUCTION 8 1.105128 0.684397 HI-REDUCTION 10 0.848758 0.684397 . . . HI-REDUCTION 148 0.203498 0.203498 HI-REDUCTION 150 0.203498 0.203498 LO-REDUCTION 152 0.203498 0.203498 HI-REDUCTION 154 0.203498 0.203498 Exiting from Nelder Mead minimizer 156 function evaluations used

Bem, esse é um optimizador, a gente já viu o do Newton aqui, mas o Nelder Mead espero um dia fazendo um post sobre ele, ja o John Nelder e Robert Wedderburn que unificaram esse monte de análise que agente vê no glm. Mas então esses são os valores de controle da optimização, la em cima a gente vê a tolerância, esses outros argumentos. Mas agora que ajustamos os parâmetros, podemos ver como é o ajuste.

> ajuste $par [1] -25.1647048 0.2062524 0.2015046 $value [1] 0.2034977 $counts function gradient 156 NA $convergence [1] 0 $message NULL

Na saída temos os valores finais de parâmetros no $par, o J no $value e o $counts
é o total de iterações que tivemos que percorrer pra chegar aqui.

Agora se você me da as notas de um aluno, eu te falo a chance de ele ser aprovado

> Nota_Exame1<-80 > Nota_Exame2<-80 > sigmoid(ajuste$par[1]+ajuste$par[2]*Nota_Exame1+ajuste$par[3]*Nota_Exame2) [1] 0.9994223

Veja que a gente monta ai em cima a função linear da regressão logística, e trazemos pro valor de probabilidade com a função sigmoid.

Beleza, quem ia bem nas provas ta dentro. Se o cara tira 8 nos dois exames de entrada, tem 99% de chance de ser aprovado, poxa quem vai bem nas provas sempre vai bem no exame de admissão, mas o que é ir bem? Podemos ir testando valores agora nessa equação

> sigmoid(ajuste$par[1]+ajuste$par[2]*90+ajuste$par[3]*40) [1] 0.8112565 > sigmoid(ajuste$par[1]+ajuste$par[2]*70+ajuste$par[3]*50) [1] 0.3425826 > sigmoid(ajuste$par[1]+ajuste$par[2]*50+ajuste$par[3]*50) [1] 0.008352106

Veja que ir bem em uma prova e mal na outra, ainda garante boas chances, mas tirar uma nota mais ou menos nas duas provas, você já esta para traz, no segundo exemplo, tirando 7 e 5 nos exames, te deixa com 34% de chance de ser admitido, quem tirou 5 e 5 não precisa nem prestar exame de admissão que já era praticamente, só uma cagada muito grande para entrar. Mas veja que testando esse valores, a gente esta estabelecendo os limites de que de acordo com sua chance você vai ser admitido ou não, dependendo das suas notas. Mas se for para dar uma resposta, esse cara vai ser admitido ou não, a gente pode pensar que 50% seria o mínimo de chance para ser admitido.
Veja que no caso de tirar 90 num exame e 40 em ouro, a pessoa tem 81% de chance de ser admitido, mas do mesmo jeito ela tem 19% de chance de não ser admitida, então pra esse cara é melhor chutar que ele vai ser admitido, prevemos que ele vai ser admitido porque ele tem mais chance de ser admitido que não admitido.
Agora vamos pegar o caboclo que tirou 70 e 50, ele tem 34% de chance de ser admitido, logo ele tem 66% de chance de não ser admitido (100%-34%), então pra ele é melhor pensar que ele não vai ser admitido.
Então nessa lógica, a gente pode traçar uma linha no 50%, que é o limite para achar que alguém vai ser admitido, e do outro lado dessa linha a gente começa a chutar que a pessoa não vão ser admitidas, porque é o que tem mais chance de acontecer.

E isso fica assim.

03

E agora temos uma linha que representa o limite de ir bem ou mal nas notas, de um lado temos os cara que esperamos que serão aprovados, do outro o que esperamos que não serão, e dividimos o espaço das possibilidade de notas em área de pessoas que achamos que serão admitidas e área de pessoas que achamos que não serão admitidas. Claro que sempre tem a galera do contra, falhamos um pouco na classificação ali perto da linha, mas quanto mais longe dela, melhor classificamos as coisas.

Esses pontos são o que usamos para fazer o modelo, mas pense que podemos adicionar novos pontos agora, usando esses mesmo parâmetros.

Legal né, e por curiosidade, vamos ajustar esse modelo usando o glm do R para ver como ficam os parâmetros.

> ajuste$par [1] -25.1647048 0.2062524 0.2015046 > glm(Status~Exame1+Exame2,data=exemplo1,family="binomial") Call: glm(formula = Status ~ Exame1 + Exame2, family = "binomial", data = exemplo1) Coefficients: (Intercept) Exame1 Exame2 -25.1613 0.2062 0.2015 Degrees of Freedom: 99 Total (i.e. Null); 97 Residual Null Deviance: 134.6 Residual Deviance: 40.7 AIC: 46.7 >

E num é que da a mesma coisa 🙂
É tão legal quando as coisas funcionam.

Bem é isso ai, esse é um exercício do curso de Machine Learning do Coursera, deem uma olhada que é bem legal. Os dados e o script está no repositório do recologia. E até o próximo post.

exemplo1<-read.table("ex2data1.txt",sep=",",stringsAsFactors =FALSE)
exemplo1[,1]<-as.numeric(exemplo1[,1])
exemplo1[,2]<-as.numeric(exemplo1[,2])
colnames(exemplo1)<-c("Exame1","Exame2","Status")
 
#Figura 1
plot(exemplo1[,1],exemplo1[,2],col=ifelse(exemplo1$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F)
legend("topright",c("Adminitido","Reprovado"),pch=19,col=c("blue","red"),bty="n")
 
#Função Sigmoid
sigmoid<- function(z) {
    return(1/(1+exp(-z)))
}
 
sigmoid(0)
 
#Figura 2
plot(seq(-6,6,0.01),sigmoid(seq(-6,6,0.01)),type="l")
 
#Parametros de entrada
X<-as.matrix(cbind(1,exemplo1[,1:2]))
theta<-matrix(0,nrow=ncol(exemplo1[,1:2])+1,ncol=1)
y<-exemplo1[,3]
 
#Função de custo
costFunction<- function(theta, X, y) {
    m <- length(y)
    J <- 0
    grad <- matrix(0,nrow=length(theta)+1,ncol=1)
    h <- sigmoid(X %*% theta)
    costPos <- t(-y) %*% log(h)
    costNeg <- (1 - t(y)) %*% log(1 - h)
    J = (1/m) * (costPos - costNeg)
    grad = (1/m) * (t(X) %*% (h - y))
    return(c(J))
}
costFunction(theta,X,y)
 
#Ajuste da regressão
ajuste<-optim(c(0,0,0),function(theta) {costFunction(theta,X,y)},method="Nelder-Mead",control = list(trace=T))
ajuste
costFunction(matrix(ajuste$par,nrow=ncol(exemplo1[,1:2])+1,ncol=1),X,y)
 
#Previsões
Nota_Exame1<-80
Nota_Exame2<-80
sigmoid(ajuste$par[1]+ajuste$par[2]*Nota_Exame1+ajuste$par[3]*Nota_Exame2)
 
sigmoid(ajuste$par[1]+ajuste$par[2]*90+ajuste$par[3]*40)
sigmoid(ajuste$par[1]+ajuste$par[2]*70+ajuste$par[3]*50)
sigmoid(ajuste$par[1]+ajuste$par[2]*50+ajuste$par[3]*50)
 
#Figura 3
plot_x<-c(min(X[,2])-2,max(X[,2])+2)
plot_y<-(-1/ajuste$par[3])*(ajuste$par[2]*plot_x + ajuste$par[1])
plot(exemplo1[,1],exemplo1[,2],col=ifelse(exemplo1$Status==0,"red","blue"),pch=19,xlab="Exame 1",ylab="Exame 2",frame=F,
     ylim=c(30,110))
lines(plot_x,plot_y,lty=3,lwd=4,col="black")
legend("topright",c("Adminitido","Reprovado","Margem de decisão"),pch=c(19,19,NA),lty=c(0,0,3),lwd=4,
       col=c("blue","red","black"),bty="n")
 
ajuste$par
glm(Status~Exame1+Exame2,data=exemplo1,family="binomial")

Machine Learning – Classificação binária bayesiana

Recentemente eu comecei a ler um livro chamado Machine Learning for Hackers. Ele é bem legal, com o intuito de ensinar por exemplos.
Bem, segundo o wikipedia, Machine learning, é um ramo da inteligencia artificial que estuda sistemas que podem aprender a partir de dados, por exemplo um sistema que pode ser treinado a distinguir mensagens de e-mail, se estas são spam ou não. Após o treinamento, o sistema pode então ser usado para classificar novas mensagens se são spam ou não. Esse treinamente é algo como delimitar as fronteiras de o que é uma classe ou outra:

01

Bem legal né. Então eu decidi reproduzir o primeiro exemplo do livro, caso alguém queira dar uma olhada, todos os exemplos do livros estão num repositório de um dos autores aqui.

O primeiro exemplo explora os dados de mensagem de um projeto chamado SpamAssassin, um filtro de spam opensource, que da pra dar uma olhada nele nesse endereço aqui, mas no repositório do github, esta tudo pronto, os e-mails que serão usados, são diretórios nos quais cada e-mail esta num arquivo txt.

Esses são arquivos de texto, então pra transformar eles em números, a gente utiliza o pacote tm para text mining,ou mineração de texto.

Bem a idéia é criar um modelo com fronteiras não lineares entre o que é spam ou não.
Por exemplo a ocorrência da palavra html é um sinal que a mensagem é um spam, normalmente spam vai ter um link pra você clicar e sei la, instalar um vírus, ver uma propaganda indesejada, abrir alguma janela. Mas não existe uma relação simples, do tipo quanto mais vezes a palavra html aparecer, mais a mensagem é um spam. Existe muitas outras palavras na mensagem que podem diminuir nossa crença de que a mensagem é um spam.

Mas a solução é probabilidade condicional, que a gente ja viu aqui. Ou melhor se a gente vai apostar na soma do resultado de dois dados de 6 lados, e ja saiu um 6 no primeiro dado, ninguém vai querer apostar que a soma vai dar 1, 3 ou 5, ou qualquer valor menor que seis correto?

Para colocar essa idéia em prática, a gente usa o algorítimo conhecido como Naive Bayes classifier. De forma geral ele funciona assim, uma fruta é considerada uma maça se ela é vermelha, arredondada e tem uma tamanho x. A Naive Bayes classifier considera esses três quesitos independentemente para a probabilidade da fruta ser ou não uma maça, independentemente da outra. Ou seja uma maça verde ainda pode ser classificada como maça, dependendo do seu tamanho e forma, e algumas irregularidades na forma são aceitáveis. Mas ele não usa exatamente o teorema de bayes, ja que apesar de multiplicar por um prior, ele não divide por uma probabilidade, por isso a probabilidade vai ficando um número pequeno, já que só temos multiplicações.
Tudo que a gente precisa é um set inicial de dados já classificados para treinar, estabelecer limites para esses parâmetros e depois podemos usar o modelo resultante em novos dados.

No caso de e-mail, algumas palavras são batatas num e-mail de spam. Por exemplo palavras como html e table.

Lembre-se que quando você le nesse post algo assim:

Naive Bayes classifier.

Na verdade o texto esta em html, o computador esta lendo:

< a href="http://en.wikipedia.org/wiki/Naive_Bayes_classifier">Naive Bayes classifier.

Nos dados de teste por exemplo, a ocorrencia dessas palavras é mais o menos assim:

Term | Spam | Ham html | 377 | 9 table | 1,182 | 43

Então podemos usar a ocorrência dessas palavras para delimitar a fronteira do que é spam ou não.

01_init_plot1

Só que mesmo mensagens que não são spam, vira e mexe teremos o envio de links, então podemos usar a ocorrência de outras palavras para corrigir um possível erro de chamar de spam uma mensagem legítima, ou o contrario, que é chamar uma mensagem legítima de spam.

Bem a parte de mineração do texto, eu acho bem complicada, além de o pacote tm ter uma estrutura de dados própria para melhor manipular essas informações, isso precisa de uma boa leitura, então pulamos para depois.

Mas se você pegar todas as mensagens, pode contar a ocorrência de certas palavras, bem com a chance de elas ocorrerem. Fazendo isso para os e-mails classificados como spam temos…

head(spam.df[with(spam.df, order(-occurrence)),]) term frequency density occurrence 2122 html 377 0.005665595 0.338 538 body 324 0.004869105 0.298 4313 table 1182 0.017763217 0.284 1435 email 661 0.009933576 0.262 1736 font 867 0.013029365 0.262 1942 head 254 0.003817138 0.246

Olha ai, em spam, quase todos os termos tem algo a ver com html. Ja para as mensagens que não são spam…

head(easyham.df[with(easyham.df, order(-occurrence)),]) term frequency density occurrence 3553 yahoo 185 0.008712853 0.180 966 dont 141 0.006640607 0.090 2343 people 183 0.008618660 0.086 1871 linux 159 0.007488344 0.084 1876 list 103 0.004850940 0.078 3240 time 91 0.004285782 0.064

A maior ocorrência, a palavra yahoo, não passa de 18% das mensagens e as frequências de ocorrências são muito baixas. Ou seja, ai em cima a gente vê a informação que permite a classificação de mensagens.

Agora que sabemos como são mensagens de spam e mensagens que não são spam, podemos implementar o classificador. Isso é mais simples que parece, se a frequência de ver html em uma mensagem é 0.30 e a chance de ver a palavra table na mensagem é 0.10, então a chance de ver os dois é 0.30 × 0.10 = 0.03. Agora um problema é quando a gente ve uma palavra que a gente nunca viu quando treinamos o classificador, porque se a gente simplesmente multiplicar por 0, 0% de chance, o resultado final do produtorio vai ser zero, ai vai ser ruim, então aqui, quando uma palavra é nova, é atribuído uma probabilidade bem baixinha do tipo 0.00001% de ocorrência dela tanto para spam como não spam, mas existem maneira mais elaboradas de lidar com esse problema.
Outra coisa é que usamos um prior, para esse partimos do principio que existe 50% de chance de cada mensagem ser spam ou não spam, meio a meio, mas isso não precisa ser necessariamente assim.

Agora se a gente, pegar um novo set de mensagens, a gente faz o processo inverso, se antes nos usamos mensagens de spam e não spam para delimitar as fronteiras do que é ou não um spam, agora é so usar essas fronteiras definidas para dar nome a uma mensagem, se ela é ou não um spam, ou mais precisamente, dar uma chance dela ser ou não um spam, a partir daqui temos que dar um limite de o que será um spam, podemos cometer, como na estatística que vemos comumente, tanto erros do tipo 1 como erros do tipo 2, mas aqui estamos comparando se a chance de ser spam é maior que a chance de ser não spam.

Nos três tipos de dados, os spam, o easyham que são mensagens fáceis de identificar como não spam, e o hardham, que são mais difíceis de identificar, por conter mais termos como html e table por exemplo e ainda sim não serem spam.

03_final_classification

Veja que a linha no meio do gráfico representa a linha de 1 pra 1, quando a probabilidade de ser spam ou não spam são iguais, então pra cima temos uma chance maior de uma mensagem ser spam e para baixo da linha a chance de não ser spam é maior que a chance de ser spam.

O classificador funciona relativamente bem, salvo os erros que passam, mas nada é perfeito, claro que num servidor de e-mail sério, não teremos algo assim tão simples, eu espero, e coisas como lista negra para e-mail ja deve dar uma boa limpada também no que é e-mail de spam. Mas visto isso, machine-learning parece ser algo bem legal, e que pode ser muito útil.

Bem, além de ler o livro, outro lugar legal para dar uma olhada é esse curso do coursera aqui sobre Machine Learning, do Andrew Ng.
Eu me inscrevi nele, e espero começar a postar alguma coisa sobre o tópico a partir de agora. Afinal de contas, curiosidade não mata ninguém. E olhando novamente esse post, ele está mais com cara de agenda que sei la o que.

Por enquanto não temos código, mas a tarefa de casa é usar esse classificador com os dados de exemplo que vem no R chamado iris :), tópico para um próximo post com certeza.