Cálculos estatísticos: polinómio de Lagrange

Dando sequência aos artigos sobre cálculos estatísticos no Kodumaro – e em outros blogs amigos –, tive a ideia de falar sobre interpolação polinomial.

Por ser um método simples, decidi usar o polinómio de Lagrange:


Não é preciso calcular os coeficientes – é possível simplificar o cálculo final transformando parte do polinómio em um vetor de coeficientes, mas é mais fácil de visualizar sem fazê-lo.

C++

Dessa vez não vou usar C, mas C++, para poder usar recursos de orientação a objetos, como encapsulamento.

Primeiro vamos criar um namespace kodumaro e uma classe Lagrange. Abra o arquivo lagrange.h:

#ifndef LAGRANGE_H_
#define LAGRANGE_H_

namespace kodumaro {
    class Lagrange {
        public:
            Lagrange(const double *, const double *, int);
            ~Lagrange() {}

            double getY(double);

        private:
            const double *vector_x;
            const double *vector_y;
            int           length;
    };
}

#endif /*LAGRANGE_H_*/

O construtor irá receber um vetor de abscissas, um vetor de ordenadas e o tamanho dos vetores. Esses vetores representarão as coordenadas dos pontos a serem interpolados.

O método getY() retornará a ordenada para uma determinada abscissa, de acordo com a interpolação de Lagrange.

Agora o arquivo lagrange.cc contendo a implementação dos métodos:

#include "lagrange.h"

using namespace kodumaro;

Lagrange::Lagrange(const double *vx, const double vy, int len):
    vector_x(vx), vetory(vy), length(len) {}

double Lagrange::getY(double x) {
    double y = 0;

    for (int i = 0; i < length; ++i) { // ciclo do somatório
        double aux = vector_y[i];
        for (int j = 0; j < length; ++j) // ciclo do produtório
            if (i != j)
                aux *= (x - vector_x[j]) / (vector_x[i] - vector_x[j]);
        y += aux;
    }

    return y;
}

E está pronto! Basta incluir lagrange.h e compilar lagrange.cc junto com seu código teste.

O uso é simples:

kodumaro::Lagrange lagr(vx, vy, len);
y = lagr.getY(x);

Common Lisp

Não seria legal fazer isso sem usar uma linguagem funcional. =)

Também vamos aproveitar orientação a objetos – sim! É possível combinar orientação a objetos e programação funcional – para guardar nossos dados.

A classe será similar à de C++, mas o atributo será uma lista de pares de coordenadas:

(defclass lagrange ()
  ((pairs :initarg :pairs :accessor lagrange-pairs)))

Agora o método para calcular a ordenada em relação a uma abscissa:

(defmethod lagrange-get-y ((a-lagrange lagrange) x)
  (reduce '+
    (loop
      for '(xi . yi) in (lagrange-pairs a-lagrange)
      for i from 0
      collect (* yi (reduce '*
        (loop
          for '(xj . yj) in (lagrange-pairs a-lagrange)
          for j from 0
          collect (if (= i j)
            1
            (/ (- x xj) (- xi xj)))))))))

A relação entre esse código e a expressão matemática acima é praticamente 1:1.

O uso é:

(let ((a-lagrange (make-instance 'lagrange :pairs coordinates)))
  (setq y (lagrange-get-y a-lagrange x)))

Smalltalk

Se foi legal usar orientação a objetos numa linguagem onde a orientação foi costurada como em uma colcha de retalhos e em uma linguagem funcional, imaginem em uma linguagem realmente orientada a objetos!

Smalltalk é uma das linguagens mãe da orientação a objetos: criada antes de C a partir de Simula, é praticamente a linguagem que define os conceitos de orientação a objetos.

Como já disse antes, tornei-me fã do Squeak, devido a sua aplicabilidade educacional e ao XO.

No Squeak, abra o System Browser, se não houver um pacote Kodumaro crie-o e, dentro dele, uma categoria Kodumaro-Interpolation.

Na categoria, crie um classe:

PointArray subclass: #Lagrange
    instanceVariableNames = ''
    classVariableNames = ''
    poolDictionaries: ''
    category: 'Kodumaro-Interpolation'

Agora crie o método para calcular a ordenada para uma abscissa:

getY: x
    | y |

    y := 0.
    1 to: self size do: [ :i | | aux |
        aux := (self at: i) y.
        1 to: self size do: [ :j |
            (i ~= j) ifTrue: [
                aux := aux * (
                    (x - (self at: j) x)
                    /
                    ((self at: i) x - (self at: j) x)
                )
            ]
        ].
        y := y + aux.
    ].

    ↑ y

E pronto! Você pode testar instanciando Lagrange – o método new: recebe o tamanho do vetor como parâmetro – e passando a mensagem getY:.

Io

Io é uma linguagem prototipada baseada principalmente em Self e Lua, mas também em outras, como Smalltalk e Lisp.

O código Io pode ser:

Point := Object clone do(
    x := 0
    y := 0
)

Lagrange := List clone do(
    // Gostei da ideia desse método add
    add := method(sx, sy,
        p := Point clone
        p x = sx
        p y = sy
        append(p)
    )

    getY := method(x,
        y := 0
        for(i, 0, size - 1,
            aux := at(i) y
            for(j, 0, size - 1,
                // Havia feito diferente antes,
                // mas acho o if assim mais claro:
                if(i != j) then(
                    aux = aux * ( \
                        (x - at(j) x) \
                        / \
                        (at(i) x - at(j) x) \
                    )
                )
            )
            y = y + aux
        )
        return y
    )
)

Vou dar um exemplo de código pois, devido a ser uma linguagem recente, não há muita documentação.

Digamos que você salvou o código acima no arquivo lagrange.io:

Io> Importer FileImporter import("lagrange")
Io> v := Lagrange clone
Io> v add(1, 10)
Io> v add(2, 100)
Io> v add(3, 1000)
Io> v getY(2)
==> 100
Io> v getY(2.5)
==> 448.75

Conclusão

Bem, desta vez deixo a conclusão ao leitor e convido a quem quiser fazer também suas implementações.

Outro método muito interessante é o polinómio de Newton.

[]'s
Cacilhas, La Batalema