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++
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 | #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:
