quarta-feira, 11 de junho de 2008

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

Paradigma funcional 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:
y=\sum_{i=0}^{n-1}\Big(y_i\prod_{j=0\\j\neq i}^{n-1}\frac{x-x_j}{x_i-x_j}\Big)

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
blog comments powered by Disqus