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 – O algoritimo Gradient descent

Vamos denovo meter o nariz na area de Machine Learning, a gente tinha começado a falar disso aqui, foi complicado, mas agora vamos ver algo mais leve e legal.

Alias, na faculdade, em geral, acho que em quase todo curso que tem calculo, todo mundo odeia calculo.

math-rage

Mas um ponto ruim, é que a gente tem que aprender um monte de coisa primeiro, e só depois usar.
Calculo sempre tem aqueles exemplo de velocidade, aceleração, etc. Isso eu acho bem chato.

Mas calculo tem sua utilidade para optimização, que é mais legal que velocidade, na minha opinião. Mas ai você deve pergunta, como isso funciona?

A gente ja comentou, por exemplo aqui, que sempre o que a gente ta fazendo, é ajustando um modelo, que tem uma parte determinística, e uma parte estocástica.

Acontece que o gradiant descent é um algoritmo de optimização, que é baseado em derivada e serve para otimizar.

Mas o que é otimizar?

Quando a gente fala em forrageamento ótimo, como aqui e aqui, a gente tem em mente coisas como ganhar o máximo de energia no menor tempo possível, ou ficar o menor tempo possível exposto, procurando comida.

Aqui é a mesma coisa, queremos normalmente achar o máximo, ou o minimo de alguma função. Mas máximo ou minimo do que? Aqui, para a regressão linear, da verossimilhança, ou likelihood.

Existe um custo, que é o tamanho do erro que vamos ter de acordo com os valores de parametros que assumimos. Custo aqui é o tamanho do desvio que temos para um dado set de parâmetros.

Está confuso certo? Vamos tentar ver numa figura. Suponha que temos uma amostra de 30 valores, e para cada amostra temos o valor de um preditor, que queremos usar ele da melhor forma possível.

Primeiro, vamos determinar aqui, qual o nosso modelo?

Bem, para a regressão linear é:

y=\alpha+\beta \cdot X

Veja que essas letras são whatever, podemos usar a mesma letra com subscritos também:

y=\theta_1+\theta_2 \cdot X

Então a gente tem essa formula, e quer achar os melhores valores de \theta_n para um conjunto de dados de resposta, dado nossos preditores.

Vamos simular, como sempre alguns valores, e olhar como podemos escolher dois valores de \theta_n e dar um parecer, se eles são bons ou ruins, qual o custo deles.

Se a gente assumir \theta_1 = 1 e \theta_2 = 1, temos a formula y=1+1 \cdot X e vemos o seguinte em relação os pontos.

01

Agora essa escolha dos valores \theta_1 = 1 e \theta_2 = 1 é ruim? Não sei, mas podemos calcular a distancia de y até o que esperaríamos de acordo com nossa formula.

> (1+1*X[1])-y[1] [1] -7.841167

Esse -7 é o quão ruim nossa previsão é, e podemos repetir isso para todas as amostras.
Mas dependendo do modelo, esses valores vão ser positivos ou negativos, para evitar isso elevamos essa operação de diminuição ao quadrado, e é só somar.

Em termos matemáticos, o que estamos falando que esse custo é:

Temos nossa hipótese h que é:

h=\theta_1+\theta_2 \cdot X^i

e o custo será:

J(\theta_1,\theta_2)={1\over{2\cdot m}}\sum\limits_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2

Porque esse custo está assim?
Bem o 1\over{2\cdot m} a gente tem que usar, para ponderar o número de amostras, senão, como temos uma somatoria, esse número so vai crescer, crescer, crescer. Como a gente ta dividindo por algo relacionado ao número de amostras, ele ta ficando como uma média.
Aqui, a gente ta dizendo, h_\theta(x^{(i)}), se usarmos essa hipótese h=\theta_1+\theta_2 \cdot X^i, teremos para um valor de X, tal previsão para y.
E no final de tudo, estamos elevando tudo ao quadrado, porque se não for assim, quando a linha estiver passando acima do ponto, da um valor positivo, quando ta abaixo um valor negativo, como tudo elevado ao quadrado é positivo, a gente consegue ver um número relativo a distancia como no gráfico.

Bem se está difícil de entender, podemos escrever isso como uma função no R.

computeCost<- function(X, y, theta) {
    m <- length(y)
    J <- 0
    J = (1/(2*m)) * sum((X%*%theta-y)^2)
    return(J)
}

A gente entra com os valores de X, y e os valores de theta que queremos e calculamos o custo. Que intuitivamente, podemos pensar como algo proporcional a somar todos os pauzinhos vermelhos no gráfico que fizemos e dividir essa soma pelo número de amostras.

Aqui a função deve estar meio confusa, mas porque estamos usando álgebra de matrizes. No R %*% multiplica matrizes. De uma olhada aqui para entender melhor porque isso da certo, mas da, e se eu explicar aqui vou enrolar demais.

Mas o ponto que quero chegar, podemos calcular esse custo, se os parâmetros estão bons ou não, mas isso so faze sentido quando comparamos uma escolha de parâmetros com outra. Aqui se a gente escolher \theta_1 = 5 e \theta_2 = 2, vemos o seguinte:

02

Essa é uma escolha muito melhor. A gente ja tentou uma busca em grid aqui, que consiste em testar sistematicamente todos os valores possiveis. Mas isso é ruim, porque a gente acaba por computar a verossimilhança para um monte de conjunto de parâmetros, valores de \theta que a gente sabe mais ou menos que vão ser ruim.

Ai que entre o gradient descent, e se a toda escolha de parâmetro, você souber pra onde ir, souber como melhorar, se vai melhorar ou não mudar?

Até agora a gente tem uma hipótese.

h=\theta_1+\theta_2 \cdot X

E a gente tem como calcular se estamos escolhendo bons valores ou não para \theta

J(\theta_1,\theta_2)={1\over{2\cdot m}}\sum\limits_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2

Quanto melhor os valores de parâmetros, menor o valor de J. E nosso objetivo e conseguir deixar J o menor valor possível.

Lembra-se que com derivada a gente consegue traçar aquelas restas tangentes, de uma olhada aqui.
Aqui esta o pulo do gato, se ao passar de um ponto para o próximo, a derivada está ficando menos negativa, menos inclinada, a gente ta indo pra baixo, se ela esta ficando positiva, a gente ta subindo, aumentado o J, indo na direção errada.

Então usando derivadas, a gente sabe se deve aumentar ou diminuir os valores de \theta para conseguir um valor menor de J.

Como colocar isso em pratica? Precisamos calcular a derivada do custo para cada parâmetro.

Mas o algoritimo vai funcionar assim:

Passo 1:
Começamos com algum valor arbitrário para \theta_1 e \theta_2.

Passo 2:
Calculamos a derivada para cada parâmetro e atual e fazemos um update do parâmetro atual para um novo valor até não conseguir mais minimizar J(\theta_1,\theta_2).

{\theta_j }:=\theta_j - \alpha \cdot \frac{\partial}{\partial \theta_j} J(\theta_1,\theta_2)

Veja que a gente multiplica por esse \alpha que é o nosso “learning rate”, quanto maior, maior o passo, ou mais brusca a mudança no \theta, e dele ta multiplicando a nossa derivada naquele ponto.

Mas antes de disso a gente vai normalizar os dados para ficar mais bonitos os gráficos depois.
Normalizar os dados significa que vamos simplesmente centralizar eles no zero. A forma dos dados não muda. Veja só.

03

Mudamos os valores no eixos, ele tem agora o centro de X e Y no zero, mas continua a mesma forma geral. Isso implica que o intercepto é zero agora.

Mas tudo bem, o nosso algoritmo vai ficar assim no R.

gradientDescent<- function(X, y, theta, alpha, num_iters) {
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=3)
    for(iter in 1:num_iters) {
        h <- ((X %*% theta)-y)
        temp1 <- theta[1,] - alpha * (1/m) * sum( h * X[,1] )
        temp2 <- theta[2,] - alpha * (1/m) * sum( h * X[,2] )
        theta[1,] = temp1
        theta[2,] = temp2
        J_history[iter,1] = temp1
        J_history[iter,2] = temp2
        J_history[iter,3] = computeCost(X, y, theta);
    }
    return(J_history)
}

Então recebemos os dados de X e y, os parâmetros \theta (isso se chama theta) iniciais, e um número de iterações que queremos realizar, quantos updates queremos fazer.

Dai a gente vai ver quantos dados temos e salvar em m, preparar uma matriz J_history para guardar os valores de \theta conforme vamos mudando, para olhar de pois, bem como os valores de J. Calculamos a derivada e fazemos simultaneamente o update para os dois valores de \theta e guardamos eles e fazemos isso até rodar o número num_iters de iterações que pedimos.

Vamos começar com uns valores de \theta bem ruins.

> theta<-matrix(c(-5,-5),ncol=1,nrow=2) > theta [,1] [1,] -5 [2,] -5 > computeCost(X, y, theta) [1] 29.82811

A gente aplica o algoritimo gradient descent

> iterations = 1500 > alpha = 0.01 > resultado<-gradientDescent(X, y, theta, alpha, iterations) > resultado [,1] [,2] [,3] [1,] -4.950000e+00 -4.942143785 29.24624141 [2,] -4.900500e+00 -4.884846846 28.67573437 [3,] -4.851495e+00 -4.828103778 28.11636476 . . . [1498,] -1.447077e-06 0.985122862 0.01427153 [1499,] -1.432606e-06 0.985122890 0.01427153 [1500,] -1.418280e-06 0.985122918 0.01427153

Aqui, a coluna 1 é o intercepto, a coluna 2 é a inclinação e a coluna 3 é o custo.
A gente começou em 29, e veja que no final, estava em 0.0143 o custo e não dava mais para sair disso.

Se a gente fizer uma regressão linear, que é outra forma de caçar os melhores parâmetros, como visto aqui.

lm(y~X[,2]) Call: lm(formula = y ~ X[, 2]) Coefficients: (Intercept) X[, 2] -2.238e-17 9.851e-01

Vemos praticamente a mesma inclinação e intercepto. Não confunda os número por causa da notação cientifica. O intercepto é um 0,000…2238, um numero muito próximo de zero, e a inclinação é
0.985, bem próximo do 0.985122918 que achamos com o gradient descent.

Se a gente for olhar como os valores do custo foram mudando, veremos o seguinte:

04

Ele foi diminuindo, mais rápido no começo, quando a inclinação da derivada era maior, depois mais devagar quando a inclinação da derivada era menor, até não mudar, quando a inclinação chegou próximo de zero.

Se comparar o modelo de regressão e os resultados do gradient descent, vemos que chegamos a mesma definição de parâmetros.

05

Sendo que a gente foi mudando os parâmetros de forma certeira, sempre minimizando o custo.

06

A gente pode fazer uma analogia dessas mudança como a escalada de uma montanha, o topo da montanha (menor custo) é onde a gente quer chegar, e dado que saímos de um ponto, vamos subindo sempre pra cima, nunca para baixo, até chegar ao cume, mas a gente começa indo rápido e diminuindo o passo ao longo da subida para não passar o cume.

07

A gente pode pegar e olhar, nessa montanha, como foram os nossos 200 primeiros passos, para ter uma ideia.

08

Muito legal certo? E um uso bem legal de derivada, conseguir andar certeiro até o maximum likelihood.
Tudo bem que para regressão linear, ninguém vai usar isso, mas a gente viu que isso funciona, e no caso de analises como o MDS (Multidimensional scaling), a gente ta fazendo algo nesse sentido, a gente tem uma matrix de distancia e esta tentando posicionar os pontos num plano, ou reta, ou sei la, mas estamos tentado posicionar e medindo se está ficando bom de acordo com uma função de custo, uma loss function.

Veja que é fácil modificar para outras funções(desde que sejam derivaveis) o gradient descent.
Peguemos o caso de uma regressão multiplica. Basta analisar para n valores de \theta.

gradientDescentMulti<-function(X, y, theta, alpha, num_iters) {
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=1+nrow(theta))
    for(iter in 1:num_iters) {
        for(j in 1:length(theta)) {
            theta[j,]<-theta[j,] - alpha * (1/m) * sum( ((X %*% theta)-y) * X[,j] )
            J_history[iter,j] = theta[j,]
        }
        J_history[iter,nrow(theta)+1] = computeCost(X, y, theta);
    }
    return(J_history)
}

Veja que eu adicionei um loop ali dentro para cada valor de \theta, ao invés das duas linhas no caso anterior, mas a ideia continua a mesma. So com mais \theta.

Aqui a nossa hipótese se torna:

h=\theta_1+\theta_2 \cdot X1^i+\theta_3 \cdot X2^i

Um modelo aditivo, com uma soma a mais, iniciamos os \theta e rodamos nosso algoritmo gradient descent

> theta<-matrix(c(0,0,0),ncol=1,nrow=3) > theta [,1] [1,] 0 [2,] 0 [3,] 0 > computeCost(X, y, theta) [1] 0.4833333 > gradientDescentMulti(X, y, theta, alpha, 1500) [,1] [,2] [,3] [,4] [1,] -1.846902e-18 0.005242762 0.005947532 0.47707775 [2,] -3.745557e-18 0.010452588 0.011853113 0.47090590 [3,] -5.621950e-18 0.015629671 0.017717052 0.46481668 . . . [1498,] -1.664308e-16 0.808733080 0.863209808 0.01529028 [1499,] -1.664447e-16 0.808733320 0.863210047 0.01529028 [1500,] -1.665603e-16 0.808733558 0.863210284 0.01529028

E temos nas três primeiras colunas os parâmetros e a última sendo o valor do custo J.

> lm(y~X[,2]+X[,3]) Call: lm(formula = y ~ X[, 2] + X[, 3]) Coefficients: (Intercept) X[, 2] X[, 3] -1.730e-16 8.088e-01 8.632e-01

Comparamos com o resultado da regressão linear e vemos que ele está funcionando perfeito novamente.
Podemos ainda avaliar a interação, sem mexer em nada na função.

A hipótese fica assim:

h=\theta_1+\theta_2 \cdot X1^i+\theta_3 \cdot X2^i+ \theta_4 \cdot X1 \cdot X2

Rodando a função vemos o seguinte.

> theta [,1] [1,] 0 [2,] 0 [3,] 0 [4,] 0 > computeCost(X, y, theta) [1] 0.4833333 > gradientDescentMulti(X, y, theta, alpha, 1500) [,1] [,2] [,3] [,4] [,5] [1,] 1.170360e-18 0.005229905 0.007172164 0.004086374 0.4738325917 [2,] 3.297842e-06 0.010411918 0.014276559 0.008114284 0.4645249067 [3,] 9.813365e-06 0.015546491 0.021313839 0.012084525 0.4554062184 . . . [1498,] 2.092109e-02 0.585063765 0.773159318 0.259232115 0.0001656415 [1499,] 2.092109e-02 0.585063777 0.773159330 0.259232111 0.0001656415 [1500,] 2.092109e-02 0.585063789 0.773159342 0.259232106 0.0001656415

E denovo conseguimos os mesmos valores de parâmetros da regressão.

> lm(y~X[,2]*X[,3]) Call: lm(formula = y ~ X[, 2] * X[, 3]) Coefficients: (Intercept) X[, 2] X[, 3] X[, 2]:X[, 3] 0.02092 0.58507 0.77316 0.25923

Bem, é isso. As vezes ver uma utilidade diferente para derivadas pode tornar menos massante estudar calculo. Eu achei bem legal quando vi esse algoritimo na aula de aprendizado de máquina nesta matéria do coursera aqui. Ela é muito legal, mas infelizmente tudo nela está na linguagem do matlab e octave, mas não é difícil passar as coisas para o R, se eu consegui, qualquer um consegue.

Bem temos o script aqui em baixo, esse foi o tópico das duas primeiras semanas, se eu conseguir, vou tentar passar a parte de modelos de classificação para o R também, e quem sabe, com sorte, eu dou conta de passar a parte de rede neural que é super legal, mas façam o curso la, que é mil vezes melhor explicado que aqui e está no começo de uma turma.

Ficamos por aqui, os scripts sempre no repositório recologia e é isso, até o próximo post.

computeCost<- function(X, y, theta) {
    m <- length(y)
    J <- 0
    J = (1/(2*m)) * sum((X%*%theta-y)^2)
    return(J)
}
 
gradientDescent<- function(X, y, theta, alpha, num_iters) {
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=3)
 
    for(iter in 1:num_iters) {
 
        h <- ((X %*% theta)-y)
 
        temp1 <- theta[1,] - alpha * (1/m) * sum( h * X[,1] )
        temp2 <- theta[2,] - alpha * (1/m) * sum( h * X[,2] )
 
        theta[1,] = temp1
        theta[2,] = temp2
 
        J_history[iter,1] = temp1
        J_history[iter,2] = temp2
 
        J_history[iter,3] = computeCost(X, y, theta);
 
    }
 
    return(J_history)
}
 
gradientDescentMulti<-function(X, y, theta, alpha, num_iters) {
 
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=1+nrow(theta))
 
    for(iter in 1:num_iters) {
 
        for(j in 1:length(theta)) {
            theta[j,]<-theta[j,] - alpha * (1/m) * sum( ((X %*% theta)-y) * X[,j] )
            J_history[iter,j] = theta[j,]
        }
 
        J_history[iter,nrow(theta)+1] = computeCost(X, y, theta);
    }
 
    return(J_history)
}
 
#gerando  dados
X<-runif(30,1,10)
y<-rnorm(30,5+2*X)
 
#figura 1
plot(y~X,main="Modelo 1+1*X",frame=F,xlim=c(0,10),ylim=c(0,30))
abline(a=1,b=1)
for(i in 1:length(X)) {
    points(c(X[i],X[i]),c(y[i],1+1*X[i]),type="l",lty=1,col="red",lwd=2)
}
 
#figura 2
plot(y~X,main="Modelo 5+2*X",frame=F,xlim=c(0,10),ylim=c(0,30))
abline(a=5,b=2)
for(i in 1:length(X)) {
    points(c(X[i],X[i]),c(y[i],5+2*X[i]),type="l",lty=1,col="blue",lwd=2)
}
 
#figura 3
par(mfrow=c(2,1))
plot(y~X,main="Dados sem normalização",frame=F)
abline(v=mean(X),h=mean(y),lty=2)
plot(scale(y)~scale(X),main="Dados normalizados",frame=F)
abline(v=mean(scale(X)),h=mean(scale(y)),lty=2)
 
#transformando em matrix
X<-matrix(c(rep(1,length(X)),scale(X)),ncol=2)
y<-matrix(scale(y),ncol=1)
 
#gerando theta inicial
theta<-matrix(c(-5,-5),ncol=1,nrow=2)
theta
 
#custo inicial
computeCost(X, y, theta)
 
#outros parametros
iterations = 1500
alpha = 0.01
 
 
#gradient descent
 
resultado<-gradientDescent(X, y, theta, alpha, iterations)
 
resultado
lm(y~X[,2])
#figura 4
plot(resultado[1:200,3],type="b")
 
#figura 5
plot(y~X[,2],frame=F)
abline(lm(y~X[,2]),col="red",lty=2)
abline(a=resultado[1500,1],b=resultado[1500,2],col="blue",lty=3)
legend("topleft",lty=c(2,3),col=c("red","blue"),legend=c("Regressão Linear","Gradient Descent"))
 
#figura 6
plot(resultado[,1],resultado[,2],type="b")
 
a<-seq(-2,2,by=0.1)
b<-seq(-2,3,by=0.1)
J<-matrix(0,nrow=length(a),ncol=length(b))
 
for(i in 1:length(a)) {
    for(j in 1:length(b)) {
        J[i,j]<--1*computeCost(X, y, matrix(c(a[i],b[j]),ncol=1,nrow=2))
 
    }
}
 
#figura 7
persp(a,b,J,theta=45,phi = 30,xlab="Intercepto",ylab="Inclinação",zlab="Menor custo")
 
#figura 8
contour(a,b,J)
points(resultado[1:300,1],resultado[1:300,2],type="b",cex=0.3)
 
##Regressão Multipla
X1<-runif(30,1,10)
X2<-runif(30,1,10)
y<-rnorm(30,5+2*X1+2*X2)
 
X<-matrix(c(rep(1,length(X1)),scale(X1),scale(X2)),ncol=3)
y<-matrix(scale(y),ncol=1)
 
theta<-matrix(c(0,0,0),ncol=1,nrow=3)
theta
 
computeCost(X, y, theta)
 
gradientDescentMulti(X, y, theta, alpha, 1500)
lm(y~X[,2]+X[,3])
 
X1<-runif(30,1,10)
X2<-runif(30,1,10)
y<-rnorm(30,5+2*X1+2*X2+3*X1*X2)
 
X<-matrix(c(rep(1,length(X1)),scale(X1),scale(X2),scale(X1)*scale(X2)),ncol=4)
y<-matrix(scale(y),ncol=1)
 
theta<-matrix(c(0,0,0,0),ncol=1,nrow=4)
theta
 
computeCost(X, y, theta)
 
gradientDescentMulti(X, y, theta, alpha, 1500)
lm(y~X[,2]*X[,3])

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.