Metodos Numericos: Regresión Polinomial

post original: https://metodosnumericosgera.wordpress.com/2017/04/06/regresion-polinomial/

Antecedentes

La aproximación polinomial es describir el comportamiento de una nube de datos a manera de una curva, en lugar de una recta. Se busca que esta recta difiera lo menos posible a la nube de datos.

Podríamos, por ejemplo, intentar aproximar la nube de datos a un polinomio de grado 2:

De tal forma que se minimice la suma de los cuadrados de los residuos, Sr:

Que, de igual manera, podría tratarse de un polinomio de m-ésimo grado:

El error estándar de estimación se calcula igual que en una regresión lineal:

El coeficiente de determinación y correlación también se calculan de la misma manera:

Justificación y propósito

Hemos hablado antes sobre la regresión lineal. Pero, ¿qué sucede cuando la nube de datos no está cerca a comportarse de manera lineal? Podríamos darnos cuenta si esto ocurre si al hacer la aproximación lineal tenemos un error estándar demasiado grande. Es aquí donde sería mejor usar regresión polinomial.

Diagrama de flujo

Codigo en C++

#include <iostream>
#include <fstream>
#include <math.h>
#include <cmath>

using namespace std;
void regresionPolinomial(double x[], double y[], int n){

 const int grado = 2;            //Cambiar el grado para cambiar el grado del polinomio resultante <------- ESTE SI SE CAMBIA DEPENDIENDO DEL PROBLEMA

 double matrix[grado+1][grado+1];
 double solutions[grado+1];

 double sumX = 0;
 double sumXY = 0;

 cout << "------------SISTEMA DE ECUACIONES-------------" << endl;
 //Inicializamos el sistema de ecuaciones
 for (int i = 0; i < grado+1; ++i)
 {

  sumXY = 0;
  for (int k = 0; k < n; ++k)
  {
   sumXY += pow(x[k],i)*y[k];
  }
  solutions[i] = sumXY; //Get the solutiona of the ecuation system

  for (int j = 0; j < grado+1; ++j)
  {
   if (i == 0 && j == 0)
   {
    matrix[i][j] = n;   //Start the matrix[0,0] as n 
   }else{
    sumX = 0;
    for (int k = 0; k < n; ++k)
    {
     sumX += pow(x[k],(j+i));
    }

    matrix[i][j] = sumX;        //Get the other values of the matrix
   }

   cout << matrix[i][j] << ", ";
  }
  cout << "= " << solutions[i] << endl;
 }

 //Ahora resolvemos el sistema de ecuaciones con algun método conocido
 cout << "-------------------GAUSS SEIDEL---------------------" << endl;
 int size = grado +1;
 double EPSILON = 0.0000001; // error iterativo como tolerancia

 double S[size]; //Array with the solutions
 double OS[size];  //Array with old solutions
 double E[size]; //Array with th error of each variable
 for (int i = 0; i < size; ++i)
 {
  S[i] = 0;  //We initialize each solution at 0;
  OS[i] = 0; //Pseudo minus infinity at each old solution;

  E[i] = 0;
 }

 bool errorCheck = false;

 int counter=0;
 
 while(!errorCheck){
  for (int i = 0; i < size; ++i)
  {
   double sum = 0;
   for (int j = 0; j < size; ++j)
   {
    if (j != i)
    {
     sum += matrix[i][j]*S[j];
    }
   }
   OS[i] = S[i];  //The old value is the actual value
   S[i] = (solutions[i] - sum)/matrix[i][i];  //The actual value changes to the substitution of the variable

   E[i] = (S[i]-OS[i]);  //Error iterativo absoluto
   //E[i] = (S[i]-OS[i])/S[i];  // Error iterativo relativo
   //E[i] = ((S[i]-OS[i])/S[i])*100;  // Error iterativo porcentual
  }

  for (int i = 0; i < size; ++i)
  {
   if(abs(E[i]) > EPSILON){
    break;
   }
   errorCheck = true;
  }
  counter++;
 }

 for (int i = 0; i < size; ++i)
 {
  cout  << " a" << i << ": " << S[i] << endl;
 }
 cout << "it took " << counter << " times to terminate" << endl;
 //-------------------Aqui termina GaussSeidel------------------------

 cout << "------------------POLINOMIO---------------" << endl;
 for (int i = size-1; i >= 1; --i)
 {
  cout << S[i] << "X^"<< i << " + ";
 }
 cout << S[0] << endl;

 cout << "------------------ERRORES------------------" << endl;
 double e2[n]; //Calcula el error para cada X
 for (int i = 0; i < n; ++i)
 {
  double Ycalculada = 0;
  for (int j = size-1; j >= 1; --j)
  {
   Ycalculada += S[j]*(pow(x[i], j));
  }
  Ycalculada += S[0];
  e2[i] = pow(y[i]-Ycalculada,2);
 }

 double sr = 0;
 double sumY = solutions[0];
 double st = 0;
 for (int i = 0; i < n; ++i)
 {
  sr += e2[i];
  st += pow(y[i] - (sumY/n),2);
 }
 cout << "Sr: " << sr << endl;

 double sxy = sqrt(sr/(n-size));
 cout << "St: " << st << endl;

 double r2 = (st-sr)/st;
 double r = sqrt(r2);
 cout << "Error Estandar (Sx/y): " << sxy << endl;
 cout << "Coeficiente de determinacion (r2): " << r2 << endl;
 cout << "Coeficiente de correlacion (r): " << r << endl;
 cout << endl;
}

int main(){
 //CON ESTE SI DA TODO BIEN
 double x[] = {0,1,2,3,4,5};
 double y[] = {2.1,7.7,13.6,27.2,40.9,61.1};
 int dataNum = 6;
 
 regresionPolinomial(x,y,dataNum);
}

Ejemplos

Valores:

x	y
0	2.1
1	7.7
2	13.6
3	27.2
4	40.9
5	61.1

Aproximación con polinomio de grado 2:

y = 1.86071x^2 + 2.35929x + 2.47857

Gráfica:

___________________________________________________________________

Valores:

x	y
10.2	69.81
11.4	82.75
11.5	81.75
12.5	80.38
13.1	85.89
13.4	75.32
13.6	69.81
15.0	78.54
15.2	81.29
15.3	99.20
15.6	86.35
16.4	110.23
16.5	106.55
17.0	85.50
17.1	90.02

Aproximación con polinomio de grado 4:

y = -0.0025011X^4 + 0.0836578X^3 + -0.21331X^2 + -11.8438X^1 + 157.18

Gráfica:

Metodos Numericos

martes, 9 de mayo de 2017

Regresión Polinomial