Los polinomios de alto grado no sobreajustan los datos
Se trata de un error común que, sin embargo, se encuentra en muchos libros de texto. En general, para especificar un modelo estadístico, es necesario especificar tanto una clase de hipótesis como un procedimiento de ajuste. Para definir la varianza del modelo ("varianza" aquí en el sentido de Bias-Variance Tradeoff), es necesario saber cómo ajustar el modelo, por lo que afirmar que una clase de hipótesis (como los polinomios) puede sobreajustar los datos es simplemente un error de categoría.
Para ilustrar este punto, considere que el método de interpolación de Lagrange es sólo una de las varias formas de ajustar un polinomio a una colección de puntos. Otro método es utilizar Polinomios de Bernstein . Dado $n$ puntos de observación $x_i=i/n$ y $y_i=f(x_i)+\epsilon_i$ la aproximación de Bernstein es, como la de Lagrange, un polinomio de grado $n-1$ (es decir, tantos grados de libertad como puntos de datos). Sin embargo, podemos ver que los polinomios de Bernstein ajustados no muestran las oscilaciones salvajes que presentan los polinomios de Lagrange. Y, lo que es más sorprendente, la varianza de los ajustes de Bernstein es la siguiente disminuye a medida que añadimos más puntos (y por tanto aumentamos el grado de los polinomios).
Como puedes ver, el polinomio de Bernstein no pasa por cada punto. Esto refleja el hecho de que utiliza un procedimiento de ajuste diferente al del polinomio de Lagrange. Pero ambos modelos tienen exactamente el mismo clase de hipótesis (el conjunto de grado $n-1$ polinomios), lo que subraya que es incoherente hablar de propiedades estadísticas de una clase de hipótesis sin especificar también la función de pérdida y el procedimiento de ajuste.
![enter image description here]()
![enter image description here]()
Como segundo ejemplo, también es posible interpolar exactamente los datos de entrenamiento utilizando un polinomio sin oscilaciones salvajes, asumiendo que estamos trabajando sobre números complejos . Si los puntos de observación $x_i$ se encuentran a lo largo del círculo unitario $x_i=e^{2\pi i\sqrt{-1}/N}$ , y luego ajustar un polinomio $y_i=\sum_d c_dx_i^d$ se reduce básicamente a calcular la transformada inversa de Fourier de $y_i$ . A continuación trazo la interpolante polinómica compleja en función del argumento complejo $i/N$ . Además, resulta que la varianza de este interpolante es constante en función de $n$ .
![enter image description here]()
El código para estas simulaciones se encuentra al final de esta respuesta.
En resumen, Es simplemente falso que la inclusión de polinomios de mayor grado aumente la varianza del modelo, o lo haga más propenso al sobreajuste
En consecuencia, cualquier explicación de por qué los interpolantes de Lagrange hacen de hecho un sobreajuste tendrá que utilizar la estructura real del modelo y el procedimiento de ajuste, y no puede basarse en argumentos generales de recuento de parámetros.
Afortunadamente, no es demasiado difícil demostrar matemáticamente por qué los interpolantes de Lagrange presentan este sobreajuste. Básicamente se trata de que el ajuste de la interpolante de Lagrange requiere la inversión de una matriz de Vandermonde, que a menudo está mal condicionada.
Más formalmente, lo que queremos hacer es demostrar que la varianza de los interpolantes de Lagrange crece muy rápidamente en función de $n$ el número de puntos de datos.
Supongamos que se nos da un número fijo de $x_i$ puntos, así como $y$ valores muestreados de $y_i\sim f(x_i)+\epsilon_i$ donde $f$ es una función suave y $\epsilon$ son muestras iid de una distribución normal unitaria. La interpolación de Lagrange toma entonces la forma $\hat{f}_y(x)=\sum_i l_i(x)y_i$ donde $l_i$ son los polinomios de Lagrange del $x_i$ . Queremos calcular $Var_y(\hat{f}(x))$ donde la varianza es con respecto a diferentes valores muestreados del $y$ s. Por la independencia del $\epsilon_i$ s, tenemos
$$Var_y(\hat{f}(x))=Var_y(\sum_i l_i(x)(f(x_i)+\epsilon_i))=\sum_i l_i(x)^2$$
Ahora, consideremos para simplificar el caso en el que los puntos de observación son $x_1,\dots, x_n=1,\dots, n$ . En general, la varianza $Var(\hat{f}(x))$ es una función de la ubicación $x$ pero consideraremos la varianza integrada $Var=\int_0^nVar(\hat{f}(x))dx$ .
Reclamación: $Var$ aumenta (al menos) como $O(e^{an}n^{-9/2})$ para algunos $a>0$ .
En general, tenemos la fórmula $l_i(x)={\frac {l(x)}{l'(x_i)(x-x_i)}}$ donde $l(x):=\prod_i (x-x_i)$ . Podemos calcular el denominador: $$l'(x_j)=\prod_{i\neq j}(j-i)=(j-1)(j-2)...(1)(-1)...(-(n-j))=(-1)^{n-j}(j-1)! (n-j)!$$
Así que la varianza en $x$ es $l(x)^2 \sum_j {\frac 1 {(j-1)!^2(n-j)!^2(x-x_j)^2}}$ . Ahora, $|x-x_j|<=n/2$ Así que ${\frac 1 {(x-x_j)^2}}\geq 4/n^2$
\begin{eqnarray*} \sum_j {\frac 1 {(j-1)!^2(n-j)!^2}}{\frac 1 {(x-x_j)^2}}&\geq &4n^{-2}n!^{-2}\sum_j {n\choose j}^2\\ & = &4n^{-2} n!^{-2}{2n\choose n}\\ &\sim& n^{-7/2} 4^n(e/n)^{2n} \end{eqnarray*}
donde utilizamos la fórmula de Stirling en la última línea.
Así que para obtener un límite inferior $\int_0^n Var(\hat{f}(x))dx$ necesitamos integrar $l(x)^2$ . Observando que $l(x)^2=n^{2n}\prod_i(x'-i/n)^2$ , donde $x'=x/n$ tenemos
$$l(x)^2=n^{2n}\prod_i (x'-i/n)^2=n^{2n}e^{2\sum_i \log |x'-i/n|}\sim n^{2n}e^{2n\int_0^1 \log |x'-y|dy}=n^{2n}e^{2n(-H(x')-1)}$$ donde $H$ es la función de entropía binaria,
$$\int_0^n l(x)^2x\sim n^{2n-1}\int_0^1 e^{2n(-H(x')-1)} dx'\geq n^{2n-1}e^{-2C(\epsilon)n}\int_{1-\epsilon}^1 dx'=n^{2n-1}e^{-2Cn}/\epsilon$$
donde $C(\epsilon)=1+H(1-\epsilon)$ y $0<\epsilon<1/2$ es arbitraria.
La varianza integrada $\int_0^n Var(\hat{f}(x))dx$ es, por tanto, al menos $O(e^{-2Cn}n^{2n-9/2}4^n(e/n)^{2n})=O(n^{-9/2}({\frac {2e}{e^C}})^{2n})$
Ahora es crucial, $C(\epsilon)<1+\log(2)$ porque la función de entropía binaria se maximiza en $1/2$ . Por lo tanto, $2ee^{-C}>1$ , por lo que obtenemos la forma reclamada para la varianza.
Edición: análisis de la varianza de Bernstein Como comparación interesante, no es demasiado difícil calcular la varianza del aproximador de Bernstein. Este polinomio está definido por $\hat{B}(x)=\sum_i b_{i,n}(x)y_i$ , donde $b_{i,n}$ son los polinomios de la base de Bernstein, como se define en el enlace de Wikipedia. Argumentando de forma similar a la anterior, la varianza con respecto a diferentes muestreos de $y$ es $\sum_i b_{i,n}(x)^2$ . Al igual que en la simulación, promediaré sobre $x\in [0,1]$ . Evidentemente, $\int_0^1 b_{i,n}(x)^2dx={n\choose i}^2B(2i+1,2n-2i+1)$ con $B$ la función beta. Tenemos que sumar esto sobre $i$ . La suma que reclamo asciende a: $$\int_0^1 Var(\hat{B}(x))dx={\frac {4^n}{(2n+1){2n\choose n}}}$$
(la prueba está abajo)
Esta fórmula sí que se ajusta a los resultados empíricos de la simulación (tenga en cuenta que para la comparación directa, tendrá que multiplicar por $.1^2$ que fue la varianza del ruido utilizada en la simulación). Por la fórmula de Stirling, la varianza del aproximador de Bernstein es $O(n^{-.5})$ .
Prueba de la fórmula de la varianza El primer paso es utilizar la identidad $B(2i+1,2n-2i+1)=(2i)!(2n-2i)!/(2n)!$ . Así que obtenemos
\begin{eqnarray*} {n\choose i}^2 B(2i+1,2n-2i+1) & = & {\frac {(n!)^2}{(i!)^2(n-i)!^2}}{\frac {(2i)!(2n-2i)!}{(2n+1)!}}\\ & = & {\frac 1 {2n+1}}{2n\choose n}^{-1}{2i\choose i}{2n-2i\choose n-i} \end{eqnarray*}
Utilizando la función generadora $\sum_{i=0}^{\infty} {2i\choose i}x^i=(1-4x)^{-1/2}$ vemos que $\sum_{i=0}^n {2i\choose i}{2n-2i\choose n-i}$ es el coeficiente de $x^n$ en $((1-4x)^{-1/2})^2=(1-4x)^{-1}$ Es decir $4^n$ . La suma de $i$ obtenemos la fórmula reclamada para la varianza.
Edición: efecto de la distribución de los puntos x en la varianza de Lagrange Como se muestra en la interesantísima respuesta de Robert Mastragostino, el mal comportamiento del polinomio de Lagrange puede evitarse mediante la elección juiciosa de $x_i$ . Esto plantea la posibilidad de que el polinomio de Lagrange se comporte de forma poco habitual para los puntos muestreados uniformemente. En principio, dados los puntos de muestra iid $x_i\sim \mu$ es posible que los polinomios de Lagrange se comporten de forma razonable para la mayoría de las elecciones de $\mu$ y la distribución uniforme es sólo una "elección desafortunada". Sin embargo, este no es el caso, y resulta que la varianza de Lagrange sigue creciendo exponencialmente para cualquier elección de $\mu$ , con la única excepción de la distribución arcoseno $d\mu(x)=\textbf{1}_{(-1,1)}\pi^{-1}(1-x^2)^{-1/2}dx$ . Nótese que en la respuesta de Robert Mastragostino, el $x_i$ los puntos se eligen como nodos de Chebyshev. Como $n\to\infty$ la densidad de estos puntos converge nada menos que a la distribución del arcoseno. Así que en cierto sentido, la situación en su respuesta es la sólo excepción a la regla del crecimiento exponencial de la varianza.
Más precisamente, tome cualquier distribución continua $\mu(dx)=p(x)dx$ apoyado en $(-1,1)$ y considerar una secuencia infinita de puntos $\{(x_i,f(x_i)+\epsilon_i)\}_{i=1,\infty}$ donde $x_i\sim \mu$ son iid, $f$ es una función continua y $\epsilon_i$ son muestras normales iid. Sea $\hat{f_n}(x)$ denotan la interpolante de Lagrange ajustada a la primera $n$ puntos. Y que $V_n=\int Var(\hat{f_n})(x)p(x)dx$ denotan la varianza media correspondiente.
Reclamación: Supongamos que (1): $p(x)>0$ para todos $x\in (-1,1)$ y (2): $\int p(x)^2\sqrt{1-x^2}dx<\infty$ . Entonces $V_n$ crece al menos como $O(e^{\epsilon n})$ para algunos $\epsilon>0$ , a menos que $p$ es la distribución del arcoseno .
Para demostrarlo, obsérvese que, como antes, $Var(\hat{f_n}(x))=\sum_{i=1}^nl_{i,n}^2$ , donde ahora $l_{i,n}$ denota el i-ésimo polinomio de la base de Lagrange con respecto a la primera $n$ puntos. Ahora,
\begin{eqnarray*} \log \int l_{i,n}^2(x)d\mu(x)&\geq& 2\int \log |l_{i,n}(x)|p(x)dx\\ & = & 2\int \sum_{j\neq i, j\leq n} \left(\log |x-x_j|-\log |x_i-x_j|\right)p(x)dx \end{eqnarray*} donde utilizamos la desigualdad de Jensen en la primera línea, y la definición de $l_{i,n}$ en el segundo.
Introduzcamos la siguiente notación: \begin{eqnarray*} f_p(x)&=&\int \log|x-y|p(y)dy\\ c_p& = & \int \log |x-y|p(x)p(y)dxdy \end{eqnarray*}
Entonces obtenemos
\begin{eqnarray*}\lim\inf n^{-1}\log \int l_{i,n}^2(x)p(x)dx&\geq& 2\lim\inf_n n^{-1}\int \sum_{j\neq i, j\leq n} \left(\log |x-x_j|-\log |x_i-x_j|\right)p(x)dx\\ & \geq & 2 \int \lim\inf_n n^{-1}\sum_{j\neq i, 1,j\leq n} \left(\log |x-x_j|-\log |x_i-x_j|\right)p(x)dx\\ & = & 2\int \left(f_p(x)-f_p(x_i)\right)p(x)dx\\ & = & 2(c_p-f_p(x_i)) \end{eqnarray*} utilizando el lema de Fatou, y luego la ley fuerte de los grandes números.
Tenga en cuenta que $f_p'(x)=H_p(x)$ , donde $H_p(x)=\int {\frac {p(y)}{x-y}}dy$ es la transformada de Hilbert de $p$ . En particular, $f_p$ es continua. Ahora, supongamos que se da el caso de que $f_p(x)$ no es constante. En este caso, existe $\epsilon>0$ así como un intervalo $[a,b]$ tal que $\inf_{x'\in [a,b]}c_p-f_p(x')\geq\epsilon/2$ . Por la suposición (1), $\int_a^b p(x)dx>0$ Así que en consecuencia, $P(\exists i: x_i\in [a,b])=1$ . Podemos reordenar un número finito de muestras sin cambiar la distribución, por lo que en WLOG podemos suponer que $x_1\in [a,b]$ . Como consecuencia $2(c_p-f_p(x_1))\geq\epsilon$ .
Combinando con lo anterior, obtenemos
$$n^{-1}\log \int l_{1,n}^2(x)p(x)dx\geq \epsilon$$ para todos $n$ suficientemente grande. Por simple reordenación, esto es equivalente a $\int l_{1,n}^2(x)p(x)dx\geq e^{n\epsilon}$ y como $V_n=\int \sum_{i=1}^n l_{i,n}^2(x)p(x)dx \geq \int l_{1,n}^2(x)p(x)dx$ obtenemos el pretendido crecimiento exponencial, suponiendo que $f_p$ no es constante .
Así que el último paso es demostrar que $f_p$ no es constante si $p$ no es la distribución del arcoseno, o, de forma equivalente, que la transformada de hilbert $H_p$ no es idéntico a cero. Bajo el supuesto (2), esto se deduce de a teorema de Tricomi . De este modo, hemos demostrado la afirmación.
Código de los polinomios de Bernstein
from scipy.special import binom
import numpy as np
import matplotlib.pyplot as plt
variances=[]
n_range=[5,10,20,50,100]
for n in n_range:
preds=[]
for _ in range(1000):
xs=np.linspace(0,1,100)
X=np.linspace(0,1,n+1)
Y=np.sin(8*X)+np.random.randn(n+1)*.1
nu=np.arange(n+1)
bern=binom(n,nu)[:,None]*(xs[None,:])**(nu[:,None])*(1-xs[None,:])**(n-nu[:,None])*Y[:,None]
pred=bern.sum(0)
preds.append(pred)
preds=np.array(preds)
variances.append(preds.var(0).mean())
variances=np.array(variances)
plt.scatter(n_range,variances)
plt.xlabel('n')
plt.ylabel('variance')
plt.ylim(0,.005)
plt.title('bernstein polynomial variance')
Código para polinomios complejos
n=64
xarg=np.linspace(0,2*np.pi,n+1,endpoint=True)[:n]
y=np.sin(2.5*xarg)+.25*np.random.randn(n)
xcomp=np.exp(1j*xarg)
xs=np.linspace(0,2*np.pi,200)
ys=np.sin(2.5*xs)
X=np.array([xcomp**i for i in range(len(xarg))])
w=np.linalg.solve(X,y)
#only need to go up to middle frequency, bc later ones are complex conjugates
interp=2*np.dot(w[1:n//2+1],np.array([np.exp(1j*xs*i) for i in range(1,n//2+1)]))+w2[0]
plt.plot(xs,interp.real,label='complex polynomial interp')
plt.plot(xs,ys,c='gray',linestyle='--',label='true function')
plt.scatter(xarg,y)
plt.xlabel('complex argument')
plt.legend()
