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()

O problema dos pesos de Bachet de Méziriac

0007r

Um mercador tinha um peso de 40 quilos que quebrou um quatro pedaços. Quando os pedaços foram pesados, ele viu que cada pedaço era um número inteiro em quilos, e que os quatro pedaços podiam ser usados para medir qualquer inteiro de 1 a 40 quilos. Quais foram os pesos em que esses pedaços foram quebrados?

E ai? Bem, primeiro a gente pode pensar na balança (que tem dois lados, daquelas antigas) como o lado do peso e o lado da medida.

No lado do peso, nos colocamos somente peças do peso, já no lado da media, nos podemos colocar o que vamos medir mais pesos adicionais. Se nos queremos pesar com o menor número de pesos possíveis, temos que colocar pesos dos dois lados na maioria dos casos. Por exemplo, para pesar 1 quilo tendo pesinhos de 2 e 3 quilos, nos colocamos 3 quilos no lado do peso e 2 quilos no lado da medida, dai se colocarmos algo no lado da medida e a balança equilibrar, bingo, o que colocamos pesa 1 quilo. Com esses dois pesos a gente pesaria medida de 1, 2, 3, e 5 quilos, mas não 4 quilos, porque não daria para combinar 2 e 3 para conseguir averiguar se algo pesa 4 quilos.

Nos vamos considerar apenas cargas e pesinhos integrais, totalizando sempre totais inteiros.

Se a gente tem uma série de pesinhos A, B, C, …, que nos permite, quando distribuídos apropriadamente sobre as escalas, pesar todas as cargas de 1 até n quilos, e se P é uma nova medida tal que seu peso p excede o total de peso n que conseguimos formar com os pesinhos anteriores por 1.

p-n=n+1

Ou isolando p

p=2n+1

é possível pesar todos as cargas integrais de 1 até p+n = 3n+1 pela adição do peso P para o conjunto de pesos A, B, C, ….

Muita calma nessa hora, primeiro pense que estamos falando que queremos pesar algo mais pesado que todos os pesinhos que temos junto, logo, precisamos de mais pesinhos, e se pegarmos um peso P que pese 2n+1, nosso conjunto de possibilidades que era 1 até n, passa a ser de 1 até 3n+1, combinando esse novo peso com os anteriores.

Primeiro, se a gente tiver somente um pesinho A, a menor quantidade que a gente quer pesar é 1, então o pesinho A tem que pesar 1, agora somente com 1 peso, se a gente quer pesar uma carga de 2, esse peso A ja não resolve, precisamos de um segundo peso B

Para conseguir pesar todos os inteiros com 2 pesos, A e B, A tem que pesar 1 quilo e B 3 ((2 \cdot 1) + 1) quilos. Com esses dois pesos a gente pode pesar 1, 2, 3 e 4 quilos.

Vamos pensar se está dando certo,

Para 1 quilo, eu ponho 1 num lado e o outro ponho a carga.
Para 2 quilos, eu ponho 3 de um lado e do do outro a carga mais o pesinho de 1 quilo
Para 3 quilos eu ponho 3 de um lado e a carga do outro
Para 4 quilos eu ponho o 3 e 1 de um lado e a carga do outro lado.

Opa funciona. Agora isso so funciona até 4 quilos, a gente pode pesar de 1 … n, sendo o n = 4 quilos, para pesar 5 quilos, a gente já precisa de mais um pesinho C.

Se a gente escolher um terceiro peso C de forma a ele pesar

c=2 \cdot 4+1 = 9

então é possível usar esses três pesos para medir todas as cargas inteiras de 1 até c+4=9+4=13.

Pensando que a carga vale x e temos então pesinhos de 1,3,9.

Para 01 usamos 1=x
Para 02 usamos 3=x+1
Para 03 usamos 3=x
Para 04 usamos 3+1=x
Para 05 usamos 9=x+3+1
Para 06 usamos 9=x+3
Para 07 usamos 9+1=x+3
Para 08 usamos 9=x+1
Para 09 usamos 9=x
Para 10 usamos 9=x+1
Para 11 usamos 9+3=x+1
Para 12 usamos 9+3=x
Para 13 usamos 9+3+1=x

Se eu não me confundi para algum número, a gente viu como com esses pesos 1,3 e 9 da para pesar tudo de 1 até 13. Agora como precisamos chegar a 40, o negocio é adicionar mais um peso ainda, escolher um quarto peso D de forma que ele pese

d=2 \cdot 13+1 = 27

Os quatro pesos A, B, C e D vão permitir pesar todas as cargas de 1 a 27+13 = 40 quilos. Não vou fazer um por um aqui, mas ja da para acreditar que vai dar certo.

Logo o peso se partiu em pedaços de 1, 3, 9 e 27 quilos.

Bem legal esse problema né? Ele foi proposto pelo matemático Claude Gaspard Bachet de Méziriac, que resolveu esse problema no seu livro famoso Problèmes plaisants et délectables qui se font par les nombres (esse link é a digitalização de uma publicação original, não entendo nada mas é legal), publicado em 1624.

Referência: Henrich Dorie 100 great problems of elementary mathematics – Their history and solutions

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")

Evolução no HIV

Bem, eu li recentemente um artigo bem legal na PNAS.

Impact of HLA-driven HIV adaptation on virulence in populations of high HIV seroprevalence

Os caras estudaram a evolução do vírus da imunodeficiência humana em Botswana, onde a epidemia de HIV começou e teve mais prevalência.

Acontece que algumas pessoa tem um gene chamado HLA-B*57, ou alelo (acho que o gene é para outras coisas, não sou muito bom em genética), mas esse carinha confere uma certa resistência ao HIV, dificultando a reprodução dele e assim tornando ele menos virulento.

Massa, nisso então o vírus vai matando todo mundo, menos os caras que tem esse HLA-B*57, e se o vírus não der um jeito, ele pode matar todos os hospedeiros possíveis e ficar sem ninguém para se replicar e se extinguir.

Então ele vai mudando para se adaptar a esse alelo do gene. So que nisso ele foi diminuindo a virulência para todas as pessoas. Porque esse é o “jeito” de sobreviver do HIV.

Um ponto que me faz refletir, é como as vezes eu penso que a medicina está evoluindo e criando soluções melhores para curar o HIV. Mas esse é um pensamento simplista demais, o mundo é muito mais dinâmico, assim como os tratamentos mudam, a propria doenca muda, ja que de certa forma o HIV é uma especie (tudo bem que indesejada por nos, e salvo as definições de especies) mudando, então a prevalência do HIV mudando ao longo do tempo pode ter haver com a evolução dos tratamentos e com a evolução do próprio HIV.

Ainda no trabalho, os caras fizeram um modelinho que sugere que como o tratamento sido distribuído mais amplamente, e principalmente para quem esta mais “ferrado”, as linhagens mais virulentas vão se reproduzindo menos, tanto pelo tratamento como menor chance de transmissão (quem esta mais doente tem menos chance de passar o pra frente, relativamente), o que vai selecionando as linhagens menos virulentas para se manter na população, o que acelera a evolução para um HIV menos virulento.

Muito legal, so que acho que não podemos somente confiar nisso, mas ha muito tempo se fala em tornar mais abrangente o tratamento ao HIV, principalmente onde a prevalência é maior, como em alguns países africanos, o duro é o custo disso, já que envolve tanto o tratamento como varias patentes.

Mas achei muito legal esse artigo e imagino que quem frequenta aqui, deve achar também, como uma leitura aleatória para a semana :).

Aqui e aqui tem duas reportagens sobre ele também, com um texto mais leve que o do artigo, sem tantos detalhes técnicos.