Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como implementar o modelo conhecido como Linear Discriminant Analysis (LDA) em python?

0 votos
248 visitas
perguntada Jun 21, 2016 em Estatística por CaioP (21 pontos)  

Considere usar dados simulados e os dados do exercício 8.8 de IZENMAN, Alan Julian. Modern multivariate statistical techniques. New York: Springer, 2008.

Compartilhe

1 Resposta

0 votos
respondida Jun 21, 2016 por CaioP (21 pontos)  

LDA é uma forma de reduzir a dimensionalidade do dados e aplicar sobre o sub-espaço gerado, um critério de separação linear.

\( y=xw \)

Nesse caso, a questão é como escolher a matriz \(w\) de modo a minimizar a perda de informação. Ou seja, o sub-espaço \(y\) contém aquilo que é mais essencial para dividir os dados nas categorias do problema.

O desenvolvimento dessa resposta segue 1, que descreve como calcular \(w\) através do seguintes cinco passos:

1) compute as médias de cada variável dos dados para cada categoria;
2) compute as matrizes de dispersão (scatter matrix) de cada categoria e entre categorias;
3) compute os auto-valores e auto-vetores das matrizes de dispersão;
4) construa \(w\) com os \(k\) auto-vetores associados aos \(k\) maiores auto-valores; e
5) use a matriz resultante para gerar o sub-espaço \(y\).

O código abaixo em python é uma implementação desses passos:

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

np.random.seed(50)

def step1_mean_each_class(X,num_cat):
    mean_vector = []
    for i in range(0,num_cat):
        mean_vector.append(np.mean(X[X[:,len(X[1,:])-1]==i,:],axis=0))
    return mean_vector

def step2a_scatter_within_class(X,num_cat):
    scat = np.zeros((len(X[1,:])-1,len(X[1,:])-1))
    m = step1_mean_each_class(X,num_cat)
    for j in range(0,num_cat):
        aux_X = X[X[:,len(X[1,:])-1]==j,:len(X[1,:])-1]
        aux_m = m[j]
        for l in range(0,len(aux_X[:,1])):
            aux_A = np.array(aux_X[l,:]-aux_m[:len(X[1,:])-1])
            aux_a = aux_A.reshape(len(aux_A),1)
            aux_b = aux_A.reshape(1,len(aux_A))
            aux_c = np.dot(aux_a,aux_b)
            scat = scat + aux_c
    return scat

def step2b_scatter_between_class(X,num_cat):
    scat = np.zeros((len(X[1,:])-1,len(X[1,:])-1))
    m = step1_mean_each_class(X,num_cat)
    m_aux = np.mean(X,axis=0)
    m_gen = m_aux[:len(m_aux)-1]
    for i in range(0,num_cat):
         aux_X = X[X[:,len(X[1,:])-1]==i,:len(X[1,:])-1]
         aux_m = m[i]
         aux_A = np.array(aux_m[:len(X[1,:])-1]-m_gen)
         aux_a = aux_A.reshape(len(aux_A),1)
         aux_b = aux_A.reshape(1,len(aux_A))
         aux_c = np.dot(aux_a,aux_b)
         scat = scat + len(aux_X[:,0])*aux_c
    return scat

def step3_eigen(X,num_cat):
    S_W = step2a_scatter_within_class(X,num_cat)
    S_B = step2b_scatter_between_class(X,num_cat)
    eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
    return eig_vals, eig_vecs

def step4_get_big_eigen(X,num_cat,num_reg):
    eig_vals, eig_vecs = step3_eigen(X,num_cat)
    eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]
    eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)
    for i in eig_pairs:
        print(i[0])
    W = np.hstack((eig_pairs[0][1].reshape(num_reg,1), eig_pairs[1][2].reshape(num_reg,1)))
    return W

def step5_new_sub(X,num_cat,num_reg):
    W = step4_get_big_eigen(X,num_cat,num_reg)
    aux_X = X[:,:len(X[0,:])-1]
    X_lda = aux_X.dot(W)
    return X_lda

Com esse código, podemos então usar uma simulação de monte carlo para gerar a base de dados e visualizar o gráfico. Ao usarmos um processo gerador de dados como em 2, poderemos checar a eficácia da implementação (deve separar perfeitamente as classes). Isso é feito com o código abaixo:

def data_gen(n,beta0,beta1,beta2,beta3,beta4,var):
    x1=np.random.normal(0,var,n) 
    x2=np.random.normal(0,var,n) 
    x3=np.random.normal(0,var,n) 
    x4=np.random.normal(0,var,n) 
    z=beta0*np.ones(n)+beta1*x1+beta2*x2+beta3*x3+beta4*x4 
    y=np.empty([n])
    for i in range(n):
        if(z[i]>beta0):
            y[i]=1
        else:
            y[i]=0
    data = np.column_stack((x1,x2,x3,x4,y))
    return data

data=data_gen(100,0,10,2,3,1,0.16)

E o gráfico é feito da seguinte forma, como em 1:

def plot_step_lda(X,num_cat,num_reg):
    label_dict = {0: 'y=0', 1: 'y=1'}
    ax = plt.subplot(111)
    X_lda = step5_new_sub(X,num_cat,num_reg)
    for label,marker,color in zip(
        range(0,2),('s', 'o'),('blue', 'red')):
        plt.scatter(x=X_lda[:,0].real[X[:,num_reg] == label],
                y=X_lda[:,1].real[X[:,num_reg] == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )
    plt.xlabel('LD1')
    plt.ylabel('LD2')
    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title('LDA')

    # hide axis ticks
    plt.tick_params(axis="both", which="both", bottom="off", top="off",  
            labelbottom="on", left="off", right="off", labelleft="on")

    # remove axis spines
    ax.spines["top"].set_visible(False)  
    ax.spines["right"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["left"].set_visible(False)    

    plt.grid()
    plt.tight_layout
    plt.show()

plot_step_lda(data,2,4)

A imagem será apresentada aqui.

Se usarmos os dados reais sugeridos no enunciado, a figura não é tão clara, como esperado.

df = pd.read_excel("path")
data = df.as_matrix()
step5_new_sub(data,2,85)

A imagem será apresentada aqui.

Note que foi feita uma opção por usar, nos códigos, apenas dois eixos na redução de dimensionalidade, mais ou menos eixos poderiam ser usados. Tal opção se deve a facilidade de obter gráficos em duas dimensões. Inclusive, como podemos ver, um dos eixos sequer é muito informativo. Assim, separação em categorias poderia ser facilmente feita se estipulássemos um critério numa redução para o subespaço de uma reta (usando o auto-vetor mais informativo).

...