Matlab u Octave. El método de mínimos cuadrados resuelve directamente el problema del ajuste polinomial. Sin embargo, la solución implica calcular la inversa de una matriz, la cual aumentará de tamaño a medida que el grado del polinomio de ajuste sea mayor. En una próxima entrada se verá método iterativo para aproximar el ajuste polinomial.
Fundamento Teórico
Dado un conjunto de $N$ datos (puntos naranjas en la fig. 1), el objetivo es encontrar un polinomio, $p(x)$ de grado $n$, tal que su gráfica (linea azul en la fig. 1) se ajuste mejor a los datos.\begin{eqnarray}
p(x) &=& a_{0} + a_{1}x + a_{2}x^{2} + \dots + a_{n}^{n}\\
p(x) &=& \sum_{i=0}^{n}a_{i}x^{i}
\end{eqnarray}
Definiremos el error en el i-ésimo dato, $e_{i}$, como
\begin{equation}
e_{i} = q_{i} - y_{i}
\end{equation}
donde $y_{i} = p(x_{i})$.
Definimos el error total, $E$, como
\begin{equation}
E = \dfrac{1}{N}\sum_{i=1}^{N}e_{i}^{2}
\end{equation}
Derivaremos $E$ respecto a cada coeficiente $a_{i}$ e igualaremos a cero esta derivada para encontrar el mínimo, o sea, el mejor ajuste. Es decir, solucionaremos el problema de ajuste mediante
Utilizando la definición de $E$
Utilizando la definición $e_{i} = q_{i} - y_{i}$ obtenemos $\frac{\partial}{\partial a_{k}}e_{i} = -\frac{\partial}{\partial a_{k}}y_{i}$, por lo tanto, reemplazando en la eq. anterior.
Utilizando la definición $y_{i} = p(x_{i}) = \sum_{j=0}^{n}a_{j}x_{i}^{j}$ obtenemos el resultado $\frac{\partial}{\partial a_{k}}y_{i} = x_{i}^{k}$, reemplazando estos hechos en la eq. anterior
Dado que debemos igualar la derivada a cero
Obtenemos finalmente
La eq. anterior puede ser escrita de forma matricial como
Utilizando las variables $S$, $A$ y $Q$ para representar a los factores de la eq. anterior.
Por lo tanto, el valor de los coeficientes $a_{i}$ y por tanto la solución al ajuste, viene dado por
Resultados Experimentales
La siguiente función en 'ajpol_mincuad.m' devuelve los coeficientes del polinomio de ajuste a un conjunto de datos.% ****************************************************************************** % * FUNC: % * ajpol_mincuad % * % * DESC: % * Ajuste de los 'N' datos en 'P' (de dimension Nx2) mediante un % * polinomio de grado 'n'. % * % * ARGS: % * P - <Nx2 double> % * n - double % * % * RETS: % * A - <nx1 double> % ****************************************************************************** function A = ajpol_mincuad(P,n) N = length(P); % 'P' es el conjunto de N puntos S = zeros(n+1); % 'S' es la matriz cuadrada de la forma SUM(x^(j+k)) Q = zeros(n+1,1); % 'Q' es el vector columna de la forma SUM(q*x^k) x = P(:,1); % Vector de componentes 'x' de los datos q = P(:,2); % Vector de componentes 'q' de los datos for k = 0:n for j = k:n S(k+1,j+1) = sum(x.^(j+k)); S(j+1,k+1) = S(k+1,j+1); % 'S' es simetrica end Q(k+1) = sum(q.*x.^k); end A = inv(S)*Q; end
Utilizaremos los siguientes datos en '
data.txt' para probar esta función:1.000000 97.62227 2.000000 97.80724 3.000000 96.62247 4.000000 92.59022 5.000000 91.23869 6.000000 95.32704 7.000000 90.35040 8.000000 89.46235 9.000000 91.72520 10.000000 89.86916 11.00000 86.88076 12.00000 85.94360 13.00000 87.60686 14.00000 86.25839 15.00000 80.74976 16.00000 83.03551 17.00000 88.25837 18.00000 82.01316 19.00000 82.74098 20.00000 83.30034 21.00000 81.27850 22.00000 81.85506 23.00000 80.75195 24.00000 80.09573 25.00000 81.07633 26.00000 78.81542 27.00000 78.38596 28.00000 79.93386 29.00000 79.48474 30.00000 79.95942 31.00000 76.10691 32.00000 78.39830 33.00000 81.43060 34.00000 82.48867 35.00000 81.65462 36.00000 80.84323 37.00000 88.68663 38.00000 84.74438 39.00000 86.83934 40.00000 85.97739 41.00000 91.28509 42.00000 97.22411 43.00000 93.51733 44.00000 94.10159 45.00000 101.91760 46.00000 98.43134 47.00000 110.4214 48.00000 107.6628 49.00000 111.7288 50.00000 116.5115 51.00000 120.7609 52.00000 123.9553 53.00000 124.2437 54.00000 130.7996 55.00000 133.2960 56.00000 130.7788 57.00000 132.0565 58.00000 138.6584 59.00000 142.9252 60.00000 142.7215 61.00000 144.1249 62.00000 147.4377 63.00000 148.2647 64.00000 152.0519 65.00000 147.3863 66.00000 149.2074 67.00000 148.9537 68.00000 144.5876 69.00000 148.1226 70.00000 148.0144 71.00000 143.8893 72.00000 140.9088 73.00000 143.4434 74.00000 139.3938 75.00000 135.9878 76.00000 136.3927 77.00000 126.7262 78.00000 124.4487 79.00000 122.8647 80.00000 113.8557 81.00000 113.7037 82.00000 106.8407 83.00000 107.0034 84.00000 102.46290 85.00000 96.09296 86.00000 94.57555 87.00000 86.98824 88.00000 84.90154 89.00000 81.18023 90.00000 76.40117 91.00000 67.09200 92.00000 72.67155 93.00000 68.10848 94.00000 67.99088 95.00000 63.34094 96.00000 60.55253 97.00000 56.18687 98.00000 53.64482 99.00000 53.70307 100.00000 48.07893 101.00000 42.21258 102.00000 45.65181 103.00000 41.69728 104.00000 41.24946 105.00000 39.21349 106.0000 37.71696 107.0000 36.68395 108.0000 37.30393 109.0000 37.43277 110.0000 37.45012 111.0000 32.64648 112.0000 31.84347 113.0000 31.39951 114.0000 26.68912 115.0000 32.25323 116.0000 27.61008 117.0000 33.58649 118.0000 28.10714 119.0000 30.26428 120.0000 28.01648 121.0000 29.11021 122.0000 23.02099 123.0000 25.65091 124.0000 28.50295 125.0000 25.23701 126.0000 26.13828 127.0000 33.53260 128.0000 29.25195 129.0000 27.09847 130.0000 26.52999 131.0000 25.52401 132.0000 26.69218 133.0000 24.55269 134.0000 27.71763 135.0000 25.20297 136.0000 25.61483 137.0000 25.06893 138.0000 27.63930 139.0000 24.94851 140.0000 25.86806 141.0000 22.48183 142.0000 26.90045 143.0000 25.39919 144.0000 17.90614 145.0000 23.76039 146.0000 25.89689 147.0000 27.64231 148.0000 22.86101 149.0000 26.47003 150.0000 23.72888 151.0000 27.54334 152.0000 30.52683 153.0000 28.07261 154.0000 34.92815 155.0000 28.29194 156.0000 34.19161 157.0000 35.41207 158.0000 37.09336 159.0000 40.98330 160.0000 39.53923 161.0000 47.80123 162.0000 47.46305 163.0000 51.04166 164.0000 54.58065 165.0000 57.53001 166.0000 61.42089 167.0000 62.79032 168.0000 68.51455 169.0000 70.23053 170.0000 74.42776 171.0000 76.59911 172.0000 81.62053 173.0000 83.42208 174.0000 79.17451 175.0000 88.56985 176.0000 85.66525 177.0000 86.55502 178.0000 90.65907 179.0000 84.27290 180.0000 85.72220 181.0000 83.10702 182.0000 82.16884 183.0000 80.42568 184.0000 78.15692 185.0000 79.79691 186.0000 77.84378 187.0000 74.50327 188.0000 71.57289 189.0000 65.88031 190.0000 65.01385 191.0000 60.19582 192.0000 59.66726 193.0000 52.95478 194.0000 53.87792 195.0000 44.91274 196.0000 41.09909 197.0000 41.68018 198.0000 34.53379 199.0000 34.86419 200.0000 33.14787 201.0000 29.58864 202.0000 27.29462 203.0000 21.91439 204.0000 19.08159 205.0000 24.90290 206.0000 19.82341 207.0000 16.75551 208.0000 18.24558 209.0000 17.23549 210.0000 16.34934 211.0000 13.71285 212.0000 14.75676 213.0000 13.97169 214.0000 12.42867 215.0000 14.35519 216.0000 7.703309 217.0000 10.234410 218.0000 11.78315 219.0000 13.87768 220.0000 4.535700 221.0000 10.059280 222.0000 8.424824 223.0000 10.533120 224.0000 9.602255 225.0000 7.877514 226.0000 6.258121 227.0000 8.899865 228.0000 7.877754 229.0000 12.51191 230.0000 10.66205 231.0000 6.035400 232.0000 6.790655 233.0000 8.783535 234.0000 4.600288 235.0000 8.400915 236.0000 7.216561 237.0000 10.017410 238.0000 7.331278 239.0000 6.527863 240.0000 2.842001 241.0000 10.325070 242.0000 4.790995 243.0000 8.377101 244.0000 6.264445 245.0000 2.706213 246.0000 8.362329 247.0000 8.983658 248.0000 3.362571 249.0000 1.182746 250.0000 4.875359
El siguiente script en '
mincuad.m' pone a prueba la funcion 'ajpol_mincuad' con los datos anteriores y con grado 20 del polinomio de ajuste.clear; clc; close all; P = load('data.txt'); % Carga los datos en 'P' n = 20; % Grado del polinomio de ajuste A = ajpol_mincuad(P,n); A = A(n+1:-1:1); % Invierte el orden de los coeficientes xval = min(P(:,1)):0.1:max(P(:,1)); yval = polyval(A,xval); plot(xval,yval,'-r'); % Plot del polinomio de ajuste hold on; plot(P(:,1),P(:,2),'ob'); % Plot de los datos
La gráfica en la salida es la siguiente:
Favor si me ayuda con la rutina de ejecutar un programa para elaboracion de presupuestos de obra con componentes de costos unitarios de los rubros del presupuesto a n68mendoza@gmail.com
ResponderEliminar