Enviar um café pro programador

Regressão Linear em Python - O método dos mínimos quadrados

 Neste tutorial, vamos aprender como aplicar o método dos mínimos quadrados para fazer uma regressão linear em seus dados experimentais.

O que é Regressão Linear

Como o nome pode sugerir, é uma técnica onde tentaremos transformar uma série de dados em uma    linha.

Ela é muito útil para tentarmos estimarmos ou prevermos novos valores, sem precisarmos realizar novamente o experimento.

Por exemplo, vamos supor que os pontos azuis abaixo sejam referente aos experimentos que você realizou em laboratório:

Método dos mínimos quadrados

No exemplo que fizemos, usamos uma mola e várias massas. Para cada massa m, temos um tamanho de mola l. Colocamos a massa no eixo X e o comprimento no eixo Y.

Note na imagem acima, que a maioria dos pontos ficarem mais ou menos alinhados. O que a regressão linear faz é tentar achar uma reta 'ótima', mais próxima desses pontos.

Assim, você pode estimar qual seria o valor de uma mola para uma determinada massa, sem precisar fazer o experimento. Obviamente, é uma aproximação.

Método dos mínimos quadrados

Existem vários métodos e maneiras de se fazer regressão, inclusive de qualquer tipo de curva que você imaginar, não somente de retas.

O método mais conhecido, é o do mínimos quadrados. Veja que existe uma distância de cada bolinha (experimento que você fez) até nossa reta otimizada por regressão linear:

Teoria dos erros método dos mínimos quadrados

O nosso objetivo é pegar essa distância, que  é o erro e minimizar o máximo possível eles, e gerar a reta de equação: 

  • y = A + B.x

Você pode achar a demonstração desse método em livros universitários que falam sobre Teoria dos Erros. As fórmulas para A e B são:

Análise de erros

Dados coletados e fórmula

No laboratório de Física da Universidade Federal do Ceará, pegamos uma mola menor e uma maior (apropriadas para as experiências), colocamos elas pendurada na vertical, e calculamos o tamanho final dela, para cada massa. Fizemos 7 experimentos em cada mola, e o resultado foi:
Experimento de Hooke da mola

Para nosso tutorial, vamos gerar o gráfico para a mola menor. Queremos provar que a lei de Hooke para o sistema massa-mola é obedecido.

Equilíbrio do sistema massa-mola:

  • −k∆x + mg = 0

Nos experimentos, medimos o tamanho l final da mola, que é obtido pela expressão acima:

  • L = Lo + (k/g)*m

Vamos fazer a regressão linear dessa fórmula, através do método dos mínimos quadrados, com a reta:

  • y = A + Bx
Onde:

y ⇔ l

x ⇔ m

Regressão Linear: Fazendo os Cálculos e gráficos

Obviamente, existe tudo isso (e muito mais) prontinho em Python, para ser usado. Bastando só inserir seus dados e ele vai fazer tudo. Mas, provavelmente seu professor vai pedir que você faça isso 'na mão', pelo menos a primeira vez.

E isso é o correto. Implementar, na unha, as coisas vai te fazer um cientista ou engenheiro bem mais capacitado, vai te fazer entender melhor como as coisas funcionam e caso precise mudar algo, vai poder mexer no código e na matemática.

É muito importante, principalmente para TCC, mestrado, doutorado etc, entender bem o que está ocorrendo, pois você pode ter que fazer algo novo, para seu trabalho.

Então vamos, lá.

Primeiro, importamos a biblioteca numpy (para os cálculos) e a matplotlib.pyplot (para gerar o gráfico). 

Depois, criamos dois arrays, um com as massas (massa,em gramas) e outro com o valor das distensões da mola (delta_x, em centímetros), medidos no laboratório. 

Seja L0 o comprimento inicial da mola, então os comprimentos finais serão dados por L0 + delta_x

Como fizemos 7 medições, N=7. E agora é só fazer os cálculos das fórmulas de A e B.
Fizemos cada uma separadamente, para melhor entender, e depois jogamos tudo para ter as variáveis A e B. Por fim, geramos o array y com os novos dos comprimentos estimados, parar gerar a reta.

Colocamos, para aqueles que desejarem, o cálculo das incertezas de A, B, y e da constante da mola. Lembrando que o parâmetro A fica em centímetros e B fica com cm/gramas.

Como a gravidade usada é 9.8 m/s², não se esqueça de passar as outras unidades para o SI.

Nossos dados experimentais, os pontos, foram plotados com a scatter (dispersão), depois usamos uma reta pontilhada para plotar a linha da regressão, e por fim um gráfico de erros.


Note que, ao final do script, coloquei vários prints, para ir exibindo os valores das variáveis e cálculos realizados, uma etapa importante para você ir acompanhando e conferindo se está tudo ok. É muito fácil e comum cometermos erros ao longo do desenvolvimento de scripts para assuntos científicos.

Nosso script fica:

import numpy as np
import matplotlib.pyplot as plt

# Massas(em gramas) e respectivas Distensões da mola deltax (em centímetros)
massa = np.array([10.7, 20.6, 30.8, 41.1, 50.9, 60.8, 71.2])
delta_x = np.array([6.4, 8.0, 8.8, 10.0, 11.0, 12.8, 13.5])

#Comprimentos totais da mola
L0 = 5.0
L = L0 + delta_x

# Seja: l = l0 + Bm a equação da reta que relaciona o tamanho
# final da mola (l), com seu tamanho inicial (l0) e o valor da massa (m)
# B representa a constante g/k (onde g é a acelaração da gravidade) e k
# a constante da mola

## Regressão linear: y = A + Bx
# x é a massa e y é o l

# Número de experimentos
N = 7

# Somatório de x[i] (massas)
S_m = massa.sum()

# Somatório de l[i] (comprimentos)
S_l = L.sum()

# Somatório dos quadrados x²[i] das massas
S_m2 = sum(i*i for i in massa)

# Somatório dos quadrados dos comprimentos l²[i]
S_l2 = sum(i*i for i in L)

# Somatório dos produtos de cada massa pelo seu respectivo comprimento final m[i]*l[i]
S_ml = 0.0
for i in range(N):
  S_ml = S_ml + massa[i]*L[i]

# Colocando tudo na fórmula
delta = N*S_m2 - (S_m*S_m)
A = ((S_m2 * S_l) - S_m * S_ml)/(delta)
B = ((N*S_ml) - S_m * S_l)/(delta)

# Reta com os novos valores, Yi
y = np.array(A + B*massa)

# Incerteza de y, A e B, com uma unidade de incerteza
residuo=0.0
y=y.round(1)
A=round(A,1)
B=round(B,1)

for i in range(N):
  residuo = residuo + np.power(y[i] - A - B*massa[i], 2)

dy = round(np.sqrt(residuo/(N-2)),1)
dA = round(dy * np.sqrt(S_m2/delta),1)
dB = round(dy * np.sqrt(N/delta),1)

# B e seu erro dB estão em centimetros/gramas = 10 metros/kilogramas
B_kgm = B*10.0
dB_kgm = dB*10.0

# Constante da mola, e seu erro propagado de B
# pois B = g/k
K = 9.8/B_kgm
dK = dB_kgm * ( 9.8/(B_kgm*B_kgm) )

# Gerando o gráfico
fig, ax = plt.subplots()
fig.suptitle('Massa(g) x Tamanho da mola (cm)')
plt.rcParams.update({'font.size': 12})

ax.title.set_text('Curso Python Progressivo')
ax.scatter(massa, L, marker='.', s=1)
ax.plot(massa, y, color='red', linestyle='--', label='reta mmq')
ax.errorbar(massa, L, yerr=dy, fmt='.k', label='dados experimentais')
ax.legend(loc="lower right")
plt.show()


print("Somatorio das massas: ", S_m)
print("Somatorio quadrado das massas: ", S_m2)
print("Somatorio comprimento: ", S_l)
print("Somatorio quadrado do scomprimentos: ", S_l2)
print("Somatorio massa pelo comprimento menores: ", S_ml)
print("Delta: ", delta)
print("A: %.2f cm"% A)
print("B: %.2f cm/g"% B)
print("Reta da mola: y = %.1f + %.1f*m" % (A, B))
print("Erro y: ", dy)
print("Erro A: ", dA)
print("Erro B: ", dB)
print("Constante: ", K)
print("Erro de k: ", dK)
print("K = %.1f +- %.1f" %(K, dK))

Gráfico da regressão linear:

Método dos mínimos quadrados em Python


Código no github:

https://github.com/jarlissonmoreira/PythonProgressivo/blob/main/Regressao%20Linear/RegressaoLinear.py


Regressão Linear usando a biblioteca SciPy

Agora vamos deixar o Python fazer (quase) tudo por nós. Para isso, importamos a biblioteca SciPy, e desde o conjunto de funções stats, onde iremos uar a função linregress para fazer a regressão linear de nossos dados.

Lembre-se: você vai precisar de dois conjuntos de valores de mesmo tamanho, para passar para a linregress. Essa função retorna um conjunto de parâmetros, como o A e B (intercept e slope), previamente calculados, além de outros, como erros padrões.

Nosso código fica:

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def func(x):
  return result.slope * x + result.intercept

# Massas(g) e respectivas Distensões da mola (cm)
massa = np.array([10.7, 20.6, 30.8, 41.1, 50.9, 60.8, 71.2])
delta_x = np.array([6.4, 8.0, 8.8, 10.0, 11.0, 12.8, 13.5])

# Comprimento total da mola
L0=5.0
L = np.array(L0 + delta_x)

# Armazenando os parametros
# Retorna a tupla com: B(coefieciente angular da reta), intercept (A, onde toca o eixo x), 
# rvalue (coeficiente Pearson de correlação), pvalue (valor-p da probabilidade), 
# stderr (erro padrão do parâmetro slope), intercept_stderr (erro padrão do parâmetro intercept)
result = stats.linregress(massa, L)

modelo = list(map(func, massa))

# Gerando o gráfico
fig, ax = plt.subplots()
fig.suptitle('Massa(g) x Tamanho da mola (cm), com SciPy')
plt.rcParams.update({'font.size': 12})
ax.title.set_text('Curso Python Progressivo')
ax.scatter(massa, L, marker='.')
ax.plot(massa, modelo, color='red', linestyle='--', label='reta mmq')
plt.show()

## Imprimindo os dados
print("B (coeficiente angular): ",result.slope)
print("A (onde incercepta o eixo x): ", result.intercept)
print("rvalue : ", result.rvalue)
print("pvalue : ", result.pvalue)
print("stderr (erro padrão do slope, B): ", result.stderr)
print("intercept_stderr (erro padrão do incercept, A): ",result.intercept_stderr)

E o gráfico:

SciPy, stats e linregress

Ao final, imprimimos os valores calculados pela linregress(). Veja que os valores de A e B realmente bateram com aqueles que calculamos antes. Porém, não deu exato, pois em nossos cálculos tivemos o cuidado de arredondar para uma casa decimal, pois nossas medições tinham apenas uma casa decimal. Sempre tomem cuidado com os algoritmos significativos.

Nenhum comentário:

Postar um comentário