En primer lugar, escribamos el modelo sin restricciones de forma que los coeficientes estén ordenados de forma coherente (en este caso, de menor a mayor grado):
unconstrained_model <- lm(y ~ x + I(x^2) +I(x^3))
En segundo lugar, la regresión por mínimos cuadrados con restricciones siempre tiene un RMSE más alto en la muestra de datos que los mínimos cuadrados sin restricciones, a menos que estos últimos cumplan las restricciones, en cuyo caso las dos soluciones son iguales. En la práctica, si los datos no respetan muy bien las restricciones, el ajuste restringido puede ser bastante malo.
En su caso, los datos están disminuyendo, lo cual es bueno. Pero también tienes 4 puntos de datos en los que la ordenada es casi constante (con respecto al rango de y, $I_y =[\min(y),\max(y)]$ ), aunque la variación correspondiente de x sea considerable, con respecto a $I_x = [\min(x),\max(x)]$ En otras palabras, tiene un gran meseta en el centro de su muestra:
df <- data.frame(x,y)
library(ggplot2)
(p <- ggplot(data=df,aes(x=x,y=y)) + geom_line()+geom_point()+geom_smooth())
![enter image description here]()
No se trata de un patrón que pueda ajustarse fácilmente mediante un polinomio monótono decreciente de tercer grado.
Dicho esto, vamos a resolver el problema de regresión restringida. En primer lugar, obtengamos ecuaciones explícitas para las restricciones en términos de los coeficientes del modelo. Ya has observado que el polinomio de tercer grado es decreciente en $I=[0,1]$ sólo si
$$3ax^2+2bx+c <0 \quad \forall x \in I $$
Tenga en cuenta que $d$ no aparece en esta desigualdad, y la razón es obvia - si nuestro modelo (el polinomio de tercer grado) es creciente o no, no depende del valor del intercepto. Puesto que $3ax^2+2bx+c$ debe ser negativo en 0, obtenemos una de las restricciones como $c<0$ . Ahora, para obtener las otras desigualdades de restricción, sólo tenemos que hacer las sustituciones
$$ t_1=x,\quad t_2=x^2$$
y observe que
$$x\in[0,1]\Rightarrow(t_1,t_2)\in [0,1]\times [0,1]$$
A continuación, se nos plantea el problema más sencillo de imponer una restricción de negatividad a un polinomio lineal (de grado uno) en dos variables:
$$f(t_1,t_2)=3at_2+2bt_1+c <0 \quad \forall (t_1,t_2) \in [0,1]\times[0,1] $$
$f$ es lineal y su dominio, el cuadrado unitario $[0,1]^2$ es convexa. Entonces si $f$ es positivo en las cuatro esquinas $(0,0)$ , $(0,1)$ , $(1,0)$ , $(1,1)$ es positivo en todas las zonas de $[0,1]^2$ . En realidad, los únicos valores posibles para $(t_1,t_2)$ son los de la mitad inferior derecha del cuadrado (el triángulo con ángulos $(0,0)$ , $(0,1)$ , $(1,1)$ ), porque
$$x \in [0,1] \Rightarrow x^2 \le x \Leftrightarrow t_2 \le t_1 $$
Por lo tanto, sólo tenemos que imponer que $f$ es negativo en las tres esquinas $(0,0)$ , $(0,1)$ , $(1,1)$ . La condición de negatividad en $(0,0)$ da de nuevo $c<0$ . Imponiéndolo también en $(1,0)$ y $(1,1)$ obtenemos finalmente las tres condiciones
$$c<0,\quad3a+c <0, \quad 3a+2b+c <0$$
Perfecto Ahora tenemos tres restricciones de desigualdad lineal que expresan la condición de que nuestro modelo es negativo en $I$ . Para resolver un problema lineal de mínimos cuadrados con restricciones de igualdad y/o desigualdad, R ofrece un paquete con una interfaz sencilla e intuitiva, limSolve .
library(limSolve)
A <- cbind(rep(1,length(x)),x,x^2,x^3)
b <- y
G <- matrix(nrow=3,ncol=4,byrow = TRUE,data = c(0, -1,-2,-3,0,-1,-2,0,0,-1,0,0))
h <- rep(0,3)
constrained_model <- lsei(A = A, B = b, G = G, H = h, type=2)
Genial, ahora hagamos predicciones para ambos modelos. lm tienen un predict pero no se dispone de un método equivalente para los resultados de lsei . Por lo tanto, utilizaremos un enfoque diferente y definiremos un my_predict función
my_predict <- function(x,coefficients){
X <- cbind(rep(1,length(x)),x,x^2,x^3)
predictions <- X%*%coefficients
}
Entonces
# compute predictions
xpred <- seq(0,1,len=100)
predictions_constrained <- my_predict(xpred,constrained_model$X)
predictions_unconstrained <- my_predict(xpred,unconstrained_model$coefficients)
df2 <- data.frame(xpred,predictions_unconstrained,predictions_constrained)
# plot results
p <- ggplot(data = df,aes(x = x, y = y,color = "data")) +
geom_point() +
geom_line(data = df2, aes(x = xpred, y = predictions_unconstrained, color = "unconstrained fit")) +
geom_line(data = df2, aes(x = xpred, y = predictions_constrained, color = "constrained fit"))
p
![enter image description here]()
Como era de esperar, el ajuste restringido es considerablemente peor que el no restringido. También podemos calcular el $R^2$ y obtener la misma imagen:
SS <-sum((y-mean(y))^2)
RSS_unconstrained <- sum((unconstrained_model$residuals)^2)
RSS_constrained <- sum((A%*%constrained_model$X-y)^2)
R2_constrained <- 1-RSS_constrained/SS
R2_unconstrained <- 1-RSS_unconstrained/SS
El modelo restringido sólo explica alrededor del 70% de la variación total, mientras que el modelo no restringido explica alrededor del 83% de la variación total. Como nota final, el hecho de que el ajuste sea peor en los datos de la muestra no significa necesariamente que el modelo restringido sea peor que el no restringido, en términos de precisión predictiva. Si, por ejemplo, el modelo ideal obedece a las restricciones, puede que el modelo restringido tenga un error de prueba menor en los datos no vistos.
En cuanto al nuevo conjunto de datos:
x <- c(0.01041667, 0.30208333, 0.61458333, 0.65625000, 0.83333333)
y <- c(772, 607, 576, 567, 550)
Repitiendo exactamente los mismos cálculos se obtiene un modelo restringido cuyo término de mayor grado es igual a cero:
> constrained_model
$X
x
759.8756 -449.1166 224.5583 0.0000
$residualNorm
[1] 0
$solutionNorm
[1] 1854.049
$IsError
[1] FALSE
$type
[1] "lsei"
Este polinomio es no creciente en $I=[0,1]$ , según se requiera, pero tiene un mínimo en $x=1$ por lo que empezaría a aumentar para $x>1$ . De hecho, la derivada es $2bx+c$ (ya que $a=0$ ), que para $x=1$ se convierte en
> 2*constrained_model$X[3]+constrained_model$X[2]
-6.82121e-13
que, dada la precisión de los cálculos implicados, es básicamente 0. El hecho de que $x=1$ es un mínimo también se desprende del gráfico: ![enter image description here]()
Este
