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 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: , e . Assim, resolvendo essas equações, temos os parâmetros ajustados.

 

3. Implementação

Como podemos notar, para um conjunto de pontos amostrados, é possível ajustar um polinômio de grau ou menor. Na implementação abaixo utilizamos um ajuste automático do grau polinomial, escolhendo o valor de 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.

 

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).