12

Preciso implementar um algoritmo de regressão linear. De preferencia que dê resultados iguais ou próximos a da função TENDENCIA (ou TREND) do Excel.

Já encontrei muito material que tenta explicar todo o conceito de regressão linear, mas não é isso o que eu quero. Gostaria de um pseudo código que seja diretamente traduzível para uma linguagem de programação imperativa.

Para quem não usou a função TENDENCIA ela funciona assim:

Dados dois vetores que representam x e y, em um par ordenado (x,y), para dos valores conhecidos, por exemplo, [ (1,10), (2,20), (3,30), (4,40) ], e um valor de x para o qual você deseja descobrir o valor de y correspondente, a função TENDENCIA encaixa os valores conhecidos em uma função da forma y=mx+b e te da o valor de y. Neste caso se eu passar o valor 5 como último argumento, a função me retorna 50.

Alguém é capaz de dar um pseudo código de um algoritmo que faça isso?

1
  • 1
    Cara, fiz isso na faculdade há mais de 10 anos, mas não consigo me recordar (nem achar referências). Também estou à procura de algo assim. Nesse blog post há uma implementação, bem simples, em Python - mas segundo o autor ela é mais "didática" que "eficiente" (afinal, existem zilhões de maneiras de se fazer regressão linear, cada uma com seus prós e contras).
    – mgibsonbr
    Commented 5/02/2014 às 17:08

7 Respostas 7

9

A resposta do @Marcos Banik pode ser traduzida por exemplo para Python com um mínimo esforço:

>>> x = [1,2,3,4]
>>> y = [10,20,30,40]
>>> m = sum(a*b for (a,b) in zip(x,y)) - sum(x) * sum(y) / len(x)
>>> m /= sum(a**2 for a in x) - (sum(x)**2)/len(x)
>>> b = sum(y)/len(y) - m * sum(x)/len(x)
>>> m*5 + b
50

Entretanto nem toda linguagem imperativa possui o mesmo nível de expressividade, por isso vou apresentar uma versão mais elaborada (em JavaScript):

var x = [1,2,3,4];
var y = [10,20,30,40];

function produto(x, y) {
    var ret = [];
    for ( var i = 0 ; i < x.length ; i++ )
        ret.push(x[i] * y[i]);
    return ret;
}

function quadrados(x) {
    var ret = [];
    for ( var i = 0 ; i < x.length ; i++ )
        ret.push(x[i] * x[i]);
    return ret;
}

function somatorio(x) {
    var ret = 0;
    for ( var i = 0 ; i < x.length ; i++ )
        ret += x[i];
    return ret;
}

function media(x) { return somatorio(x) / x.length; }

var m = somatorio(produto(x,y)) - somatorio(x) * somatorio(y) / x.length;
m /= somatorio(quadrados(x)) - somatorio(x)*somatorio(x) / x.length;

var b = media(y) - m * media(x);

console.log(m*5 + b);

Exemplo no jsFiddle.

8

Montei um pequeno exemplo em Java baseado em alguns exemplos que achei na Web.

Utilizei os dados que você forneceu...

public class Main {

    public static void main(String[] args) {
        double[] x = {1, 2, 3, 4};
        double[] y = {10, 20, 30, 40};
        System.out.println(trend(y, x, 5));
    }

    public static double trend(double[] known_y, double[] known_x, double new_x)
    {
        double[] values = LeastSquaresFitLinear(known_y, known_x);
        return (values[0] * new_x) + values[1];
    }

    public static double[] LeastSquaresFitLinear(double[] known_y, double[] known_x)
    {
        double M, B;
        if (known_y.length != known_x.length)
        {
            return new double[]{0,0};
        }

        int numPoints = known_y.length;

        double x1, y1, xy, x2, J;

        x1 = y1 = xy = x2 = 0.0;
        for (int i = 0; i < numPoints; i++)
        {
            x1 = x1 + known_x[i];
            y1 = y1 + known_y[i];
            xy = xy + known_x[i] * known_y[i];
            x2 = x2 + known_x[i] * known_x[i];
        }

        M = B = 0;
        J = ((double)numPoints * x2) - (x1 * x1);

        if (J != 0.0)
        {
            M = (((double)numPoints * xy) - (x1 * y1)) / J;
            B = ((y1 * x2) - (x1 * xy)) / J;
        }
        return new double[]{M,B};
    }

}
1
  • 1
    Perfeito! A solução em Java ajuda já que é quase uma lingua franca e o algoritmo funciona, não só para os valores apresentados como para outros valores, de maneira muito próxima, senão identica, a função TENDENCIA/TREND do Excel. Commented 5/02/2014 às 18:06
5

Eu não sei qual é o algorítimo que a função TENDENCIA utiliza, mas se for apenas uma variável de entrada x e uma de saída y, a resposta para m e b pelo métodos dos mínimos quadrados é: (código em R)

m <- (sum(x*y) - sum(x) * sum(y) / n) / ( sum(x^2) - ( sum(x) )^2 ) / n )
b <- mean(y) - m * mean(x)

n é o número de elementos dos vetores x e y, sum é a função que retorna a soma dos vetores, x * y é o vetor com o produto dos elementos de mesmo índice. x^2 é o vetor onde todos os seus elementos estão elevados ao quadrado

você pode encontrar essa solução na wikipedia

6
  • O que é j? E esse x*y significa o que em R? (produto escalar - acho que não - produto cruzado, etc. Em J, x*y seria um vetor com os produtos dos elementos individuais de x e y, mas não creio que seja a mesma coisa em R)
    – mgibsonbr
    Commented 5/02/2014 às 17:11
  • j foi um erro de edição (deveria ser y). x*y significa a mesma coisa em J e R pelo o que entendi. Commented 5/02/2014 às 17:17
  • 1
    Sim.Na verdade em R é melhor utilizar a função lm para calcular regressões, já que ela retorna o intervalo de confiança, residuais e outros parâmetros importantes na análise da qualidade da regressão. Commented 5/02/2014 às 17:24
  • 1. Eu pergunto explicitamente "em uma linguagem imperativa" pois sei que apareceriam respostas em linguagem funcionais já que problemas assim são do domínio delas. Mesmo assim obrigado pela informação. 2. A ordem de precedência das operações são iguais as matemáticas? A substração será a última operação realizada? Não consegui reproduzir o seu exemplo no Excel assumindo essa precedência. Commented 5/02/2014 às 17:29
  • 1
    @RaulAlmeida De minha parte, achei a resposta facilmente "traduzível" para uma linguagem imperativa (entretanto, minha linguagem preferida é Python - que engloba várias características funcionais no seu design - além de ter experiência prévia em J), mas entendo que pra nem todo mundo é assim. Postei uma resposta separada onde expresso esse mesmo método num formato puramente imperativo.
    – mgibsonbr
    Commented 5/02/2014 às 17:54
4

Utilizando o mesmo algoritmo da resposta do @Marcos Banik, escrevi um pequeno algoritmo em C baseado no método dos mínimos quadrados:

void lms(double *x, double *y, int n, double *m, double *b)
{
    int i;
    double sumYX = 0.;
    double sumX = 0.;
    double sumY = 0.;
    double sumX2 = 0.;
    double sum2X = 0.;
    for(i = 0; i < n; i++) {
        sumYX += x[i] * y[i];
        sumX += x[i];
        sumY += y[i];
        sumX2 += x[i] * x[i];
    }
    sum2X = sumX * sumX;

    *m = (sumYX - (sumX * sumY) / (double)n) / (sumX2 - sum2X / (double)n);
    *b = sumY / (double)n - *m * sumX / (double)n;
}

O algoritmo é exatamente o mesmo utilizado na resposta dele. O laço for inicial é feito para calcular todos as somatórias e posteriormente seus resultados são utilizados para achar m e b. Acredito que isso seja facilmente traduzível para sua linguagem de preferência. :)

Com m e b você tem a equação da reta. Para calcular qual o valor de y, basta fazer uma função:

double trend(double m, double b, double x)
{
    return m*x + b;
}

Para chamar as funções:

int main(void) 
{
    double x[] = {1., 2., 3., 4.};
    double y[] = {10., 20., 30., 40.};
    double m, b;
    double nx = 5.;
    double ny;

    lms(x, y, 4, &m, &b);
    ny = trend(m, b, nx);

    printf("m: %lf \nb: %lf \nnx: %lf \nny: %lf\n", m, b, nx, ny);

    return 0;
}
2

Para quem possa interessar aqui está a implementação que eu estava procurando.

Ela está escrita na "linguagem" Clipper (compilador x/Harbour) e funciona igual a TENDENCIA para diversos valores que testei.

Function Main()

   LOCAL x := { 1,2,3,4 }
   LOCAL y := { 10,20,30,40 }

   ? Trend(y, x, 5)

   Return NIL

Function Trend( known_y, known_x, new_x )

   LOCAL v := LSFL(known_y, known_x)

   Return (v[1] * new_x) + v[2]

Function LSFL(known_y, known_x)

   LOCAL M,B
   LOCAL n
   LOCAL x1, y1, xy, x2, J
   LOCAL i

   IF Len( known_y ) != Len( known_x )
      Return { 0, 0 }
   ENDIF

   n := Len( known_y )

   x1 := 0; y1 := 0; xy := 0; x2 := 0

   For i := 1 To n
      x1 := x1 + known_x[i]
      y1 := y1 + known_y[i]
      xy := xy + known_x[i] * known_y[i]
      x2 := x2 + known_x[i] * known_x[i]
   Next

   M := 0
   B := 0

   J := ( n * x2 ) - ( x1 * x1 )

   IF !Empty( J )
      M := (( n * xy ) - ( x1 * y1 )) / J
      B := ( (y1 * x2) - ( x1 * xy )) / J
   ENDIF

   Return { M, B }
2

Você precisa de y=mx+b e para isso basta os dados x. y servira apenas para descobrir outras variaveis, mas nao para definir a equação de regressão, coeficiente de correlação, media, mediana e moda dos dados.

Escopo generico: "qual o ano (x) em que o crescimento populacional absoluto será zero?"

Primeiro precisamos descobrir a taxa de crescimento (%) e se ela é positiva ou negativa, o que define TREND ou a inclinação do gráfico gerado pela ER.

Com base nos dados existentes (anos/x) definimos a coeficiente de correlação ou precisão das respostas para os novos y.

No meio do caminho definimos moda, mediana, média etc para sustentar a solução do problema que permite trabalhar qualquer matriz de x:y dependendo da TREND, o sinal de bx ou mx + ou -.

y = a + bx + c ou y=mx+b

O código completo que gera a equação de regressão e coeficiente de correlação está em The State of Art de James Holmes nos fontes do capitulo 8.

1

Segue um algorítimo que fiz para calculo do MMQ com base numa vídeo-aula que assisti no Youtube:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>


void main(){

    float x[] = {-1,0,1,2};
    float y[] = {-1,1,3,5};
    float sx;
    float xy;
    float totalx = 0;
    float totaly = 0;
    float totalsx = 0;
    float totalxy = 0;
    int n = sizeof(x)/sizeof(float);

    int cont;

    for(cont = 0;cont<n;cont++){
        totalx = totalx + x[cont];
        totaly = totaly + y[cont];
        totalsx = totalsx + pow(x[cont],2);
        totalxy = totalxy + (x[cont]*y[cont]);
    }

    // Passo1
    float vbb = n;
    float vba = totalx;
    float somvba = totaly;
    float b2p1 = totalx;
    float a2p1 = totalsx;
    float soma2p1 = totalxy;

    //Passo 2
    float b1p2 = vbb / vbb;
    float a1p2 = vba / vbb;
    float soma1p2 = somvba / vbb;
    float b2p2 = b2p1-(b1p2*totalx);
    float a2p2 = a2p1-(a1p2*totalx);
    float soma2p2 = soma2p1-(soma1p2*totalx);

    // Passo 3
    float b1p3 = b1p2;
    float a1p3 = a1p2;
    float soma1p3 = soma1p2;
    float a2p3 = a2p2 / a2p2;
    float soma2p3 = soma2p2 / a2p2;

    float afinal = soma2p3 / a2p3;
    float bfinal = soma1p3 - a1p3 * afinal;

    printf("\nValor final de A= %.5f e valor de B= %.5f\n\n", afinal, bfinal);
    system("pause >> log");
}
1
  • 1
    Parece um código interessante, poderia editar a pergunta e explicar como ele funciona?
    – Syzoth
    Commented 13/10/2015 às 0:13

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .