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

Como implementar o modelo de regressão conhecido como Principal Components Regression (PCR) em python?

+1 voto
213 visitas
perguntada Jun 17, 2016 em Ciência da Computação por Saulo (426 pontos)  

Como implementar o modelo de regressão conhecido como Principal Components Regression (PCR) em python? Exemplifique o uso do seu modelo considerando a questão 5.10 do livro "Modern multivariate statistical techniques" - Alan Julian Izenman que versa sobre o uso da base de dados "detergent" e PCR, PLSR e RR implementando explicitamente esses modelos de regressão. Interprete os resultados e documente cuidadosamente o seu material e implementação.

Compartilhe

4 Respostas

+1 voto
respondida Jun 17, 2016 por Saulo (426 pontos)  

Introdução

Análise de componentes principais (PCA) é uma técnica que pode ser aplicada em diversas situações, tais como: redução de dimensão, compressão de dados, extração de características e visualização de dados. Ela pode ser definida de duas formas equivalentes: (1) PCA é a projeção ortogonal dos dados em um espaço linear de dimensão menor tal que a variância dos dados projetados é maximizada; (2) PCA é projeção linear que minimiza a média do custo médio projetado (média do quadrado da distância entre os pontos dos dados e suas projeções). Mais detalhes e demonstrações podem ser obtidos em [2, seção 12.1, p. 561-570] e [1, seção 7.2, p. 196-215].

O modelo de regressão conhecido como PCR (Principal Components Regression PCR ou Regressão sobre os componentes principais) utiliza a análise de componentes principais para fazer a regressão. Dessa forma, consegue reduzir a dimensão da regressão excluindo as dimensões que contribuem para causar problema de colinearidade. Vamos explicar a seguir a técnica mostrada em [1, p. 131-133].

Referencial teórico

Considere o modelo \(\mathbf{Y} = \mathbf{X}^{T} \boldsymbol{\beta} + \mathbf{e}\), onde \(\mathbf{Y}_{N \times 1} = (y_1, \ldots, y_N)^T\) é o vetor da variável dependente, \(\mathbf{X}_{N \times K} = (\mathbf{x}_1, \ldots, \mathbf{x}_K)^T\) é matriz de regressores, \(N\) é o total de observações, \(K\) é o número de variáveis independentes, e \(\boldsymbol{\beta} \in \mathbb{R}^{K}\) é o vetor de parâmetros, os coeficientes da regressão, a ser estimado. Suponha que \(\mathbf{X}^{T}\mathbf{X}\) tenha posto \(r\), \(1 \leq r \leq R = \min{\{N,K\}}\), de forma que os \(R-r\) autovalores de \(\mathbf{X}^{T}\mathbf{X}\) são iguais a zero.

Para fazer a regressão nas componentes principais, inicialmente, precisamos encontrar os autovalores e autovetores de \(\mathbf{X}^{T}\mathbf{X}\). Para encontrar os autovalores e os autovetores, vamos usar a decomposição em valores singulares ("Singular Value Decomposition" ou SVD), pois é garantido encontrar uma solução em tempo finito para as técnicas padrões de PCA baseadas em SVD, além dos autovalores já estarem ordenados, conforme observado por [2, p. 593] em outro contexto. Conforme mostrado em [1, p. 138], SVD nos fornece que \(\mathbf{X} = \mathbf{U} \boldsymbol{\Lambda}^{1/2} \mathbf{V}^T\), onde \(\boldsymbol{\Lambda} = diag[\lambda_j]\) e \(\mathbf{X}^T \mathbf{X} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^T\). Os autovalores ordenados de \(\mathbf{X}^T \mathbf{X}\) são dados por \(\{\lambda_j\}\), com \(\lambda_1 \geq \cdots \geq \lambda_r \geq 0\). Os autovetores \(\mathbf{v}_j\) são as colunas de \(\mathbf{V} = (\mathbf{v}_1, \ldots, \mathbf{v}_r)\) associadas aos autovalores de \(\lambda_j\), e correspondem à \(j\)-ésima direção da componente principal. \(\mathbf{X}\mathbf{v}_j\) é a \(j\)-ésima componente principal.

Em seguida, precisamos criar a matriz \(\mathbf{Z}_t = \mathbf{X} \mathbf{V}_t = [\mathbf{X}\mathbf{v}_1, \ldots, \mathbf{X}\mathbf{v}_t]\) com \(t\) componentes principais, \(1 \leq t \leq r\), onde \(\mathbf{V}_t\) é uma matriz com \(t\) colunas de \(\mathbf{V}\).

Finalmente, vamos fazer a regressão de \(\mathbf{Y}\) em \(\mathbf{Z}_t\) para obter os coeficientes \(\hat{\boldsymbol{\beta}}_{pcr}^{(t)} \in \mathbb{R}^t\), isto é, \(\hat{\boldsymbol{\beta}}_{pcr}^{(t)} = (\mathbf{Z}_t^T \mathbf{Z}_t)^{-1}\mathbf{Z}_t^T\mathbf{Y}\). Logo, a saída estimada (calculada no espaço reduzido) será dada por \(\hat{y}_{pcr}^{(t)} = \tilde{\mathbf{z}}^T \hat{\boldsymbol{\beta}}_{pcr}^{(t)}\), onde \(\tilde{\mathbf{z}} = \tilde{\mathbf{x}}\mathbf{V}_t \in \mathbb{R}^{t}\) é o vetor no espaço reduzido e \(\tilde{\mathbf{x}} \in \mathbb{R}^K\) é o vetor no espaço original do problema. Alternativamente, podemos lidar com o problema no espaço original. Neste caso, sendo \(\hat{\boldsymbol{\beta}}_{pcr}^{*(t)} \in \mathbb{R}^K\) os parâmetros no espaço original, temos que os coeficientes serão dados por \(\hat{\boldsymbol{\beta}}_{pcr}^{*(t)} = \mathbf{V}_t \hat{\boldsymbol{\beta}}_{pcr}^{(t)}\) e a saída estimada será dada por \(\hat{y}_{pcr}^{*(t)} = \tilde{\mathbf{x}}^T \hat{\boldsymbol{\beta}}_{pcr}^{*(t)}\). Os erros de ambas abordagens são iguais.

Metodologia

Vamos exemplificar o exposto usando o exercício 5.10 de [1, p. 156-157]. Os resultados serão mostrados passo-a-passo, mas o código-fonte completo encontra-se no final deste post. Usaremos a base de dados "detergent" (detergentyx.txt, que pode ser obtida em Master in Statistics Course). Esta base de dados possui 5 variáveis de saída, \(y1, \ldots, y5\), correspondendo a 4 compostos em uma solução aquosa (a quinta variável é a quantidade de água na solução), sendo que a soma destes 4 é igual a 1. As variáveis de entrada, \(x1, \ldots, x1168\), correspondem a valores de frequências espectrais gravadas como absorbâncias em \(K=1168\) frequências igualmente espaçadas em um intervalo \(3100-759 cm^{-1}\). Os dados consistem de \(N=12\) amostras de preparação de detergente. A questão pede para fazer o gráfico dos 12 espectros da absorbâncias e aplicar PCR nos dados usando cada uma das 4 primeiras variáveis de \(Y\) em separado.

A solução contará com os seguintes passos:

  1. Será lido o arquivo contendo os dados.
  2. Será feito um gráfico dos autovalores e uma medida de ajuste para auxiliar na verificação dos resultados. A medida de ajuste do modelo (chamada de "goodness-of-fit measure") é igual a \(\frac{\lambda_{t+1} + \cdots \lambda_{K}}{\lambda_{1} + \cdots \lambda_{K}}\), proposto em [1, p. 201], e corresponde a proporção que as variáveis de entrada são explicadas pelas \(r-t\) componentes principais.
  3. Serão feitas diversas simulações para obter os modelos PCR para cada saída \(y\) variando a quantidade de componentes principais. A comparação das saídas para as diferentes quantidades de componentes principais será feita por meio da visualização da saída calculada e por meio do cálculo do MSE e do \(R^2\) não-centrado, cuja fórmula pode ser obtida em qualquer livro de econometrica, como em [3, p. 19-20].

Existem diversos métodos para determinar qual a quantidade de componentes ideais. Por exemplo, [1, p. 132 e 208] afirma que a regra de Kaiser sugere que mantenha apenas as componentes principais cujos autovalores sejam maiores do que um. Apesar desta regra ser popular, ela é bastante controversa, pois existe evidência de que o corte em um é muito elevado. Alternativamente, podemos considerar o valor de corte como 0.7. Existe também a possibilidade de determinar \(t\) por validação cruzada.

+1 voto
respondida Jun 17, 2016 por Saulo (426 pontos)  

Resultados

Os resultados estão divididos em duas partes. A primeira estima os modelos PCR calculando explicitamente os autovalores e autovetores. A segunda utiliza uma biblioteca para fazer isto.

Usando a função SVD

  • Leitura do arquivo de dados:

    import pandas as pd
    import copy
    detergent_df = pd.read_table('detergentyx.txt')
    
    columns_x = list()
    for column in detergent_df.columns:
        if('x' in column):
            columns_x.append(column)
    
    X_df = copy.deepcopy(detergent_df[columns_x])
    N, K = X_df.shape
    
  • Plotando os espectros das variáveis, conforme solicitado na questão:

    fig, ax = plt.subplots()
    ax.plot(np.linspace(1, K, K), X_df.transpose())
    ax.set_xlim(1, K)
    plt.xlabel("Frequencia")
    plt.ylabel("Absorbancia")
    plt.show()
    

A imagem será apresentada aqui.

  • Calculando o posto, que mostra que posto(\(\mathbf{X}\))=12:

    r = np.linalg.matrix_rank(X_df)
    print "Posto(X) = ", r
    
  • Considerando \(t=r=12\), vamos mostrar e plotar os autovalores e a medida de ajuste do modelo:

    U, s, V = np.linalg.svd(X_df, full_matrices=False)
    
    autovalores = np.sqrt(s[0:r])
    figure_handle = 1
    
    fitness_measure = np.zeros(len(autovalores))
    print " "*17, "medida de\n", " i   autovalor    ajuste"
    for i in range(len(autovalores)):
    fitness_measure[i] = np.sum(autovalores[(i+1):len(autovalores)]) / np.sum(autovalores)
    print "{:2d} {:12.8f} {:12.8f}".format(i+1, autovalores[i], fitness_measure[i])
    
    import matplotlib.pyplot as plt
    plt.figure(figure_handle)
    plt.clf()
    plt.plot(range(1, len(autovalores)+1), autovalores, 'o:')
    #    plt.axes().set_yscale('log')
    plt.axes().set_ylim(0, 1.01*np.max(autovalores))
    plt.axes().set_xlim(1, r)
    plt.xlabel('i')
    plt.ylabel(r'$\lambda_i$')
    
    plt.figure(figure_handle+1)
    plt.clf()
    plt.plot(range(1, len(fitness_measure)+1), fitness_measure, 'o:')
    plt.axes().set_xlim(1, r)
    plt.xlabel('t')
    plt.ylabel(r'$\frac{\lambda_{t+1} + \cdots \lambda_{K}}{\lambda_{1} + \cdots \lambda_{K}}$')
    

Os autovalores e a medida de ajuste do modelo são dados por:

              medida de
 i   autovalor    ajuste
 1   8.05042364   0.45694811
 2   1.62657140   0.34722560
 3   1.40767214   0.25226923
 4   1.08235694   0.17925742
 5   0.91464375   0.11755892
 6   0.36440713   0.09297736
 7   0.28433771   0.07379699
 8   0.26350662   0.05602181
 9   0.23356801   0.04026617
10   0.20819034   0.02622242
11   0.19918446   0.01278617
12   0.18954736   0.00000000

Se considerarmos a regra de Kaiser, no mínimo devemos ter \(t=4\) ou \(t=5\) componentes para termos uma boa estimativa, dependendo se considerarmos o corte como sendo um ou 0.7.

O gráfico dos autovalores é:
A imagem será apresentada aqui.

O gráfico da medida de ajuste é:
A imagem será apresentada aqui.

  • Agora vamos calcular os autovalores e autovetores usando a função svd da biblioteca numpy.linalg. Por exemplo, considerando \(t=8\), \(\boldsymbol{\Lambda}\), \(\mathbf{V}\) e \(\mathbf{Z}_t\) podem ser calculados por:

    t = 8
    U, s, V = np.linalg.svd(X_df, full_matrices=False)
    Lambda = np.diag(np.sqrt(s[0:t]))
    V = V.transpose()[:,0:t]
    Zt = np.matmul(X_df, V)
    

E, para cada saída \(\mathbf{y}\), podemos calcular \(\hat{\boldsymbol{\beta}}_{pcr}^{(t)}\), \(\hat{\boldsymbol{\beta}}_{pcr}^{*(t)}\), \(\hat{\mathbf{Y}}_{pcr}^{(t)}\) e \(\hat{\mathbf{Y}}_{pcr}^{*(t)}\):

for saida in range(4):
    ystr = 'y'+str(saida+1)
    y = detergent_df[ystr]

    beta_pcr = np.matmul(np.matmul(np.linalg.inv(np.matmul(Zt.transpose(), Zt)), Zt.transpose()), y)
    beta_pcr_estrela = np.matmul(V[:,0:t], beta_pcr[0:t])

    result_df = pd.DataFrame()
    result_df[ystr+'_hat'] = np.matmul(Zt, beta_pcr)
    result_df[ystr+'*_hat'] = np.matmul(X_df, beta_pcr_estrela)

    erro_pcr = y - result_df[ystr+'_hat']

    result_df[ystr] = y
    result_df[ystr] = y
    R2 = np.dot(result_df[ystr+'_hat'], result_df[ystr+'_hat']) / np.dot(y, y)
    MSE = np.sum(np.square(erro_pcr))/(float(N)-1.0)

    print "Saida: ", ystr
    print result_df
    print ""
    print "R2 = ", R2
    print "MSE = ", MSE

Fazendo agora o mesmo, porém, para cada valor de \(t\) da lista \(\{1, 2, 4, 8, 12\}\), temos:

    lista_t = list([1, 2, 4, 8, 12])
    total_t = len(lista_t)

    tabela_beta_pcr = dict()
    tabela_beta_pcr_estrela = dict()
    tabela_y_calculado_pcr = dict()
    tabela_y_calculado_pcr_estrela = dict()
    tabela_erro_pcr = dict()
    tabela_erro_pcr_estrela = dict()
    tabela_beta_R2 = dict()
    tabela_beta_MSE = dict()

    conferencia = 0.0

    total_saidas = 4

    for indice_t in range(total_t):
    t = lista_t[indice_t]

    U, s, V = np.linalg.svd(X_df, full_matrices=False)
    Lambda = np.diag(np.sqrt(s[0:t]))
    V = V.transpose()[:,0:t]

    Zt = np.matmul(X_df, V)

    for saida in range(total_saidas):
        ystr = 'y'+str(saida+1)
        y = detergent_df[ystr]

        if (not tabela_beta_pcr.has_key(t)):
            tabela_beta_pcr[t] = np.zeros((t, total_saidas))
            tabela_beta_pcr_estrela[t] = np.zeros((K, total_saidas))
            tabela_y_calculado_pcr[t] = np.zeros((N, total_saidas))
            tabela_y_calculado_pcr_estrela[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr_estrela[t] = np.zeros((N, total_saidas))

        beta_pcr = np.matmul(np.matmul(np.linalg.inv(np.matmul(Zt.transpose(), Zt)), Zt.transpose()), y)
        tabela_beta_pcr[t][:,saida] = beta_pcr
        y_calculado_pcr = np.matmul(Zt, beta_pcr)
        tabela_y_calculado_pcr[t][:,saida] = y_calculado_pcr
        erro_pcr = y - y_calculado_pcr
        tabela_erro_pcr[t][:,saida] = erro_pcr

        beta_pcr_estrela = np.matmul(V[:,0:t], beta_pcr[0:t])
        tabela_beta_pcr_estrela[t][:,saida] = beta_pcr_estrela
        y_calculado_pcr_estrela = np.matmul(X_df, beta_pcr_estrela)
        tabela_y_calculado_pcr_estrela[t][:,saida] = y_calculado_pcr_estrela
        tabela_erro_pcr_estrela[t][:,saida] = y - y_calculado_pcr_estrela

        R2 = np.dot(y_calculado_pcr, y_calculado_pcr) / np.dot(y, y)
        tabela_beta_R2[(saida,t)] = R2
        MSE = np.sum(np.square(erro_pcr))/(float(N)-1.0)
        tabela_beta_MSE[(saida,t)] = MSE

        conferencia += np.sum(np.square(y_calculado_pcr-y_calculado_pcr_estrela))

    print "Conferencia se as duas saidas calculadas sao iguais = ", conferencia

    for saida in range(total_saidas):
    ystr = 'y'+str(saida+1)
    y = detergent_df[ystr]
    print "\n"
    print "Saida: y%d" % (saida+1)
    print ""

    print " "*10, "y",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print "  ycalc(t={:2d})".format(t),
    print ""

    for k in range(len(y)):
        print "{:12.6f}".format(y[k]),
        for indice_t in range(total_t):
            t = lista_t[indice_t]
            y_calculado_pcr = tabela_y_calculado_pcr[t]
            print " {:12.6f}".format(y_calculado_pcr[k,saida]),
        print ""

    print ""
    print " "*8, "R^2",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_R2[(saida,t)]),
    print ""

    print " "*8, "MSE",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_MSE[(saida,t)]),
    print ""

    print "\n"

Obtemos o seguinte resultado. Considerando \(y1\), a saída calculada \(\hat{y}_{pcr}^{(t)}\), \(R^2\) e \(MSE\) para cada valor de \(t\) são dados por:

Saida: y1

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
   20.002000     19.403282     16.778657     19.874215     19.977004     20.002000 
   14.918000     19.714532     19.422673     17.188864     14.913503     14.918000 
   19.954000     18.866698     19.457761     19.170299     19.964602     19.954000 
   24.645000     19.443239     23.146123     26.475495     24.647522     24.645000 
   24.913000     20.280792     20.003790     23.353754     24.951649     24.913000 
   14.868000     18.592264     15.818994     12.669047     14.874155     14.868000 
   25.197000     19.737467     22.865253     23.580380     25.183339     25.197000 
   15.124000     20.184657     16.921716     17.472600     15.139870     15.124000 
   19.917000     20.466943     23.718642     19.997403     19.972368     19.917000 
   21.392000     19.485616     19.248306     19.546207     21.338097     21.392000 
   17.652000     20.086186     20.524888     18.596645     17.559838     17.652000 
   16.685000     19.094982     17.214837     17.169972     16.744908     16.685000 

     R^2    0.96692614    0.98023551    0.99391454    0.99999565    1.00000000 
     MSE   14.36545222    8.58459839    2.64318582    0.00188981    0.00000000 

Para \(y2\), temos:

Saida: y2

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
    6.989000      9.906081      8.326370      7.234526      7.023600      6.989000 
   13.096000     10.064985      9.889321     10.740789     12.982780     13.096000 
    9.906000      9.632135      9.987884     10.471424      9.915938      9.906000 
   13.090000      9.926481     12.155174     11.435174     13.093331     13.090000 
    7.080000     10.354082     10.187360      8.558967      6.974023      7.080000 
    6.894000      9.492027      7.822849      9.281619      6.894732      6.894000 
   10.073000     10.076695     11.959248     11.662922     10.133738     10.073000 
   10.084000     10.305002      8.341100      7.844055     10.179133     10.084000 
   13.220000     10.449119     12.406253     13.316493     13.259940     13.220000 
    7.908000      9.948116      9.805283      9.644362      7.911111      7.908000 
   12.119000     10.254728     10.518775     11.014240     12.097933     12.119000 
    9.487000      9.748682      8.617060      8.911060      9.480203      9.487000 

     R^2    0.95202347    0.97023637    0.97732850    0.99996822    1.00000000 
     MSE    5.51648330    3.42231028    2.60683583    0.00365380    0.00000000 

Para \(y3\), temos:

Saida: y3

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
    3.020000      5.417296      4.661054      2.899539      3.035241      3.020000 
    5.294000      5.504195      5.420101      6.601911      5.296262      5.294000 
    5.389000      5.267484      5.437789      5.064400      5.386964      5.389000 
    3.008000      5.428452      6.495375      3.962439      3.003286      3.008000 
    5.242000      5.662292      5.582479      4.304030      5.235249      5.242000 
    7.116000      5.190864      4.391793      5.695634      7.113305      7.116000 
    6.951000      5.510599      6.411817      6.066832      6.956447      6.951000 
    3.456000      5.635452      4.695290      4.807203      3.444451      3.456000 
    9.382000      5.714265      6.651187      9.334013      9.373149      9.382000 
    6.283000      5.440283      5.371906      5.280861      6.281637      6.283000 
    6.425000      5.607959      5.734364      7.091084      6.447449      6.425000 
    4.065000      5.331220      4.789488      4.424145      4.057625      4.065000 

     R^2    0.90272263    0.91595660    0.97536964    0.99999720    1.00000000 
     MSE    3.52776788    3.04783733    0.89322100    0.00010165    0.00000000 

Para \(y4\), temos:

Saida: y4

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
    6.346000      6.060365      6.221583      6.310676      6.320264      6.346000 
    6.025000      6.157580      6.175507      5.813597      6.046456      6.025000 
    4.004000      5.892770      5.856463      4.059188      4.007786      4.004000 
    4.004000      6.072845      5.845395      3.814380      4.008275      4.004000 
    8.363000      6.334444      6.351459      8.539617      8.404535      8.363000 
    4.094000      5.807054      5.977402      4.258808      4.097314      4.094000 
    6.066000      6.164743      5.972619      6.199614      6.039748      6.066000 
    8.192000      6.304417      6.504844      7.937963      8.181386      8.192000 
    7.962000      6.392586      6.192850      7.969043      7.975963      7.962000 
    6.164000      6.086081      6.100658      6.370681      6.149460      6.164000 
    7.172000      6.273661      6.246714      7.055653      7.136879      7.172000 
    4.697000      5.964071      6.079559      4.776006      4.720735      4.697000 

     R^2    0.95053959    0.95104562    0.99939209    0.99998756    1.00000000 
     MSE    2.13187841    2.11006698    0.02620254    0.00053604    0.00000000 

Como se pode perceber, à medida que acrescentamos as componentes principais, \(R^2 \rightarrow 1\) e \(MSE \rightarrow 0\).

+1 voto
respondida Jun 17, 2016 por Saulo (426 pontos)  

Usando a biblioteca sklearn

Podemos usar também a análise de componentes principais de uma biblioteca. Usaremos a classe sklearn.decomposition.PCA da biblioteca sklearn.

O gráfico do percentual da variância explicada associada à cada componente principal é dado por:

    from sklearn import decomposition
    pca = decomposition.PCA(n_components=r)
    pca.fit(X_df)
    plt.figure(3)
    plt.clf()
    plt.plot(pca.explained_variance_, linewidth=2)
    plt.axis('tight')
    plt.xlabel('Componente i')
    plt.ylabel('Variancia explicada')

A imagem será apresentada aqui.

Podemos perceber que para \(t>4\), acrescenta-se muito pouco na explicação da variância.

Para efetuar os cálculos usando a biblioteca, apenas precisamos alterar a forma de cálculo de \(\mathbf{V}\). Ou seja, para um determinado \(t\), vamos fazer o seguinte:

    pca = decomposition.PCA(n_components=t)
    pca.fit(X_df)
    Lambda = np.diag(pca.explained_variance_[0:t])
    V = pca.components_.transpose()[:,0:t]

Logo, para fazer as mesmas estatísticas para \(t \in \{1, 2, 4, 8, 12\}\), para cada saída, fazemos:

    from sklearn import decomposition

    lista_t = list([1, 2, 4, 8, 12])
    total_t = len(lista_t)

    tabela_beta_pcr = dict()
    tabela_beta_pcr_estrela = dict()
    tabela_y_calculado_pcr = dict()
    tabela_y_calculado_pcr_estrela = dict()
    tabela_erro_pcr = dict()
    tabela_erro_pcr_estrela = dict()
    tabela_beta_R2 = dict()
    tabela_beta_MSE = dict()

    conferencia = 0.0

    total_saidas = 4

    for indice_t in range(total_t):
    t = lista_t[indice_t]

    pca = decomposition.PCA(n_components=t)
    pca.fit(X_df)
    Lambda = np.diag(pca.explained_variance_[0:t])
    V = pca.components_.transpose()[:,0:t]

    Zt = np.matmul(X_df, V)

    for saida in range(total_saidas):
        ystr = 'y'+str(saida+1)
        y = detergent_df[ystr]

        if (not tabela_beta_pcr.has_key(t)):
            tabela_beta_pcr[t] = np.zeros((t, total_saidas))
            tabela_beta_pcr_estrela[t] = np.zeros((K, total_saidas))
            tabela_y_calculado_pcr[t] = np.zeros((N, total_saidas))
            tabela_y_calculado_pcr_estrela[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr_estrela[t] = np.zeros((N, total_saidas))

        beta_pcr = np.matmul(np.matmul(np.linalg.inv(np.matmul(Zt.transpose(), Zt)), Zt.transpose()), y)
        tabela_beta_pcr[t][:,saida] = beta_pcr
        y_calculado_pcr = np.matmul(Zt, beta_pcr)
        tabela_y_calculado_pcr[t][:,saida] = y_calculado_pcr
        erro_pcr = y - y_calculado_pcr
        tabela_erro_pcr[t][:,saida] = erro_pcr

        beta_pcr_estrela = np.matmul(V[:,0:t], beta_pcr[0:t])
        tabela_beta_pcr_estrela[t][:,saida] = beta_pcr_estrela
        y_calculado_pcr_estrela = np.matmul(X_df, beta_pcr_estrela)
        tabela_y_calculado_pcr_estrela[t][:,saida] = y_calculado_pcr_estrela
        tabela_erro_pcr_estrela[t][:,saida] = y - y_calculado_pcr_estrela

        R2 = np.dot(y_calculado_pcr, y_calculado_pcr) / np.dot(y, y)
        tabela_beta_R2[(saida,t)] = R2
        MSE = np.sum(np.square(erro_pcr))/(float(N)-1.0)
        tabela_beta_MSE[(saida,t)] = MSE

        conferencia += np.sum(np.square(y_calculado_pcr-y_calculado_pcr_estrela))

    print "Conferencia se as duas saidas calculadas sao iguais = ", conferencia

    for saida in range(total_saidas):
    ystr = 'y'+str(saida+1)
    y = detergent_df[ystr]
    print "\n"
    print "Saida: y%d" % (saida+1)
    print ""

    print " "*10, "y",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print "  ycalc(t={:2d})".format(t),
    print ""

    for k in range(len(y)):
        print "{:12.6f}".format(y[k]),
        for indice_t in range(total_t):
            t = lista_t[indice_t]
            y_calculado_pcr = tabela_y_calculado_pcr[t]
            print " {:12.6f}".format(y_calculado_pcr[k,saida]),
        print ""

    print ""
    print " "*8, "R^2",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_R2[(saida,t)]),
    print ""

    print " "*8, "MSE",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_MSE[(saida,t)]),
    print ""

    print "\n"

Finalmente, usando a biblioteca, temos os seguintes resultados. Para \(y1\), temos:

Saida: y1

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
   20.002000     18.291587     16.917360     19.350422     19.806822     20.002000 
   14.918000     19.451146     19.343660     17.236763     14.884519     14.918000 
   19.954000     18.615550     19.102209     18.096622     19.817120     19.954000 
   24.645000     20.858656     22.776147     25.050709     24.733301     24.645000 
   24.913000     21.048151     20.480461     24.212126     25.014848     24.913000 
   14.868000     16.534639     15.476114     11.754349     14.956918     14.868000 
   25.197000     21.471271     22.849294     23.771639     25.184696     25.197000 
   15.124000     19.125846     17.298946     18.096178     15.194420     15.124000 
   19.917000     22.621543     23.917122     21.726103     19.926762     19.917000 
   21.392000     19.512870     19.300331     19.599649     21.361838     21.392000 
   17.652000     20.579757     20.645965     19.351370     17.577765     17.652000 
   16.685000     17.780957     17.019062     16.235083     16.804594     16.685000 

     R^2    0.97743059    0.98052857    0.99161896    0.99997699    1.00000000 
     MSE    9.80290020    8.45731070    3.64026095    0.00999610    0.00000000 

Para \(y2\), temos:

Saida: y2

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
    6.989000      9.336359      8.215602      7.554280      7.090948      6.989000 
   13.096000      9.928218      9.840558     11.370014     12.999669     13.096000 
    9.906000      9.501715      9.898612     11.035295      9.982041      9.906000 
   13.090000     10.646637     12.210455     12.942508     13.056189     13.090000 
    7.080000     10.743359     10.280377      7.493094      6.946871      7.080000 
    6.894000      8.439580      7.576296      9.249922      6.857592      6.894000 
   10.073000     10.959327     12.083179     11.071095     10.132934     10.073000 
   10.084000      9.762179      8.272244      8.093781     10.159857     10.084000 
   13.220000     11.546447     12.603061     12.098041     13.278584     13.220000 
    7.908000      9.959724      9.786387      9.086000      7.909486      7.908000 
   12.119000     10.504282     10.558278     10.795871     12.081389     12.119000 
    9.487000      9.075724      8.454358      9.706089      9.452404      9.487000 

     R^2    0.96192054    0.96970421    0.98440958    0.99995174    1.00000000 
     MSE    4.37848848    3.48349885    1.79263202    0.00554903    0.00000000 

Para \(y3\), temos:

Saida: y3

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
    3.020000      5.107420      4.682583      2.836312      3.008149      3.020000 
    5.294000      5.431194      5.397965      5.814788      5.300216      5.294000 
    5.389000      5.197877      5.348326      5.044861      5.375010      5.389000 
    3.008000      5.824203      6.416987      3.022932      3.013031      3.008000 
    5.242000      5.877114      5.701615      5.060157      5.244201      5.242000 
    7.116000      4.616840      4.289602      6.295955      7.121185      7.116000 
    6.951000      5.995258      6.421268      6.661929      6.954616      6.951000 
    3.456000      5.340363      4.775584      4.123396      3.455755      3.456000 
    9.382000      6.316440      6.716963      9.739111      9.366765      9.382000 
    6.283000      5.448429      5.382724      5.920300      6.294399      6.283000 
    6.425000      5.746328      5.766795      6.890349      6.437905      6.425000 
    4.065000      4.964840      4.729304      4.041146      4.059545      4.065000 

     R^2    0.91270937    0.91625548    0.99464966    0.99999748    1.00000000 
     MSE    3.16559807    3.03699848    0.19403031    0.00009121    0.00000000 

Para \(y4\), temos:

Saida: y4

       y   ycalc(t= 1)   ycalc(t= 2)   ycalc(t= 4)   ycalc(t= 8)   ycalc(t=12) 
    6.346000      5.693057      6.577745      6.544830      6.366096      6.346000 
    6.025000      6.053956      6.123153      5.186635      6.041223      6.025000 
    4.004000      5.793886      5.480589      4.623406      4.030849      4.004000 
    4.004000      6.492029      5.257603      3.864918      3.990760      4.004000 
    8.363000      6.551007      6.916469      8.654235      8.387583      8.363000 
    4.094000      5.146226      5.827674      5.210300      4.082835      4.094000 
    6.066000      6.682698      5.795567      6.551026      6.044054      6.066000 
    8.192000      5.952710      7.128816      7.079472      8.164867      8.192000 
    7.962000      7.040708      6.206651      7.346316      7.987092      7.962000 
    6.164000      6.073167      6.209994      6.830302      6.130318      6.164000 
    7.172000      6.405224      6.362602      6.495385      7.148770      7.172000 
    4.697000      5.534128      6.024614      4.987641      4.714898      4.697000 

     R^2    0.95412197    0.96706010    0.98879080    0.99998706    1.00000000 
     MSE    1.97746773    1.41979948    0.48314694    0.00055764    0.00000000 

Chegamos a mesma conclusão, à medida que acrescentamos as componentes principais, \(R^2 \rightarrow 1\) e \(MSE \rightarrow 0\).

Comparação

Não é possível determinar qual método apresenta os melhores resultados. Em algumas situações (para um determinado valor de \(t\) e para alguma saída), a função "numpy.linalg.svd" apresenta melhores resultados (maior \(R^2\) e menor \(MSE\)), mas em outras a classe "sklearn.decomposition.PCA" apresenta melhores resultados.

+1 voto
respondida Jun 17, 2016 por Saulo (426 pontos)  

Referências

[ 1 ] Izenman, Alan Julian. "Modern multivariate statistical techniques". Springer, 2008.
[ 2 ] Bishop, Christopher M. "Pattern Recognition and Machine Learning". Springer, 2006.
[ 3 ] Fumio, Hayashi. "Econometrics". Princeton University Press, 2000.

Anexo: Código-fonte completo

    import numpy as np
    import matplotlib.pyplot as plt

    # ****************************************************************************
    # Leitura dos dados
    # ****************************************************************************
    import pandas as pd
    import copy
    detergent_df = pd.read_table('detergentyx.txt')

    columns_x = list()
    for column in detergent_df.columns:
    if('x' in column):
        columns_x.append(column)

    X_df = copy.deepcopy(detergent_df[columns_x])
    N, K = X_df.shape

    # ****************************************************************************
    # Plotando os espectros
    # ****************************************************************************
    fig, ax = plt.subplots()
    ax.plot(np.linspace(1, K, K), X_df.transpose())
    ax.set_xlim(1, K)
    plt.xlabel("Frequencia")
    plt.ylabel("Absorbancia")
    plt.show()


    # ****************************************************************************
    # Calculando o posto(X)
    # ****************************************************************************
    r = np.linalg.matrix_rank(X_df)
    print "Posto(X) = ", r
    print ""

    # ****************************************************************************
    # Mostrando e plotando os autovalores e a metrica
    # ****************************************************************************
    U, s, V = np.linalg.svd(X_df, full_matrices=False)

    autovalores = np.sqrt(s[0:r])
    figure_handle = 1

    fitness_measure = np.zeros(len(autovalores))
    print " "*17, "medida de\n", " i   autovalor    ajuste"
    for i in range(len(autovalores)):
    fitness_measure[i] = np.sum(autovalores[(i+1):len(autovalores)]) / np.sum(autovalores)
    print "{:2d} {:12.8f} {:12.8f}".format(i+1, autovalores[i], fitness_measure[i])

    import matplotlib.pyplot as plt
    plt.figure(figure_handle)
    plt.clf()
    plt.plot(range(1, len(autovalores)+1), autovalores, 'o:')
    #    plt.axes().set_yscale('log')
    plt.axes().set_ylim(0, 1.01*np.max(autovalores))
    plt.axes().set_xlim(1, r)
    plt.xlabel('i')
    plt.ylabel(r'$\lambda_i$')

    plt.figure(figure_handle+1)
    plt.clf()
    plt.plot(range(1, len(fitness_measure)+1), fitness_measure, 'o:')
    #    plt.axes().set_ylim(0, 1.01*np.max(autovalores))
    plt.axes().set_xlim(1, r)
    plt.xlabel('t')
    plt.ylabel(r'$\frac{\lambda_{t+1} + \cdots \lambda_{K}}{\lambda_{1} + \cdots \lambda_{K}}$')


    # ****************************************************************************
    # Calculando para t = 8
    # ****************************************************************************
    print "***********************************************************************"
    print "Resultado para t = 8"
    print "***********************************************************************"

    t = 8
    U, s, V = np.linalg.svd(X_df, full_matrices=False)
    Lambda = np.diag(np.sqrt(s[0:t]))
    V = V.transpose()[:,0:t]
    Zt = np.matmul(X_df, V)

    for saida in range(4):
    ystr = 'y'+str(saida+1)
    y = detergent_df[ystr]

    beta_pcr = np.matmul(np.matmul(np.linalg.inv(np.matmul(Zt.transpose(), Zt)), Zt.transpose()), y)
    beta_pcr_estrela = np.matmul(V[:,0:t], beta_pcr[0:t])

    result_df = pd.DataFrame()
    result_df[ystr+'_hat'] = np.matmul(Zt, beta_pcr)
    result_df[ystr+'*_hat'] = np.matmul(X_df, beta_pcr_estrela)

    erro_pcr = y - result_df[ystr+'_hat']

    result_df[ystr] = y
    result_df[ystr] = y
    R2 = np.dot(result_df[ystr+'_hat'], result_df[ystr+'_hat']) / np.dot(y, y)
    MSE = np.sum(np.square(erro_pcr))/(float(N)-1.0)

    print "Saida: ", ystr
    print result_df
    print ""
    print "R2 = ", R2
    print "MSE = ", MSE


    # ****************************************************************************
    # Calculando para uma lista de t
    # Usando SVD
    # http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.svd.html
    # ****************************************************************************
    print "***********************************************************************"
    print "Resultados para PCS usando numpy.linalg.svd"
    print "***********************************************************************"

    lista_t = list([1, 2, 4, 8, 12])
    total_t = len(lista_t)

    tabela_beta_pcr = dict()
    tabela_beta_pcr_estrela = dict()
    tabela_y_calculado_pcr = dict()
    tabela_y_calculado_pcr_estrela = dict()
    tabela_erro_pcr = dict()
    tabela_erro_pcr_estrela = dict()
    tabela_beta_R2 = dict()
    tabela_beta_MSE = dict()

    conferencia = 0.0

    total_saidas = 4

    for indice_t in range(total_t):
    t = lista_t[indice_t]

    U, s, V = np.linalg.svd(X_df, full_matrices=False)
    Lambda = np.diag(np.sqrt(s[0:t]))
    V = V.transpose()[:,0:t]

    Zt = np.matmul(X_df, V)

    for saida in range(total_saidas):
        ystr = 'y'+str(saida+1)
        y = detergent_df[ystr]

        if (not tabela_beta_pcr.has_key(t)):
            tabela_beta_pcr[t] = np.zeros((t, total_saidas))
            tabela_beta_pcr_estrela[t] = np.zeros((K, total_saidas))
            tabela_y_calculado_pcr[t] = np.zeros((N, total_saidas))
            tabela_y_calculado_pcr_estrela[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr_estrela[t] = np.zeros((N, total_saidas))

        beta_pcr = np.matmul(np.matmul(np.linalg.inv(np.matmul(Zt.transpose(), Zt)), Zt.transpose()), y)
        tabela_beta_pcr[t][:,saida] = beta_pcr
        y_calculado_pcr = np.matmul(Zt, beta_pcr)
        tabela_y_calculado_pcr[t][:,saida] = y_calculado_pcr
        erro_pcr = y - y_calculado_pcr
        tabela_erro_pcr[t][:,saida] = erro_pcr

        beta_pcr_estrela = np.matmul(V[:,0:t], beta_pcr[0:t])
        tabela_beta_pcr_estrela[t][:,saida] = beta_pcr_estrela
        y_calculado_pcr_estrela = np.matmul(X_df, beta_pcr_estrela)
        tabela_y_calculado_pcr_estrela[t][:,saida] = y_calculado_pcr_estrela
        tabela_erro_pcr_estrela[t][:,saida] = y - y_calculado_pcr_estrela

        R2 = np.dot(y_calculado_pcr, y_calculado_pcr) / np.dot(y, y)
        tabela_beta_R2[(saida,t)] = R2
        MSE = np.sum(np.square(erro_pcr))/(float(N)-1.0)
        tabela_beta_MSE[(saida,t)] = MSE

        conferencia += np.sum(np.square(y_calculado_pcr-y_calculado_pcr_estrela))

    print "Conferencia se as duas saidas calculadas sao iguais = ", conferencia

    for saida in range(total_saidas):
    ystr = 'y'+str(saida+1)
    y = detergent_df[ystr]
    print "\n"
    print "Saida: y%d" % (saida+1)
    print ""

    print " "*10, "y",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print "  ycalc(t={:2d})".format(t),
    print ""

    for k in range(len(y)):
        print "{:12.6f}".format(y[k]),
        for indice_t in range(total_t):
            t = lista_t[indice_t]
            y_calculado_pcr = tabela_y_calculado_pcr[t]
            print " {:12.6f}".format(y_calculado_pcr[k,saida]),
        print ""

    print ""
    print " "*8, "R^2",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_R2[(saida,t)]),
    print ""

    print " "*8, "MSE",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_MSE[(saida,t)]),
    print ""

    print "\n"

    # ****************************************************************************
    # Plotando o grafico da quantidade de variancia explicada para cada componente
    # ****************************************************************************
    from sklearn import decomposition
    pca = decomposition.PCA(n_components=r)
    pca.fit(X_df)
    plt.figure(3)
    plt.clf()
    plt.plot(pca.explained_variance_, linewidth=2)
    plt.axis('tight')
    plt.xlabel('Componente i')
    plt.ylabel('Variancia explicada')


    # ****************************************************************************
    # Calculando para uma lista de t
    # Usando a biblioteca sklearn
    # http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html
    # ****************************************************************************
    print "***********************************************************************"
    print "Resultados para PCR usando sklearn.decomposition.PCA"
    print "***********************************************************************"

    from sklearn import decomposition

    lista_t = list([1, 2, 4, 8, 12])
    total_t = len(lista_t)

    tabela_beta_pcr = dict()
    tabela_beta_pcr_estrela = dict()
    tabela_y_calculado_pcr = dict()
    tabela_y_calculado_pcr_estrela = dict()
    tabela_erro_pcr = dict()
    tabela_erro_pcr_estrela = dict()
    tabela_beta_R2 = dict()
    tabela_beta_MSE = dict()

    conferencia = 0.0

    total_saidas = 4

    for indice_t in range(total_t):
    t = lista_t[indice_t]

    pca = decomposition.PCA(n_components=t)
    pca.fit(X_df)
    Lambda = np.diag(pca.explained_variance_[0:t])
    V = pca.components_.transpose()[:,0:t]

    Zt = np.matmul(X_df, V)

    for saida in range(total_saidas):
        ystr = 'y'+str(saida+1)
        y = detergent_df[ystr]

        if (not tabela_beta_pcr.has_key(t)):
            tabela_beta_pcr[t] = np.zeros((t, total_saidas))
            tabela_beta_pcr_estrela[t] = np.zeros((K, total_saidas))
            tabela_y_calculado_pcr[t] = np.zeros((N, total_saidas))
            tabela_y_calculado_pcr_estrela[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr[t] = np.zeros((N, total_saidas))
            tabela_erro_pcr_estrela[t] = np.zeros((N, total_saidas))

        beta_pcr = np.matmul(np.matmul(np.linalg.inv(np.matmul(Zt.transpose(), Zt)), Zt.transpose()), y)
        tabela_beta_pcr[t][:,saida] = beta_pcr
        y_calculado_pcr = np.matmul(Zt, beta_pcr)
        tabela_y_calculado_pcr[t][:,saida] = y_calculado_pcr
        erro_pcr = y - y_calculado_pcr
        tabela_erro_pcr[t][:,saida] = erro_pcr

        beta_pcr_estrela = np.matmul(V[:,0:t], beta_pcr[0:t])
        tabela_beta_pcr_estrela[t][:,saida] = beta_pcr_estrela
        y_calculado_pcr_estrela = np.matmul(X_df, beta_pcr_estrela)
        tabela_y_calculado_pcr_estrela[t][:,saida] = y_calculado_pcr_estrela
        tabela_erro_pcr_estrela[t][:,saida] = y - y_calculado_pcr_estrela

        R2 = np.dot(y_calculado_pcr, y_calculado_pcr) / np.dot(y, y)
        tabela_beta_R2[(saida,t)] = R2
        MSE = np.sum(np.square(erro_pcr))/(float(N)-1.0)
        tabela_beta_MSE[(saida,t)] = MSE

        conferencia += np.sum(np.square(y_calculado_pcr-y_calculado_pcr_estrela))

    print "Conferencia se as duas saidas calculadas sao iguais = ", conferencia

    for saida in range(total_saidas):
    ystr = 'y'+str(saida+1)
    y = detergent_df[ystr]
    print "\n"
    print "Saida: y%d" % (saida+1)
    print ""

    print " "*10, "y",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print "  ycalc(t={:2d})".format(t),
    print ""

    for k in range(len(y)):
        print "{:12.6f}".format(y[k]),
        for indice_t in range(total_t):
            t = lista_t[indice_t]
            y_calculado_pcr = tabela_y_calculado_pcr[t]
            print " {:12.6f}".format(y_calculado_pcr[k,saida]),
        print ""

    print ""
    print " "*8, "R^2",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_R2[(saida,t)]),
    print ""

    print " "*8, "MSE",
    for indice_t in range(total_t):
        t = lista_t[indice_t]
        print " {:12.8f}".format(tabela_beta_MSE[(saida,t)]),
    print ""

    print "\n"
...