3.3.3 Aproximação de Funções — Regressões — Regressão Polinomial

1. Introdução

Nos dois últimos posts apresentamos o desenvolvimento para ajustar dados em uma equação da reta através de mínimos quadrados. Contudo, existem muitos eventos em que o modelo obedece à um comportamento polinomial. Para tais modelos é necessário adaptar o ajuste para uma função polinomial de grau superior. A imagem abaixo ilustra um caso em que a regressão linear não ajusta adequadamente os dados ao comportamento do evento observado.

 

2. Regressão Polinomial

A regressão polinomial pode ser vista como uma generalização da regressão linear. Para isso, podemos ver a regressão linear simples como o ajuste de um polinômio de grau um. Assim, ao invés de ajustarmos a função


\(y = \alpha_0 + \alpha_1 x + \epsilon \)

utilizamos

\(y = \alpha_0 + \alpha_1 x + \alpha_2 x^2 + \cdots + \alpha_m x ^m + \epsilon \)

Para ajustarmos os parâmetros dessa função, basta que resolvamos um sistema de \(m + 1 \) equações lineares simultâneas, similarmente ao desenvolvimento da regressão linear múltipla. Assim, no caso da regressão polinomial, o erro padrão pode ser formulado como

\(s = \sqrt{\dfrac{S_r}{n – (m + 1)}}\)
.

Como exemplo, vamos tomar o ajuste dos dados apresentados no gráfico do início desse texto. Desta figura, podemos notar que o comportamento dos pontos segue uma função parabólica, o que sugere que devemos ajustar um polinômio de segundo grau:


\(y = \alpha_0 + \alpha_1 x + \alpha_2 x^2 + \epsilon \)

Nesse caso, a soma dos quadrados dos resíduos é

\(S_r = \sum_{i=1}^{n} \left( y_i – \alpha_0 – \alpha_1 x_i – \alpha_2 x_{i}^2 \right)^2 \)

Assim como fizemos para regressão linear, tomamos a derivada da equação acima com relação à cada um dos coeficientes desconhecidos do polinômio:


\(\dfrac{dS_r}{d \alpha_0} = -2 \sum_{i=1}^{n} \left( y_i – \alpha_0 – \alpha_1 x_i – \alpha_2 x_i^2 \right) \)
\(\dfrac{dS_r}{d \alpha_1} = -2 \sum_{i=1}^{n} x_i \left( y_i – \alpha_0 – \alpha_1 x_i – \alpha_2 x_i^2 \right) \)
\(\dfrac{dS_r}{d \alpha_2} = -2 \sum_{i=1}^{n} x_i^2 \left( y_i – \alpha_0 – \alpha_1 x_i – \alpha_2 x_i^2 \right) \)

Igualando estas equações à zero e obtemos o seguinte conjunto de equações:

\(n \alpha_0 + \alpha_1 \sum_{i=1}^{n} x_i + \alpha_2 \sum_{i=1}^{n} x_i^2 = \sum_{i=1}^{n} y_i \)
\(\alpha_0 \sum_{i=1}^{n} x_i + \alpha_1 \sum_{i=1}^{n} x_i^2 + \alpha_2 \sum_{i=1}^{n} x_i^3 = \sum_{i=1}^{n} y_i x_i \)
\(\alpha_0 \sum_{i=1}^{n} x_i^2 + \alpha_1 \sum_{i=1}^{n} x_i^3 + \alpha_2 \sum_{i=1}^{n} x_i^4 = \sum_{i=1}^{n} y_i x_i^2 \)

Dessas três últimas equações lineares, temos apenas três incógnitas: \(\alpha_0 \) , \(\alpha_1 \) e \(\alpha_2 \) . Assim, resolvendo essas equações, temos os parâmetros ajustados.

 

3. Implementação

Como podemos notar, para um conjunto de \(n = m + 1 \) pontos amostrados, é possível ajustar um polinômio de grau \(m \) ou menor. Na implementação abaixo utilizamos um ajuste automático do grau polinomial, escolhendo o valor de \(m \) que minimize o erro padrão.

Note que na função abaixo, utilizamos o método de Cholesky para resolução do sistema linear. Poderíamos utilizar qualquer outro método, não sendo esse escolhido obrigatório para o ajuste correto dos dados.

def polinomial_regression(points, order = -1):
    """
    Return a FUNCTION with parameters adjusted by polinomial regression
 
    f = polinomial_regression(points)
 
    INPUT:
      * points: list of tuples of sampled points (x_i, y_i)
      * order: order of adjusted polynom (-1 to auto-adjust)
 
    return: a polinomial function f(x) = a_0 + a_1 * x + ... + a_m * x^m
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Sep 2011
    """
    def prepare_linear_system(x, y, order):
        A = [[0.0 for i in xrange(order)] for j in xrange(order)]
        B = [0.0 for i in xrange(order)]
        for row in xrange(order):
            for col in xrange(row, order):
                c = col + row
                xc = [i ** c for i in x]
                s = sum(xc)
                A[row][col] = s
                A[col][row] = s
        for col in xrange(order):
            bc = [(xi ** col) * yi for (xi, yi) in zip(x, y)]
            B[col] = sum(bc)
        return (A, B)
 
    def select_degree(points, order):
        if order < 1:
            return len(points)
        else:
            return order + 1
 
    n = select_degree(points, order)
    (x, y) = zip(*points)
    (A, b) = prepare_linear_system(x, y, n)
    alphas = cholesky(A, b)
    fun = lambda t: sum([t * a for a in alphas])
    return fun

Esta função está disponível em http://www.sawp.com.br/code/regression/polinomial_regression.py assim como um exemplo de como utilizá-la.

 

4. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

 

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001).

[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall (2006).