As manipulação de matrizes no Python são feitas comumente com o uso da biblioteca NumPy.

Vamos criar um script em Python com o cabeçalho:

#!/usr/bin/env python3

# -*- coding: utf8 -*-

E “carregar” a biblioteca NumPy com o comando import numpy as np.

A matriz “Y” pode ser criada de duas formas:

>>> Y_exp = np.array([[24], [40], [60], [70], [77], [86], [91], [86], [84]])
>>> Y_exp
array([[24],
       [40],
       [60],
       [70],
       [77],
       [86],
       [91],
       [86],
       [84]])

Mas também pode ser criada como uma lista (1 dimensão) e posteriormente transformada em uma matriz (2 dimensões) com o comando:

>>> Y = np.array([24, 40, 60, 70, 77, 86, 91, 86, 84])
>>> Y.shape
(9,)
>>> Y_exp = Y_exp.reshape(-1,1)
>>> Y_exp
array([[24],
       [40],
       [60],
       [70],
       [77],
       [86],
       [91],
       [86],
       [84]])
>>> Y_exp.shape
(9, 1)

Um exemplo de matriz de dados “X” pode ser criada com os comandos:

  >>> X = np.array([[1, 30], [1, 35], [1, 40], [1, 45], [1, 50], [1, 55], [1, 60], [1, 65], [1, 70]])
>>> X
array([[ 1, 30],
       [ 1, 35],
       [ 1, 40],
       [ 1, 45],
       [ 1, 50],
       [ 1, 55],
       [ 1, 60],
       [ 1, 65],
       [ 1, 70]])

Mas também pode ser “montada” com o comando:

X = np.array([[30, 35, 40, 45, 50, 55, 60, 65, 70]])
>>> X
array([[30, 35, 40, 45, 50, 55, 60, 65, 70]])
>>> X = X.T
>>> X
array([[30],
       [35],
       [40],
       [45],
       [50],
       [55],
       [60],
       [65],
       [70]])
>>> X = np.hstack([np.ones(X.shape), X])
>>> X
array([[ 1., 30.],
       [ 1., 35.],
       [ 1., 40.],
       [ 1., 45.],
       [ 1., 50.],
       [ 1., 55.],
       [ 1., 60.],
       [ 1., 65.],
       [ 1., 70.]])

E a sua transposta “X_t”:

 >>> X_t = X.T
>>> X_t
array([[ 1,  1,  1,  1,  1,  1,  1,  1,  1],
       [30, 35, 40, 45, 50, 55, 60, 65, 70]])

O produto matricial “X^tX” (X_t_prod_matricial_X) é obtido pelo comando:

>>> X_t_prod_matricial_X = np.dot(X_t,X)
>>> X_t_prod_matricial_X
array([[    9,   450],
       [  450, 24000]])

E a matriz inversa do produto de X^t por X é:

>>> inv_X_t_dot_X = np.linalg.inv(X_t_prod_matricial_X)
>>> inv_X_t_dot_X
array([[ 1.77777778e+00, -3.33333333e-02],
       [-3.33333333e-02,  6.66666667e-04]])

O segundo termo da equação S.2 ((X^tY) pode ser calculado pelo comando:

>>> X_t_Y_exp = np.dot(X_t,Y_exp)
>>> X_t_Y_exp
array([  618, 33180])
>>> X_t_Y_exp.shape
(2,)

Podemos observar que o resultado do produto matricial de X^t por Y é uma estrutura no formato de lista (unidimensinal). Para transformar em uma matriz de 2 linhas por 1 coluna, usamos o “método” reshape.

>>> X_t_Y_exp = X_t_Y_exp.reshape(-1,1)
array([[  618],
       [33180]])
>>> X_t_Y_exp.shape
(2, 1)

E finalmente a matriz “B” com os parâmetros da equação de regressão é calculada com o comando:

>>> B = np.dot(inv_X_t_dot_X, X_t_Y_exp)
>>> B
array([[-7.33333333],
       [ 1.52      ]])
>>> B.shape
(2, 1)

Portanto, de acordo com a equação S.2, as etapas para o cálculo da matriz de parâmetros “b” é representada na figura S.5

Figura S.5. Etapas para o cálculo matricial dos parâmetros do modelo (b₀ e b₁) aplicando as matrizes da figura S.4 na equação S.2.

A matriz “B” representa os valores de b₀ e b₁ que são as “estimativas” de β₀ e β₁, respectivamente.

E o modelo empírico, descrito na equação S.1, para descrever a variação do rendimento da reação com a temperatura, pode ser escrito como na equação S.3.

A equação S.3 também pode ser representada na forma matricial (Ŷ = Xb) como mostra a figura S.4.

Equação S.4. Equação matricial para o cálculo da matriz dos “rendimentos previstos” (Ŷ) pela multiplicação da matriz das temperaturas (X) pela matriz dos parâmetros (b).

E a implementação do cálculo com o Python através do comando (representando a matriz Ŷ pela variável Y_est) :

>>> Y_est = np.dot(X,B)
>>> Y_est
array([[38.26666667],
       [45.86666667],
       [53.46666667],
       [61.06666667],
       [68.66666667],
       [76.26666667],
       [83.86666667],
       [91.46666667],
       [99.06666667]])

A diferença entre os valores previstos (Y_est) e os valores experimentais (Y_exp) correspondem aos “resíduos” (E).

>>> E = Y_exp - Y_est
>>> E
array([[-14.26],
       [ -5.86],
       [  6.54],
       [  8.94],
       [  8.34],
       [  9.74],
       [  7.14],
       [ -5.46],
       [-15.06]])

Para avaliar a “qualidade do ajuste” de qualquer modelo deve-se verificar se os “resíduos” gerados pelo modelo são pequenos.

O desvio de uma resposta individual (y_i) em relação à média das respostas observadas (y_medio) (Figura S.6) é formado pelos componentes:

(ŷ − y_medio) - desvio de previsão feito pelo modelo

(y_i − ŷ) - diferença entre o valor observado e o valor previsto (deve ser pequena em um bom modelo)

(y_i − y_medio) = (y_i − ŷ) + (ŷ − y_medio)

Figura S.6. Decomposição do desvio de uma observação em relação à média global. (Fonte: Como fazer experimentos, 2001)

Elevando a equação ao quadrado e fazendo o somatório em todos os pontos:

Σ(y_i − y_medio)² = Σ[(y_i − ŷ) + (ŷ − y_medio)]²

chega-se à equação:

Σ(y_i − y_medio)² = Σ(y_i − ŷ)² + 2 × Σ(y_i − ŷ) × (ŷ − y_medio) + Σ(ŷ − y_medio)²

Mas como:

Σ(y_i − ŷ) × (ŷ − y_medio) = 0

chega-se à equação:

Σ(y_i − y_medio)² = Σ(y_i − ŷ)² + Σ(ŷ − y_medio)²

Esses somatórios são chamados de “somas quadráticas” (SQ), e por isso a equação pode ser representada como:

SQ em torno da média = SQ residual + SQ devido à regressão

ou

SQ_T = SQ_r + SQ_R

Ou seja, parte da variação total das observações (y_i) em torno da média (y_medio) é devida à regressão ( SQ_R = Σ(ŷ − y_medio)²) e parte é devida aos resíduos (SQ_r = Σ(y_i − ŷ)² )

E portanto quanto maior for o percentual de contribuição devido à regressão na soma quadrática total melhor será o “ajuste do modelo”, cujo índice é chamado de “coeficiente de determinação” (R²):

R² = SQ_R / SQ_T

O valor máximo de R² é 1, quando toda a variação em torno da média for devida ao modelo de regressão (ou explicada pelo modelo de regressão), ou seja, SQ_r = Σ(y_i − ŷ)² = 0.

O número de graus de liberdade da soma quadrática total (ν_T) deve ser igual à soma dos graus de liberdade da regressão (ν_R) e do resíduo (ν_r):

ν_T = ν_r + ν_R

O número de graus de liberdade da regressão ν_R é igual ao número de parâmetros p (β₀, β₁ ...) menos 1:

ν_R = p − 1

E o número de graus de liberdade da soma quadrática residual (ν_r) é a diferença entre o número de observações (n) e o número de parâmetros (p):

ν_r = n − p

E como o número de graus de liberdade da soma quadrática total (ν_T) é n − 1 temos a relação geral:

ν_T = ν_r + ν_R

(n − 1) = (n − p) + (p − 1)

No caso de uma regressão com dois parâmetros (β₀ e β₁) teríamos:

(n − 1) = (n − 2) + 1

Esses parâmetros são importantes para montar a “Tabela de Análise de Variância” (ANOVA - Analysis of Variance), que fornece informações sobre os níveis de variabilidade em um modelo de regressão e formam a base para “testes de significância”.

Para montar a tabela de ANOVA dividimos as Somas Quadráticas SQ_R (Σ(ŷ − y_medio)²) e SQ_r (Σ(y_i − ŷ)²) pelos respectivos Graus de Liberdade (ν_R = p − 1) e (ν_r = n − p) e obtemos as respectivas Médias Quadráticas:

Admitindo que os erros seguem uma distribuição normal podemos usar a tabela de ANOVA para fazer uma análise da “Significância Estatística da Regressão”.

Quando β₁ = 0, isto é, quando não há relação entre x e y, pode-se demonstrar que a razão entre as médias quadráticas MQ_R e MQ_r segue uma distribuição “F”:

MQ_R / MQ_r ≈ F_{p-1, n-p}

Onde p − 1 é o número de graus de liberdade da Média Quadrática devida à Regressão (ν_R - MQ_R) e n − p é o número de graus de liberdade da Média Quadrática Residual (ν_r).

Para avaliar a significância da regressão é feito o Teste de Hipóteses:

H₀ - Hipótese nula: β₀ = β₁ = β_i = 0

H₁ - Hipótese alternativa: β₀, β₁ ... ≠ 0

Pela “Hipótese nula” NÃO existe relação linear estatísticamente significativa entre x e y.

E pela “Hipótese alternativa” existe uma relação linear estatísticamente significativa entre x e y.

O teste consiste basicamente em “rejeitar” ou “não rejeitar” a “Hipótese nula”.

A análise é feita calculando-se a razão MQ_R/MQ_r e comparando o resultado com o valor de F tabelado para p − 1 e n − p graus de liberdade, no nível de confiança desejado. (Não confundir que na tabela de distribuição F a confiança de 95% é indicada pela percentagem de 5%)

Se a razão MQ_R/MQ_r for maior do que F_{p − 1, n − p}

MQ_R/MQ_r > F _{p − 1, n − p}

deve-se descartar a possibilidade de que β₀ = β₁ = β_i = 0, ou seja, a hipótese nula é “rejeitada”.

Ou seja, teremos evidência de que existe uma relação linear estatísticamente significativa entre x e y.

Quanto maior a razão MQ_R/MQ_r em relação a F_{p − 1, n − p}, mais evidente se torna a relação.

A Média Quadrática Residual MQ_r pode, também, ser interpretada como uma medida aproximada do erro médio quadrático (s²) associado ao uso da equação de regressão para prever a resposta y_i correspondente a um dado valor x_i.

Quando X_a = X_médio, o erro padrão da estimativa ŷ assume seu valor mínimo. E à medida que nos afastamos desse ponto, em qualquer direção, o erro vai aumentando. Quanto mais longe estivermos de X_a = X_médio, mais incertas serão as previsões feitas a partir da regressão.

E portanto os limites do intervalo de confiança da previsão ŷ variam ao longo da reta de regressão formando “hipérboles”.

>>> SQ_R = np.sum( ( Y_est - np.mean(Y_exp) )**2 )
>>> SQ_R
3465.5999999999985

>>> SQ_r = np.sum( (Y_exp - Y_est)**2 )
>>> SQ_r
832.4000000000004

>>> SQ_T = np.sum( (Y_exp - np.mean(Y_exp) )**2 )
>>> SQ_T
4298.0

>>> GL_R = b.size - 1
>>> GL_R
1

>>> GL_r = X.shape[0] - b.size
>>> GL_r
7

>>> MQ_R = SQ_R / GL_R
>>> MQ_R
3465.5999999999985

>>> MQ_r = SQ_r / GL_r
>>> MQ_r
118.91428571428578

>>> F_calc = MQ_R / MQ_r
>>> F_calc
29.14368092263332

>>> F_tab = f.ppf(q=0.95, dfn=GL_R, dfd=GL_r)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'f' is not defined

>>> from scipy.stats import f
>>> F_tab = f.ppf(q=0.95, dfn=GL_R, dfd=GL_r)
>>> F_tab
5.591447851220738

>>> if F_calc > (10*F_tab):
       print ("F_calc > 10 * F_tab")
       print ("Deve-se rejeitar a hipótese nula (b0 = b1 =  bi = 0).")
       print ("Ou seja, temos evidência de que existe uma relação linear estatisticamente significativa entre x e y. (b0, b1... != 0)")
    else:
       print ("Não se rejeita a hipótese nula (b0 = b1 =  bi = 0).")
       print ("Ou seja, NÃO temos evidência de que existe uma relação linear estatisticamente significativa entre x e y.")

Não se rejeita a hipótese nula (b0 = b1 =  bi = 0).
Ou seja, NÃO temos evidência de que existe uma relação linear estatisticamente significativa entre x e y.
>>>

Se o verdadeiro valor médio de y_i é β₀ + β₁x_i, é de se esperar que observações repetidas no mesmo ponto x_i se distribuam simetricamente em torno de β₀ + β₁x_i, com desvios positivos sendo tão freqüentes quanto desvios negativos, de tal maneira que a média dos erros ε_i seja zero.

Num dado x_i os erros em y_i se distribuirão com uma certa variância σ_i², que em princípio também varia com x_i.

Mas para se calcular o intervalo de confiança das estimativas dos parâmetros (β₀, β₁ ...) serão consideradas as seguintes hipóteses: (Fonte: Como fazer experimentos, 2001)

Estas três hipóteses sobre o comportamento dos erros aleatórios podem ser resumidas nas expressões:

ε_i ≊ N(0, σ²) e Cov(ε_i,ε_j) = 0

Com base nessas suposições pode-se dar prosseguimento ao cálculo do “intervalo de confiança” para os parâmetros do modelo.

Os intervalos de confiança dos parâmetros podem ser calculados pela equação S.5

Admitindo que o valor de s², a variância residual em torno da regressão corresponde à Média Quadrática Residual (MQ_r calculada na tabela ANOVA), seja uma boa estimativa de σ², podemos obter uma estimativa do erro padrão de b_i.

Como também se admite que os erros se distribuem normalmente, pode-se usar a distribuição de Student para testar a significância do valor estimado para b_i.

O número de graus de liberdade do valor de t é n − p, que é o número de graus de liberdade da estimativa s² (MQ_s, e conseqüentemente também do erro padrão.

O uso de matrizes facilita sobremaneira o cálculo dos “erros padrão” dos parâmetros. E isso é feito pelo cálculo da “Matriz de Covariância” dos parâmetro como mostra a figura S.7

Figura S.7. Matriz de covariância dos parâmetros (b₀, b₁ ...) do modelo de regressão.

E essa matriz pode ser calculada pela equação S.6

Essa é outra equação muito importante, que se aplica ao ajuste por mínimos quadrados de qualquer modelo “linear nos parâmetros”, dentro das suposições que foram feitas no início desta seção.

Multiplicando a matriz (X^tX)⁻¹ pela Média Quadrática Residual (MQ_r), como estimativa de σ², e tirando a raiz quadrada dos elementos da diagonal principal, obtemos os erros padrão dos parâmetros (b₀, b₁ ...).

>>> matriz_covar = inv_X_t_dot_X * MQ_r
>>> matriz_covar
array([[ 2.11403175e+02, -3.96380952e+00],
       [-3.96380952e+00,  7.92761905e-02]])

>>> V_B = np.diag(matriz_covar)
>>> V_B
array([2.11403175e+02, 7.92761905e-02])

>>> erro_padrao_B = np.sqrt(V_B)
>>> erro_padrao_B
array([14.53971027,  0.28156028])

>>> from scipy.stats import t
>>> t_crit = t.ppf(0.95, GL_r)
>>> t_crit
1.894578605061305

i = 0

for par_b, erro_b in zip(B, erro_padrao_B):
    print("Parâmetro: %f e Erro: %f" % (par_b.item(), erro_b))
    intervalo = t_crit * erro_b
    print ("Intervalo: ", intervalo)
    b_max = par_b.item() + intervalo
    b_min = par_b.item() - intervalo
    print("Intervalo de confiança do parâmetro b%d (%f), limite superior (%f) e limite inferior (%f)" % (i, par_b.item(), b_max, b_min))
    print ("Isto significa que há 95%s de probabilidade  de que o verdadeiro valor do parâmetro b%d esteja entre %f e %f" % ('%', i, b_min, b_max))
    if abs(intervalo) > abs(par_b.item()):
        print ("Como esses dois limites têm sinais contrários, e como nenhum valor num intervalo de confiança é mais provável do que outro, pode ser que o verdadeiro valor de b%d  seja zero. Ou seja, o valor de b%d = %f não é estatisticamente significativo e portanto não existe evidência suficiente para manter o termo b%d no modelo." % (i, i, par_b.item(), i))
    print()
    i += 1

Parâmetro: -7.333333 e Erro: 14.539710
Intervalo:  27.546623992708057
Intervalo de confiança do parâmetro b0 (-7.333333), limite superior (20.213291) e limite inferior (-34.879957)
Isto significa que há 95% de probabilidade  de que o verdadeiro valor do parâmetro b0 esteja entre -34.879957 e 20.213291
Como esses dois limites têm sinais contrários, e como nenhum valor num intervalo de confiança é mais provável do que outro, pode ser que o verdadeiro valor de b0  seja zero. Ou seja, o valor de b0 = -7.333333 não é estatisticamente significativo e portanto não existe evidência suficiente para manter o termo b0 no modelo.

Parâmetro: 1.520000 e Erro: 0.281560
Intervalo:  0.5334380798399798
Intervalo de confiança do parâmetro b1 (1.520000), limite superior (2.053438) e limite inferior (0.986562)
Isto significa que há 95% de probabilidade  de que o verdadeiro valor do parâmetro b1 esteja entre 0.986562 e 2.053438

Como pode ser visto no gráfico da figura S.2 não e possivel traçar uma reta que passe por todos os pontos experimentais.

Por isso vamos modelar a influência da temperatura sobre o rendimento com um modelo quadrático como mostra a equação S.7

O ajuste deste novo modelo aos valores experimentais também é feito por meio da equação matricial S.2. Mas agora as matrizes precisam ser expandidas para se ajustarem à equação quadrática S.7 como mostra a figura S.8.

Figura S.8. Equação matricial representando o sistema de 9 equações do modelo empírico “quadrático” que descreve a relação rendimento e temperatura descrita na tabela S.1 (Fonte: Como fazer experimentos, 2001)

Neste caso a matriz “X” pode ser montada incluindo a coluna correspondente a “t²” com os comandos:

>>> X = np.array([[30, 35, 40, 45, 50, 55, 60, 65, 70]])
>>> X
array([[30, 35, 40, 45, 50, 55, 60, 65, 70]])
>>> X = X.T
>>> X
array([[30],
       [35],
       [40],
       [45],
       [50],
       [55],
       [60],
       [65],
       [70]])
>>> X = np.hstack([np.ones(X.shape), X, X**2])
>>> X
array([[1.000e+00, 3.000e+01, 9.000e+02],
       [1.000e+00, 3.500e+01, 1.225e+03],
       [1.000e+00, 4.000e+01, 1.600e+03],
       [1.000e+00, 4.500e+01, 2.025e+03],
       [1.000e+00, 5.000e+01, 2.500e+03],
       [1.000e+00, 5.500e+01, 3.025e+03],
       [1.000e+00, 6.000e+01, 3.600e+03],
       [1.000e+00, 6.500e+01, 4.225e+03],
       [1.000e+00, 7.000e+01, 4.900e+03]])
>>> 

Aplicando a equação S.2, conforme os comandos descritos anteriormente, obtemos a matriz “b” com os parâmetros da equação de regressão:

>>> Y_exp = np.array([[24], [40], [60], [70], [77], [86], [91], [86], [84]])
>>> X_t = X.T
>>> X_t_prod_matricial_X = np.dot(X_t,X)
>>> inv_X_t_dot_X = np.linalg.inv(X_t_prod_matricial_X)
>>> X_t_Y_exp = np.dot(X_t,Y_exp)
>>> b = np.dot(inv_X_t_dot_X, X_t_Y_exp)
>>> b
array([[-1.58242424e+02],
       [ 7.98753247e+00],
       [-6.46753247e-02]])

O que significa que o modelo quadrático estima os rendimentos por meio da equação S.8.

A figura S.9 mostra o gráfico da equação S.8, onde se observa que o ajuste é muito melhor do que o modelo linear (Figura S.2).

Figura S.9. Gráfico com os resultados experimentais da tabela S.1 e a curva de regressão do modelo quadrático.(Gráfico feito com o programa Calc)

>>> Y_est = np.dot(X, b)
>>> SQ_R = np.sum( ( Y_est - np.mean(Y_exp) )**2 )
>>> SQ_T = np.sum( (Y_exp - np.mean(Y_exp) )**2 )
>>> SQ_r = np.sum( (Y_exp - Y_est)**2 )
>>> GL_R = b.size - 1
>>> GL_r = X.shape[0] - b.size
>>> MQ_R = SQ_R / GL_R
>>> MQ_r = SQ_r / GL_r
>>> F_calc = MQ_R / MQ_r
>>> from scipy.stats import f
>>> F_tab = f.ppf(q=0.95, dfn=GL_R, dfd=GL_r)
>>> if F_calc > (10*F_tab):
...     print ("F_calc > 10 * F_tab")
...     print ("Deve-se rejeitar a hipótese nula (b0 = b1 =  bi = 0).")
...     print ("Ou seja, temos evidência de que existe uma relação linear estatisticamente significativa entre x e y. (b0, b1... != 0)")
... else:
...     print ("Não se rejeita a hipótese nula (b0 = b1 =  bi = 0).")
...     print ("Ou seja, NÃO temos evidência de que existe uma relação linear estatisticamente significativa entre x e y.")
... 
F_calc > 10 * F_tab
Deve-se rejeitar a hipótese nula (b0 = b1 =  bi = 0).
Ou seja, temos evidência de que existe uma relação linear estatisticamente significativa entre x e y. (b0, b1... != 0)

#Matriz de Covariância (matriz_covar):
>>> matriz_covar = inv_X_t_dot_X * MQ_r
#A diagonal da Matriz de Covariância corresponde às variâncias das estimativas dos parâmetros. (V(b0), V(b1 e V(b2) ...)
>>> V_b = np.diag(matriz_covar)
#Erros padrão dos parâmetros:
>>> erro_padrao_b = np.sqrt(V_b)
#O valor t crítico (t_crit) é calculado para o número de Graus de Liberdade da Soma Quadrática residual (GL_r) e no nível de 95% de confiança:
>>> from scipy.stats import t
>>> t_crit = t.ppf(0.95, GL_r)

>>> i = 0
>>> for par_b, erro_b in zip(b, erro_padrao_b):
...     print("Parâmetro: %f e Erro: %f" % (par_b.item(), erro_b))
...     intervalo = t_crit * erro_b
...     print ("Intervalo: ", intervalo)
...     b_max = par_b.item() + intervalo
...     b_min = par_b.item() - intervalo
...     print("Intervalo de confiança do parâmetro b%d (%f), limite superior (%f) e limite inferior (%f)" % (i, par_b.item(), b_max, b_min))
...     print ("Isto significa que há 95%s de probabilidade  de que o verdadeiro valor do parâmetro b%d esteja entre %f e %f" % ('%', i, b_min, b_max))
...     if abs(intervalo) > abs(par_b.item()):
...         print ("Como esses dois limites têm sinais contrários, e como nenhum valor num intervalo de confiança é mais provável do que outro, pode ser que o verdadeiro valor de b%d  seja zero. Ou seja, o valor de b%d = %f não é estatisticamente significativo e portanto não existe evidência suficiente para manter o termo b%d no modelo." % (i, i, par_b.item(), i))
...     print()
...     i += 1
... 

Parâmetro: -158.242424 e Erro: 11.672005
Intervalo:  22.680810531882145
Intervalo de confiança do parâmetro b0 (-158.242424), limite superior (-135.561614) e limite inferior (-180.923235)
Isto significa que há 95% de probabilidade  de que o verdadeiro valor do parâmetro b0 esteja entre -180.923235 e -135.561614

Parâmetro: 7.987532 e Erro: 0.488315
Intervalo:  0.9488844520661023
Intervalo de confiança do parâmetro b1 (7.987532), limite superior (8.936417) e limite inferior (7.038648)
Isto significa que há 95% de probabilidade  de que o verdadeiro valor do parâmetro b1 esteja entre 7.038648 e 8.936417

Parâmetro: -0.064675 e Erro: 0.004852
Intervalo:  0.009428537605566101
Intervalo de confiança do parâmetro b2 (-0.064675), limite superior (-0.055247) e limite inferior (-0.074104)
Isto significa que há 95% de probabilidade  de que o verdadeiro valor do parâmetro b2 esteja entre -0.074104 e -0.055247

>>> 

Matplotlib é uma biblioteca Python especializada no desenvolvimento de gráficos bi e tridimensionais.

A primeira etapa é carregar a biblioteca:

import matplotlib.pyplot as plt

E executar os comandos:

>>> Y_exp = np.array([[24], [40], [60], [70], [77], [86], [91], [86], [84]])
>>> X = np.array([[30, 35, 40, 45, 50, 55, 60, 65, 70]])
>>> X = X.T
>>> X = np.hstack([np.ones(X.shape), X, X**2])
>>> X_t = X.T
>>> X_t_prod_matricial_X = np.dot(X_t,X)
>>> inv_X_t_dot_X = np.linalg.inv(X_t_prod_matricial_X)
>>> X_t_Y_exp = np.dot(X_t,Y_exp)
>>> b = np.dot(inv_X_t_dot_X, X_t_Y_exp)
>>> Y_est = np.dot(X,b)

>>> X_col_2 = X[:,1]

>>> import matplotlib.pyplot as plt

>>> plt.scatter(X_col_2, Y_exp, c='blue')
>matplotlib.collections.PathCollection object at 0x7f677420b160>
>>> plt.plot(X_col_2, Y_est, c='r')
[>matplotlib.lines.Line2D object at 0x7f678ad62ef0>]
>>> plt.xlabel('Temperatura (C)')
>matplotlib.text.Text object at 0x7f67750a33c8>
>>> plt.ylabel('Rendimento (%)')
>matplotlib.text.Text object at 0x7f67750bc1d0>
>>> plt.title('Gráfico com dados experimentais do rendimento em função da temperatura ')
>matplotlib.text.Text object at 0x7f6774260ac8>
>>> plt.show()

Vamos lembrar que os valores de temperatura (variável independente) foram organizados, tanto no modelo linear quanto no modelo quadrático, na segunda coluna da matriz (array) X. Por isso a segunda coluna foi extraída (slice) com o comando X_col_2 = X[:,1] e armazenada na variável X_col_2.

Com os comandos acima foi possível obter o gráfico dos dados experimentais e da regressao do modelo quadrático como mostra a figura S.10.

Figura S.10. Gráfico gerado com a biblioteca Matplotlib a partir dos dados da tabela S.1 e os dados gerados pelo modelo de regressão quadrático.

Mas também pode ser usado o método plot para exibir pontos ao invés de linhas.

plt.plot(X_col_2, Y_exp, 's', c='blue')

O terceiro argumento “s” define que seja usado o símbolo quadrado (square) para representar os pontos experimentais.

Também podemos criar o gráfico dos “resíduos” (y_exp - y_est) com os comandos:

#Ajuste do tamanho
plt.figure(figsize=(10,5))

#Cria subplot com separação horizontal
plt.subplot(1,2,1)

plt.plot(X_col_2, Y_exp, 's', c='red')

plt.plot(X_col_2, Y_est, c='blue')

plt.xlabel('Temperatura (C)')

plt.ylabel('Rendimento (%)')

plt.title('Dados experimentais')

plt.subplot(1,2,2)

plt.axis([20, X_col_2[-1]+10, -10, 10])

plt.plot(X_col_2, E, 's', c='blue')

plt.hlines(y=0, xmin=20, xmax=X_col_2[-1]+10, color='red')

#Aumentar espaço entre os gráficos
plt.subplots_adjust(wspace=0.4)

plt.show()

Os dois gráficos são exibidos na mesma janela como mostra a figura: S.11.

Figura S.11. No lado esquerdo o gráfico gerado a partir dos dados da tabela S.1 e os dados gerados pelo modelo de regressão quadrático. E no lado direito os dados dos resíduos.

Também podemos adicionar um terceiro gráfico mostrando a “correlação” entre os resultados observados e simulados como mostra a figura S.12.

Figura S.12. No lado esquerdo o gráfico gerado a partir dos dados da tabela S.1 e os dados gerados pelo modelo de regressão quadrático. No centro os dados dos resíduos. E à direita o gráfico com os dados observados e simulados

Para exibir o terceiro gráfico com os resultados simulados e observados usamos os comandos adicionais:

plt.subplot(1,3,3)

plt.plot(Y_exp, Y_est, 's', c='red')

#Dica
#https://stackoverflow.com/questions/40516661/adding-line-to-scatter-plot-using-pythons-matplotlib


x_lim = plt.xlim()
y_lim = plt.ylim()

plt.plot(x_lim, y_lim, c='blue')

plt.xlabel('Experimental')

plt.ylabel('Simulado')

#Aumentar espaço entre os gráficos
plt.subplots_adjust(wspace=0.4)

Para usar a biblioteca é necessário primeiro “carregar” com o comando:

>>> import pandas as pd

Criamos um arquivo arq_Temp_Rend.csv contendo os valores de Temperatura e Rendimento:

T,R
30,24
35,40
40,60
45,70
50,77
55,86
60,91
65,86
70,84

Para abrir o arquivo usamos o comando read_csv:

data = pd.read_csv('arq_Temp_Rend.csv', sep=',')

Mas data é uma variável do tipo “DataFrame”:

>>> type(data)
<class 'pandas.core.frame.DataFrame'>

Por isso fiz primeiro a conversão para o tipo “array” para realizar as operações matriciais com NumPy.

X = data.values[:,0]

X = X.reshape(-1,1)

Y_exp = data.values[:,1]

Y_exp = Y_exp.reshape(-1,1)

Programa exemplo com todos os códigos usados nessa seção.

Controlling an organic synthesis robot with machine learning to search for new reactivity

Notes On Using Data Science & Machine Learning