Dados pares de observações $(x_j, y_j), 1 \leq j \leq m, x_j \in \mathbb{R}^n, m > n$, queremos encontrar uma função $g(x; \theta), \theta \in \mathbb{R}^k$ tal que $g(x_j; \theta) = y_j$ para uma escolha adequada dos parâmetros $\theta$. A função $g(x; \theta)$ é o modelo utilizado e precisa ser definido previamente.
Idealmente, o modelo $g(x; \theta)$ seria flexível o suficiente para representar a transformação desejada. Na prática, porém, o modelo $g(x; \theta)$ pode somente aproximar o valor correto $y_j$ a partir de $x_j$. Desta forma, podemos definir o erro de aproximação do par $(x_j, y_j)$ como
\[ r_j(\theta) = g(x_j; \theta) - y_j \]
e minimizar a função objetivo
\[ f(\theta) = \frac{1}{2}\sum^m_j [r_j(\theta)]^2 = \frac{1}{2}\sum^m_j [ g(x_j; \theta) - y_j ]^2. \]
Ao minimizar $f(\theta)$ encontramos o vetor de parâmetros $\theta$ que "melhor" ajusta o modelo $g(x; \theta)$ aos dados $(x_j, y_j)$. Neste caso, "melhor" é definido como o vetor $\theta$ que minimiza a soma do quadrado dos erros. Por isto o nome de método de Mínimos Quadrados.
Justificativa estatística
Suponha que os $r_j(\theta)$ são independentes e identicamente distribuídos com média 0, variância $\sigma^2$ e função densidade de probabilidade $g_\sigma$. Sob estas condições, a probabilidade de observar os $r_j(\theta)$ correspondentes ao conjunto de dados $(x_j, y_j)$, dadas a variância $\sigma$ e o vetor de parâmetros $\theta$ é:\[p( (x, y) ; \theta ; \sigma) = \prod_j^m g_\sigma( r_j(\theta) ) = \prod_j^m g_\sigma( g(x_j; \theta) - y_i ) \]
Fixadas as observações $(x_j, y_j)$ e $\sigma$, o valor "mais provável" de $\theta$ é o que maximiza $p( (x, y) ; \theta ; \sigma)$. Este estimador para $\theta$ é chamado de Estimador de Máxima Verossimilhança.
Suponha, então, que a distribuição dos $r_j$ é normal. Portanto,
\[ g_\sigma(r) = \frac{1}{\sqrt{2\pi \sigma^2}} \mathrm{exp}(-\frac{r^2}{2\sigma^2}), \]
e então,
\[p( (x, y) ; \theta ; \sigma) = \frac{1}{(2\pi\sigma^2)^{m / 2}} \mathrm{exp}(-\frac{1}{2\sigma^2} \sum_j^m [g(x_j; \theta) - y_i )]^2 ). \]
Para maximizar $p( (x, y) ; \theta ; \sigma)$ (dado $\sigma^2$ fixo) é necessário minimizar $\sum_j^m [g(x_j; \theta) - y_i )]^2$. Portanto, a estimativa de $\theta$ obtida pelo método de Mínimos Quadrados equivale ao estimador de máxima verossimilhança de $\theta$ (supondo que os erros $r_j$ são i.i.d. com distribuição normal de média 0 e variância fixa $\sigma^2$).
Modelos Lineares
Diversos problemas podem ser modelados utilizando funções $g(x; \theta)$ lineares em $\theta$. Podemos, portanto, escrever os resíduos como $r_j(\theta) = x_j^T\theta - y_j$.Seja $J \in \mathbb{R}^{m\times n}$ tal que
\[ J =
\begin{pmatrix}
x_1^T \\
\dots \\
x_m^T\\
\end{pmatrix}, \]
então o problema de mínimos quadrados se reduz a minimizar a função
\[ f(\theta) = \frac{1}{2} || J\theta - y ||_2^2 \]
Modelos lineares são vantajosos pois a função acima pode ser minimizada globalmente. Isto não é necessariamente verdade para modelos não lineares em $\theta$.
Note que a função $g(x; \theta)$ deve ser linear em $\theta$, mas não necessariamente em $x$. Logo, a função
\[ g(x; \theta) = {x^{(1)}}^2 \theta_{(1)} + x^{(1)}x^{(2)}\theta_{(2)} \]
é um modelo válido para a estimação linear, pois ela equivale a construir um novo conjunto de dados $(x'_j, y_j)$, onde $x'_j = ({x^{(1)}_j}^2, x_j^{(1)}x_j^{(2)})$, e fazer a minimização linear "comum".
Exemplos usando Python
A biblioteca Scipy possui tanto a implementação do Método de Mínimos Quadrados para o caso linear quanto para o caso geral. Neste texto mostrarei como utilizar a função numpy.linalg.lstsq para estimar modelos lineares.Iremos estimar duas funções diferentes. Ambas podem ser estimadas usando modelos lineares. A função $f_1$ é um polinômio de segundo grau e a $f_2$ é uma função de duas variáveis.
\[ f_1(x) = 2x^2 + 5x + 7 \]
\[ f_2(x_1, x_2) = 0.3x^3 - 0.2y^2 + 3y - xy + 5 \]
Estimação de polinômios
Utilizaremos três modelos para estimar $f_1$:
\[ g_1(x) = p_1x + p_2 \]
\[ g_2(x) = p_1x + p_2x^2 + p_3 \]
\[ g_3(x) = p_1x + p_2x^2 + p_3x^3 + p_4 \]
É esperado que o modelo $g_1$ possua erro grande, pois um polinômio de grau 1 não consegue aproximar um polinômio de grau 2. Por outro lado, tanto $g_2$ quanto $g_3$ deverão obter boas aproximações, sendo que em $g_3$ o componente $x^3$ deverá ser próximo de 0.
Primeiramente, iremos gerar nosso conjunto de observações de $f_1$. Em cada observação também adicionaremos um pequeno erro.
x_data = [1, 4, 7, 8, -2, 5, 6, 12, -3, -5, -10]
y_data = [2 * x*x + 5*x + 7 - (random.random() * 5 - 2.5) for x in x_data]
Agora, iremos definir a matrix $J$ para cada um dos modelos e aplicar a função lstsq para estimar os parâmetros de cada modelo.
m = len(x_data)
x2 = np.array([x*x for x in x_data]).reshape((m, 1))
x3 = np.array([x*x*x for x in x_data]).reshape((m, 1))
x_data = np.reshape(x_data, (m, 1))
y_data = np.reshape(y_data, (m, 1))
ones = np.ones(m).reshape(m, 1)
J1 = np.hstack([x_data, ones])
J2 = np.hstack([x2, x_data, ones])
J3 = np.hstack([x3, x2, x_data, ones])
theta_g_1 = np.linalg.lstsq(J1, y_data )
theta_g_2 = np.linalg.lstsq(J2, y_data )
theta_g_3 = np.linalg.lstsq(J3, y_data )
Os resultados obtidos podem ser vistos nas figuras abaixo:
Como esperado, a estimação linear foi bem ruim e a quadrática deu resultados quase perfeitos. O polinômio de grau 3 também teve resultados bons e apresentou um coeficiente bem pequeno para o termo $x^3$.
Estimação (linear) de funções de várias variáveis
\[ g_1(x, y) = \theta_1 x + \theta_2 y + \theta_3 xy + \theta_4\]
\[ g_2(x, y) = \theta_1 x^2 + \theta_2 x y + \theta_3 y^2 + \theta_4 x + \theta_5 y + \theta_6\]
\[ g_3(x, y) = \theta_1 x^3 + \theta_2 x^2 y + \theta_3 x y^2 + \theta_4 y^3 + \theta_5 x^2 + \theta_6y^2 + \theta_7xy + \theta_8x + \theta_9y + \theta_{10}\]
Primeiramente, geramos as observações usando o código abaixo. Note que é necessário um número muito maior de pontos para estimar funções mais complexas.
data = np.array([(random.randint(-50, 50), random.randint(-50, 50)) for i in range(100)])
y_data = [0.3*x**3 - 0.2*y**2 + 3*y - x*y + 5 - (random.random() * 5 - 2.5) for (x, y) in data]
Após, geramos os vetores com as componentes de cada modelo e as matrizes $J_1, J2$ e $J_3$.
m = data.shape[0]
x3 = np.array([x**3 for (x, y) in data]).reshape((m, 1))
x2 = np.array([x**2 for (x, y) in data]).reshape((m, 1))
x1 = np.array([x for (x, y) in data]).reshape((m, 1))
y3 = np.array([y**3 for (x, y) in data]).reshape((m, 1))
y2 = np.array([y**2 for (x, y) in data]).reshape((m, 1))
y1 = np.array([y for (x, y) in data]).reshape((m, 1))
x2y = np.array([x**2 * y for (x, y) in data]).reshape((m, 1))
y2x = np.array([x * y**2 for (x, y) in data]).reshape((m, 1))
xy = np.array([x*y for (x, y) in data]).reshape((m, 1))
ones= np.ones(m).reshape((m, 1))
J_1 = np.hstack([x1, y1, xy, ones])
J_2 = np.hstack([x2, xy, y2, x1, y1, ones])
J_3 = np.hstack([x3, x2y, y2x, y3, x2, y2, xy, x1, y1, ones])
theta_g_1, err1, _, _ = np.linalg.lstsq(J_1, y_data)
theta_g_2, err2, _, _ = np.linalg.lstsq(J_2, y_data)
theta_g_3, err3, _, _ = np.linalg.lstsq(J_3, y_data)
Os resultados obtidos foram os seguintes:
Erro: 3.73790848e+09
Erro: 3.64574253e+09
Erro: 204.59062094
Uma desvantagem do método de Mínimos Quadrados é que a minimização feita dificilmente zera algum $\theta_i$, pois ele pode super ajustar a curva utilizando pequenos coeficientes em componentes que podem ser irrelevantes. Desta maneira, o método não é capaz de detectar quais parâmetros/características são relevantes para um dado problema.
Comentários finais
A estimação utilizando o método de Mínimos Quadrados otimiza os parâmetros de um modelo de modo a minimizar a soma dos erros ao quadrado. Este método encontra mínimos globais se o modelo utilizado for linear (nos parâmetros a serem estimados) e é vantajoso quando se conhece ou se tem alguma intuição sobre a forma analítica do modelo original dos dados.Atualização: veja a parte 2 deste texto, onde apresento um método de solução direta do problema e uma formulação regularizada chamada Ridge Regression que evita overfitting.
