Se você não sabe ainda o que é um modelo de resposta binária, dê uma olhada aqui.
Esses modelos já estão implementados em Python, veja por exemplo, aqui. Também já estão implementados em R ou Matlab.
A resposta que darei aqui implementa explicitamente o classificador em Python focando o entendimento e aprendizado. Outras pessoas poderiam dar uma resposta explicitando rapidez de impletamentação utilizando as ferramentas já prontas do Python, R ou Matlab.
Utilizaremos a notação para o modelo de resposta binária considerada nessa outra resposta.
Para testar nossa implementação, vamos usar dados gerados utilizando simulações Monte-Carlo. A forma mais interessante de gerar amostras de simulações monte carlo para modelos de resposta binária é:
1) Crie \(n\) valores de um vetor bidimensional de regressores com uma distribuição pre-especificada.
2) Para cada um desses \(n\) valores, calcule \(F(z_i)\), onde \(z_i=\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}\), usando os \(\beta\)s reais (lembre-se que estamos fazendo o papel da natureza quando usamos a simulação de Monte-Carlo) e a distribuição acumulada \(F\)
desejada (nesse caso usamos normal, mas poderiamos ter usado a distribuição logística).
3) Sorteie \(n\) valores de de uma variável aleatória uniforme \(w_i,\;\;i=1,\cdots,n\) no intervalo \([0,1]\);
4) Gere os \(n\) valores de \(y_i\) fazendo o seguinte teste: Se \(w_i\in [0,F(x_i'\beta)]\) então \(y_i=1\). Em caso contrário, \(y_i=0,\;\;i=1,\cdots,n\).
Entenda os passos (3) e (4) com cuidado. Na vida real, muitas vezes temos características que nos sugerem que vamos tomar um determinado tipo de atitude e tomamos uma atitude diferente daquela prevista (erros de modelagem). Então esses passos nos dizem que quanto maior a probabilidade medida por \(F\), maior a chance de eu me comportar como esperado.
Apenas para comparação vamos também gerar valores de \(y_i\) que refletem perfeitamente as características individuais. Substitua os passos (3) e (4) pelo passo (3\(^\prime\)):
3\(^\prime\)) Gerar os $n$ valores de \(y_i\) fazendo o seguinte teste: Se \(F(x_i'\beta)>1/2\) então \(y_i=1\). Em caso contrário, \(y_i=0,\;\;i=1,\cdots,n\).
A figura abaixo apresenta essas duas simulações acima. Se usarmos o passo (3\(^\prime\)) geramos dados que são perfeitamente separáveis e em caso contrário geramos dados que refletem mais situações reais e não são perfeitamente separáveis.
A reta usada para classificação é exatamente
\[z=\beta_0 + \beta_1x_1 +\beta_2 x_2=0\]
Note que a reta contínua utiliza os parâmetros reais e a reta tracejada usa os parâmetros estimados.
Note também que podemos provar geometricamente que a distância de um ponto \((x_1,x_2)\) a reta usada para classificação nos dá uma medida de quanto ela atende bem as características daquela região do classificador.
O código utilizado para essa implementação em Python está aqui (se você ainda não sabe Python e gostaria de aprender, siga por aqui):
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import fmin
import numpy as np
import math as math
np.random.seed(50)
plt.figure(num=None, figsize=(8, 4), dpi=80, facecolor='w', edgecolor='k')
def geraClassificador(beta0,beta1,beta2,x1,parameter):
minx1=np.min(x1)
maxx1=np.max(x1)
minx2=(beta0+beta1*minx1)/(-beta2)
maxx2=(beta0+beta1*maxx1)/(-beta2)
if(parameter):
plt.plot([minx1,maxx1],[minx2,maxx2],'k-')
else:
plt.plot([minx1,maxx1],[minx2,maxx2],'y--')
def geracaoDados(n,beta0,beta1,beta2):
x1=np.random.normal(1,0.64,n) # Primeiro regressor
x2=np.random.normal(1,0.36,n) # Segundo regressor
z=beta0*np.ones(n)+beta1*x1+beta2*x2 # Variavel latente
w=np.random.uniform(0,1,n)
y=np.empty([n])
yP=np.empty([n])
for i in range(n):
if(stats.norm.cdf(z[i])>0.5):
yP[i]=1
else:
yP[i]=0
if(w[i]< stats.norm.cdf(z[i])):
y[i]=1
else:
y[i]=0
return yP,y,x1,x2
def probitEst(betas,*data):
(beta0,beta1,beta2)=betas
y,x1,x2=data
n=np.size(y)
mLike=0;
for i in range(n):
z=beta0+beta1*x1[i]+beta2*x2[i]
F=stats.norm(0, 1).cdf(z)
if((F != 0) and (F != 1)):
mLike=mLike+ y[i]*math.log(F) + (1-y[i])*math.log(1-F);
return -mLike
def estimacaoRegBin(y,x1,x2):
data=(y,x1,x2)
beta00=0.5
beta10=0.5
beta20=0.5
betas = fmin(probitEst, [beta00,beta10,beta20], data,xtol=1e-4)
return betas
beta, alpha, r_value, p_value, std_err = stats.linregress(x,y)
return alpha,beta
if __name__ == '__main__':
n=100 # numero de observacoes
beta0=2.5
beta1=-1.5
beta2=-1
[yP,y,x1,x2]=geracaoDados(n,beta0,beta1,beta2)
axes=plt.subplot(121)
plt.hold(True)
for i in range(n):
if(yP[i]==1):
plt.plot(x1[i],x2[i],u'ro',markersize=5, markeredgewidth=0)
else:
plt.plot(x1[i],x2[i],u'bo',markersize=5, markeredgewidth=0)
axes.set_xlabel('$x_1$')
axes.set_ylabel('$x_2$')
axes.set_title('Perfeitamente Separavel')
geraClassificador(beta0,beta1,beta2,x1,True)
[hatBeta0,hatBeta1,hatBeta2]=estimacaoRegBin(yP,x1,x2)
print hatBeta0,hatBeta1,hatBeta2
geraClassificador(hatBeta0,hatBeta1,hatBeta2,x1,False)
axes=plt.subplot(122)
plt.hold(True)
for i in range(n):
if(y[i]==1):
plt.plot(x1[i],x2[i],u'ro',markersize=5, markeredgewidth=0)
else:
plt.plot(x1[i],x2[i],u'bo',markersize=5, markeredgewidth=0)
axes.set_xlabel('$x_1$')
axes.set_ylabel('$x_2$')
axes.set_title('Nao Perfeitamente Separavel')
geraClassificador(beta0,beta1,beta2,x1,True)
[hatBeta0,hatBeta1,hatBeta2]=estimacaoRegBin(y,x1,x2)
print hatBeta0,hatBeta1,hatBeta2
geraClassificador(hatBeta0,hatBeta1,hatBeta2,x1,False)
Talvez alguém possa contribuir com o mesmo código em R ou Matlab em uma outra resposta. Mais detalhes sobre modelos de resposta binária e essa implementação você pode encontrar nessa referência.