Regresión Cuantílica GAMLSS
2023-10-04
#1.- INTRODUCCIÓN
La predicción de una variable Y en función de uno o varios predictores X es un problema de aprendizaje supervisado que puede resolverse con múltiples métodos de Machine Learning y aprendizaje estadístico. Algunos de ellos consideran que la relación entre Y y X es únicamente lineal (regresión lineal, GLM), mientras que otros permiten incorporar relaciones no lineales o incluso interacciones entre predictores (SVM, Random Forest, Boosting).
De una forma u otra, todos ellos tratan de inferir la relación entre X e Y para obtener información sobre la distribución condicional de la variable respuesta en función de variables predictoras. Sin embargo, la gran mayoría de los modelos de regresión únicamente modelan la media de la variable respuesta E(Y|X=x), asumiendo que el resto de características de la distribución (dispersión, asimetría…) son constantes. Esto supone una limitación importante a la hora de modelar distribuciones complejas, sobretodo si se pretende predecir intervalos de confianza o cuantiles.
#2.- GAMLSS
Los modelos aditivos generalizados para posición, escala y forma GAMLSS (Generalized Additive Models for Location, Scale and Shape), son modelos de regresión semi-paramétricos. Paramétricos en cuanto a que requieren asumir que la variable respuesta sigue una determinada distribución paramétrica (normal, beta, gamma…) y semi porque los parámetros de esta distribución pueden ser modelados, cada uno de forma independiente, siguiendo funciones no paramétricas (lineales, aditivas o no lineales). Esta versatilidad hace de los GAMLSS una herramienta adecuada para modelar variables que siguen todo un abanico de distribuciones (no normales, asimétricas, con varianza no constante…).
Los modelos GAMLSS asumen que la variable respuesta tiene una función de densidad definida por hasta 4 parámetros (μ,σ,ν,τ) que determinan su posición (p.ej. media), escala (p.ej. desviación estándar) y forma (p. ej. skewness y kurtosis), y que cada uno de ellos puede variar independientemente de los otros en función de los predictores. Estos modelos aprenden por lo tanto hasta 4 funciones, donde cada una establece la relación entre uno de los parámetros y las variables predictoras. Por ejemplo, en la regresión gaussiana, la variable respuesta depende de dos parámetros: la media μ y de la desviación típica σ. En lugar de asumir que σ es constante (como hacen los modelos LM y GAM), los modelos GAMLSS modelan ambos parámetros en función de las variables predictoras.
Los modelos GAMLSS son capaces de caracterizar la distribución completa, permitiendo generar intervalos probabilísticos y predicción de cuantiles, sin tener que asumir que la varianza es constante ni que las relaciones son únicamente lineales.
#3.- DATOS
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'tidyr' was built under R version 4.2.3## Warning: package 'readr' was built under R version 4.2.3## Warning: package 'purrr' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## Warning: package 'forcats' was built under R version 4.2.3## Warning: package 'lubridate' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors#simulación de una variable no uniforme en el rango de X
#=======================================================
set.seed(12345)
n <- 2000
x <- runif(min = 0, max = 24, n = n) %>% sort()
y <- rnorm(
n,
mean = 10,
sd = 1 + 1.5*(4.8 < x & x < 7.2) + 4*(7.2 < x & x < 12) + 1.5*(12 < x & x < 14.4) + 2*(x > 16.8)
)# Calculamos el cuantil 10
#================================
cuantil_10 <- qnorm(
p=0.1,
mean = 10,
sd = 1 + 1.5*(4.8 < x & x < 7.2) + 4*(7.2 < x & x < 12) + 1.5*(12 < x & x < 14.4) + 2*(x > 16.8)
)
# Calculamos el cuantil 90
#================================
cuantil_90 <- qnorm(
p=0.9,
mean = 10,
sd = 1 + 1.5*(4.8 < x & x < 7.2) + 4*(7.2 < x & x < 12) + 1.5*(12 < x & x < 14.4) + 2*(x > 16.8)
)datos <- data.frame(consumo = y, hora = x, cuantil_10, cuantil_90)head(datos)## consumo hora cuantil_10 cuantil_90
## 1 11.677512 0.01237827 8.718448 11.28155
## 2 10.079474 0.02727808 8.718448 11.28155
## 3 9.143573 0.06556776 8.718448 11.28155
## 4 9.221223 0.10210288 8.718448 11.28155
## 5 9.619064 0.11792478 8.718448 11.28155
## 6 8.102642 0.12558777 8.718448 11.28155summary(datos)## consumo hora cuantil_10 cuantil_90
## Min. :-5.631 Min. : 0.01238 Min. :3.592 Min. :11.28
## 1st Qu.: 8.441 1st Qu.: 6.49881 1st Qu.:6.155 1st Qu.:11.28
## Median : 9.965 Median :12.36921 Median :6.155 Median :13.84
## Mean : 9.955 Mean :12.33160 Mean :6.491 Mean :13.51
## 3rd Qu.:11.392 3rd Qu.:18.33743 3rd Qu.:8.718 3rd Qu.:13.84
## Max. :23.805 Max. :23.98961 Max. :8.718 Max. :16.41# Restricción de no negatividad
#=======================================
datos <- datos %>%
filter(consumo >= 0)# Restricción para pertenencia al cuantil 10-90
# ===================================================
datos <- datos %>%
mutate(dentro_intervalo = ifelse(
consumo > cuantil_10 & consumo < cuantil_90,
TRUE,
FALSE
)
)str(datos)## 'data.frame': 1990 obs. of 5 variables:
## $ consumo : num 11.68 10.08 9.14 9.22 9.62 ...
## $ hora : num 0.0124 0.0273 0.0656 0.1021 0.1179 ...
## $ cuantil_10 : num 8.72 8.72 8.72 8.72 8.72 ...
## $ cuantil_90 : num 11.3 11.3 11.3 11.3 11.3 ...
## $ dentro_intervalo: logi FALSE TRUE TRUE TRUE TRUE FALSE ...head(datos)## consumo hora cuantil_10 cuantil_90 dentro_intervalo
## 1 11.677512 0.01237827 8.718448 11.28155 FALSE
## 2 10.079474 0.02727808 8.718448 11.28155 TRUE
## 3 9.143573 0.06556776 8.718448 11.28155 TRUE
## 4 9.221223 0.10210288 8.718448 11.28155 TRUE
## 5 9.619064 0.11792478 8.718448 11.28155 TRUE
## 6 8.102642 0.12558777 8.718448 11.28155 FALSEsummary(datos)## consumo hora cuantil_10 cuantil_90
## Min. : 0.648 Min. : 0.01238 Min. :3.592 Min. :11.28
## 1st Qu.: 8.465 1st Qu.: 6.46013 1st Qu.:6.155 1st Qu.:11.28
## Median : 9.974 Median :12.41007 Median :6.155 Median :13.84
## Mean :10.021 Mean :12.34468 Mean :6.505 Mean :13.49
## 3rd Qu.:11.404 3rd Qu.:18.36864 3rd Qu.:8.718 3rd Qu.:13.84
## Max. :23.805 Max. :23.98961 Max. :8.718 Max. :16.41
## dentro_intervalo
## Mode :logical
## FALSE:373
## TRUE :1617
##
##
## p <- ggplot() +
geom_point(
data = datos,
aes(x = hora, y = consumo),
alpha = 0.3,
color = "blue") +
geom_line(
data = datos,
aes(x = hora , y = 10, color = "media"),
linetype = "solid",
size = 1) +
scale_color_manual(name = "",
breaks = c("media"),
values = c("media" = "red")) +
labs(title = "Evolución del consumo eléctrico a lo largo del día",
x = "Hora del día",
y = "Consumo eléctrico (MWh)") +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(face = "bold"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p <- p +
scale_x_continuous(breaks = seq(0, 24, 2),
labels = seq(0, 24, 2))
p
La media del consumo eléctrico es la misma que durante todo el día, consumo = 10 Mwh; sin embargo, su dispersión no es constante (heterocedasticidad). Véase el resultado de predecir el consumo medio en función de la hora del día con un modelo Random Forest.
El valor predicho es muy próximo a la media real, es decir, el modelo es bueno prediciendo el consumo medio esperado. Ahora, imagínese que la compañía encargada de suministrar la electricidad debe de ser capaz de provisionar, en un momento dado, con hasta un 50% de electricidad extra respecto al promedio. Esto significa un máximo de 15 Mwh. Estar preparado para suministrar este extra de energía implica gastos de personal y maquinaría, por lo que la compañía se pregunta si es necesario estar preparado para producir tal cantidad durante todo el día, o si, por lo contrario, podría evitarse durante algunas horas, ahorrando así gastos.
Un modelo que predice únicamente el promedio no permite responder a esta pregunta, ya que tanto para las 2h de la mañana como para las 8h, el consumo promedio predicho es en torno a 10 Mwh, sin embargo, la probabilidad de que se alcancen consumos de 15 Mwh a las 2h es prácticamente nula mientras que esto ocurra a las 8h sí es razonable.
Una forma de describir la dispersión de una variable es el uso de cuantiles. El cuantil de orden τ (0<τ<1) de una variable continua Y es el valor yτ tal que, una proporción τ de valores de la población, es menor o igual que dicho valor (Prob(Y≤yτ)=τ). Por ejemplo, el cuantil de orden 0.36 deja un 36% de valores por debajo y el cuantil de orden 0.50 el 50% (se corresponde con la mediana de la distribución).
Dado que los datos se han simulado empleando distribuciones normales, se conoce el valor de los cuantiles teóricos para cada X. Se muestra de nuevo el mismo gráfico pero esta vez añadiendo los cuantiles 0.1 y 0.9.
Si como resultado del modelo, además de la predicción de la media, se predice también el valor de los cuantiles, se dispone de una caracterización mayor de la distribución de la variable respuesta, y con ello se puede responder a más preguntas. Por ejemplo, en el caso de la energía, se tendría cierta seguridad al decir que, durante los intervalos de 0h a 4h y de 12h a 14h, es poco probable que se alcancen consumos de 15 Mwh
Otros casos en los que conocer la distribución de cuantiles puede ser útil son:
*Identificación de regiones en las que la variable respuesta tiene mayor dispersión en torno a su media.
*Entrenar modelos que predicen la mediana (cuantil 0.5) en lugar de la media. Estos modelos son más robustos frente a outliers.
*Detectar anomalías, identificando aquellas observaciones que están fuera de un determinado intervalo cuantílico.
En este ejemplo se describe cómo los modelos GAMLSS permiten caracterizar toda la distribución y por lo tanto predecir los cuantiles.
4.- MODELO
El primer paso para ajustar un modelo GAMLSS es identificar qué tipo de distribución paramétrica sigue la variable respuesta. En este caso, como los datos han sido simulados, se sabe que siguen una distribución normal. Aunque para este primer ejemplo ilustrativo se asume como cierto, en la práctica, esta información se desconoce, por lo que es necesario un primer estudio de la distribución.
hist(datos$consumo)
qqnorm(datos$consumo)
qqline(datos$consumo)
library(gamlss)## Warning: package 'gamlss' was built under R version 4.2.3## Loading required package: splines## Loading required package: gamlss.data##
## Attaching package: 'gamlss.data'## The following object is masked from 'package:datasets':
##
## sleep## Loading required package: gamlss.dist## Warning: package 'gamlss.dist' was built under R version 4.2.3## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## Loading required package: parallel## ********** GAMLSS Version 5.4-18 **********## For more on GAMLSS look at https://www.gamlss.com/## Type gamlssNews() to see new features/changes/bug fixes.modelo <- gamlss(
formula = consumo ~ pb(hora),
sigma.formula = ~ pb(hora),
family = NO,
data = datos,
ccontrol = gamlss.control(trace = FALSE)
)## GAMLSS-RS iteration 1: Global Deviance = 8985.282
## GAMLSS-RS iteration 2: Global Deviance = 8980.298
## GAMLSS-RS iteration 3: Global Deviance = 8980.308
## GAMLSS-RS iteration 4: Global Deviance = 8980.308summary(modelo)## ******************************************************************
## Family: c("NO", "Normal")
##
## Call: gamlss(formula = consumo ~ pb(hora), sigma.formula = ~pb(hora),
## family = NO, data = datos, ccontrol = gamlss.control(trace = FALSE))
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.9549493 0.0591163 168.396 <2e-16 ***
## pb(hora) -0.0008836 0.0053652 -0.165 0.869
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.481390 0.032502 14.81 <2e-16 ***
## pb(hora) 0.028841 0.002298 12.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1990
## Degrees of Freedom for the fit: 18.02599
## Residual Deg. of Freedom: 1971.974
## at cycle: 4
##
## Global Deviance: 8980.308
## AIC: 9016.36
## SBC: 9117.232
## ******************************************************************#5.- PREDICCIÓN
# Se predice todo el rango de hora para representar los cuantiles
grid_predictor <- seq(0, 24, length.out = 2500)
predicciones <- predictAll(
modelo,
newdata = data.frame(hora = grid_predictor),
type = "response"
)## new prediction
## new predictionpredicciones <- as.data.frame(predicciones)
predicciones <- bind_cols(data.frame(hora = grid_predictor), predicciones)
predicciones %>% head (n=5)## hora mu sigma
## 1 0.000000000 9.954907 0.9048709
## 2 0.009603842 9.954899 0.9065713
## 3 0.019207683 9.954891 0.9082748
## 4 0.028811525 9.954882 0.9099789
## 5 0.038415366 9.954874 0.9116814Una vez predichos los parámetros que caracterizan a la distribución en función del predictor, se puede calcular la probabilidad de cada observación o el intervalo que acumula un determinado porcentaje de probabilidad (intervalo cuantílico). Todas las distribuciones de los paquetes gamlss y gamlss.dist disponen de funciones d, p, q y r para calcular probabilidad, densidad, cuantiles y generar valores aleatorios.
# Cálculo de los cuantiles teóricos para establecer el intervalo central que acumula
# un 90% de probabilidad empleando los parámetros predichos.
predicciones <- predicciones %>%
mutate(
cuantil_10_pred = purrr::pmap_dbl(
.l = list(mu = mu, sigma = sigma),
.f = function(mu, sigma) {qNO(p = 0.1, mu, sigma)}
),
cuantil_90_pred = purrr::pmap_dbl(
.l = list(mu = mu, sigma = sigma),
.f = function(mu, sigma) {qNO(p = 0.9, mu, sigma)}
)
)p <- ggplot() +
geom_ribbon(
data = datos,
aes(x = hora, ymin = cuantil_10, ymax = cuantil_90),
fill = "red",
alpha = 0.2) +
geom_point(
data = datos,
aes(x = hora, y = consumo),
alpha = 0.2,
color = "black") +
geom_line(
data = predicciones,
aes(x = hora, y = cuantil_10_pred, color = "prediccion_cuantil"),
size = 0.7) +
geom_line(
data = predicciones,
aes(x = hora, y = cuantil_90_pred, color = "prediccion_cuantil"),
size = 0.7) +
geom_line(
data = predicciones,
aes(x = hora, y = mu, color = "prediccion_media"),
size = 1) +
scale_color_manual(name = "",
breaks = c("prediccion_cuantil", "prediccion_media"),
values = c("prediccion_cuantil" = "blue",
"prediccion_media" = "firebrick")) +
labs(title = "Evolución del consumo eléctrico a lo largo del día",
subtitle = "Intervalo real entre cuantiles 0.1 y 0.9 sombreado en rojo",
x = "Hora del día",
y = "Consumo eléctrico (MWh)") +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(face = "bold"))
p <- p +
scale_x_continuous(breaks = seq(0, 24, 2),
labels = seq(0, 24, 2))p
Si, por ejemplo, se desea conocer la probabilidad de que a las 8h el
consumo supere los 15 Mwh, primero se predicen los parámetros de la
distribución a para (hora=8) y después se calcula la
probabilidad de consumo≥15 con su función de
distribución Normal.
prediccion <- predictAll(
modelo,
newdata = data.frame(hora=8),
type = "response"
)## new prediction
## new predictionprediccion## $mu
## [1] 9.947882
##
## $sigma
## [1] 4.607569
##
## attr(,"family")
## [1] "NO" "Normal"# Se calcula la probabilidad
probabilidad_consumo <- pNO(
q = 15,
mu = prediccion$mu,
sigma = prediccion$sigma,
lower.tail = FALSE
)
probabilidad_consumo## [1] 0.1364339Acorde al modelo, la probabilidad de que a las 8h el consumo sea igual o superior a 15 Mwh es del 13.6%.
Probabilidad h = 0
prediccion0 <- predictAll(
modelo,
newdata = data.frame(hora = 0),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo0 <- pNO(
q = 15,
mu = prediccion0$mu,
sigma = prediccion0$sigma,
lower.tail = FALSE
)
probabilidad_consumo0## [1] 1.234218e-08prediccion1 <- predictAll(
modelo,
newdata = data.frame(hora = 1),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo1 <- pNO(
q = 15,
mu = prediccion1$mu,
sigma = prediccion1$sigma,
lower.tail = FALSE
)
probabilidad_consumo1## [1] 5.477145e-07prediccion2 <- predictAll(
modelo,
newdata = data.frame(hora = 2),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo2 <- pNO(
q = 15,
mu = prediccion2$mu,
sigma = prediccion2$sigma,
lower.tail = FALSE
)
probabilidad_consumo2## [1] 1.305351e-07prediccion3 <- predictAll(
modelo,
newdata = data.frame(hora = 3),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo3 <- pNO(
q = 15,
mu = prediccion3$mu,
sigma = prediccion3$sigma,
lower.tail = FALSE
)
probabilidad_consumo3## [1] 2.431742e-08prediccion4 <- predictAll(
modelo,
newdata = data.frame(hora = 4),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo4 <- pNO(
q = 15,
mu = prediccion4$mu,
sigma = prediccion4$sigma,
lower.tail = FALSE
)
probabilidad_consumo4## [1] 3.249909e-06prediccion5 <- predictAll(
modelo,
newdata = data.frame(hora = 5),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo5 <- pNO(
q = 15,
mu = prediccion5$mu,
sigma = prediccion5$sigma,
lower.tail = FALSE
)
probabilidad_consumo5## [1] 0.001423056prediccion6 <- predictAll(
modelo,
newdata = data.frame(hora = 6),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo6 <- pNO(
q = 15,
mu = prediccion6$mu,
sigma = prediccion6$sigma,
lower.tail = FALSE
)
probabilidad_consumo6## [1] 0.02525045prediccion7 <- predictAll(
modelo,
newdata = data.frame(hora = 7),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo7 <- pNO(
q = 15,
mu = prediccion7$mu,
sigma = prediccion7$sigma,
lower.tail = FALSE
)
probabilidad_consumo7## [1] 0.08066608prediccion8 <- predictAll(
modelo,
newdata = data.frame(hora = 8),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo8 <- pNO(
q = 15,
mu = prediccion8$mu,
sigma = prediccion8$sigma,
lower.tail = FALSE
)
probabilidad_consumo8## [1] 0.1364339prediccion9 <- predictAll(
modelo,
newdata = data.frame(hora = 9),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo9 <- pNO(
q = 15,
mu = prediccion9$mu,
sigma = prediccion9$sigma,
lower.tail = FALSE
)
probabilidad_consumo9## [1] 0.1433312prediccion10 <- predictAll(
modelo,
newdata = data.frame(hora = 10),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo10 <- pNO(
q = 15,
mu = prediccion10$mu,
sigma = prediccion10$sigma,
lower.tail = FALSE
)
probabilidad_consumo10## [1] 0.1446728prediccion11 <- predictAll(
modelo,
newdata = data.frame(hora = 11),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo11 <- pNO(
q = 15,
mu = prediccion11$mu,
sigma = prediccion11$sigma,
lower.tail = FALSE
)
probabilidad_consumo11## [1] 0.1378342prediccion12 <- predictAll(
modelo,
newdata = data.frame(hora = 12),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo12 <- pNO(
q = 15,
mu = prediccion12$mu,
sigma = prediccion12$sigma,
lower.tail = FALSE
)
probabilidad_consumo12## [1] 0.07187003prediccion13 <- predictAll(
modelo,
newdata = data.frame(hora = 13),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo13 <- pNO(
q = 15,
mu = prediccion13$mu,
sigma = prediccion13$sigma,
lower.tail = FALSE
)
probabilidad_consumo13## [1] 0.02786572prediccion14 <- predictAll(
modelo,
newdata = data.frame(hora = 14),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo14 <- pNO(
q = 15,
mu = prediccion14$mu,
sigma = prediccion14$sigma,
lower.tail = FALSE
)
probabilidad_consumo14## [1] 0.005508971prediccion15 <- predictAll(
modelo,
newdata = data.frame(hora = 15),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo15 <- pNO(
q = 15,
mu = prediccion15$mu,
sigma = prediccion15$sigma,
lower.tail = FALSE
)
probabilidad_consumo15## [1] 1.443351e-05prediccion16 <- predictAll(
modelo,
newdata = data.frame(hora = 16),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo16 <- pNO(
q = 15,
mu = prediccion16$mu,
sigma = prediccion16$sigma,
lower.tail = FALSE
)
probabilidad_consumo16## [1] 6.424849e-07prediccion17 <- predictAll(
modelo,
newdata = data.frame(hora = 17),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo17 <- pNO(
q = 15,
mu = prediccion17$mu,
sigma = prediccion17$sigma,
lower.tail = FALSE
)
probabilidad_consumo17## [1] 0.003023241prediccion18 <- predictAll(
modelo,
newdata = data.frame(hora = 18),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo18 <- pNO(
q = 15,
mu = prediccion18$mu,
sigma = prediccion18$sigma,
lower.tail = FALSE
)
probabilidad_consumo18## [1] 0.0386389prediccion19 <- predictAll(
modelo,
newdata = data.frame(hora = 19),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo19 <- pNO(
q = 15,
mu = prediccion19$mu,
sigma = prediccion19$sigma,
lower.tail = FALSE
)
probabilidad_consumo19## [1] 0.04613378prediccion20 <- predictAll(
modelo,
newdata = data.frame(hora = 20),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo20 <- pNO(
q = 15,
mu = prediccion20$mu,
sigma = prediccion20$sigma,
lower.tail = FALSE
)
probabilidad_consumo20## [1] 0.03533945prediccion21 <- predictAll(
modelo,
newdata = data.frame(hora = 21),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo21 <- pNO(
q = 15,
mu = prediccion21$mu,
sigma = prediccion21$sigma,
lower.tail = FALSE
)
probabilidad_consumo21## [1] 0.03236875prediccion22 <- predictAll(
modelo,
newdata = data.frame(hora = 22),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo22 <- pNO(
q = 15,
mu = prediccion22$mu,
sigma = prediccion22$sigma,
lower.tail = FALSE
)
probabilidad_consumo22## [1] 0.03808713prediccion23 <- predictAll(
modelo,
newdata = data.frame(hora = 23),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo23 <- pNO(
q = 15,
mu = prediccion23$mu,
sigma = prediccion23$sigma,
lower.tail = FALSE
)
probabilidad_consumo23## [1] 0.048114prediccion24 <- predictAll(
modelo,
newdata = data.frame(hora = 24),
type = "response"
)## new prediction
## new predictionprobabilidad_consumo24 <- pNO(
q = 15,
mu = prediccion24$mu,
sigma = prediccion24$sigma,
lower.tail = FALSE
)
probabilidad_consumo24## [1] 0.04166023probabilidad_consumo0## [1] 1.234218e-08probabilidad_consumo1## [1] 5.477145e-07probabilidad_consumo2## [1] 1.305351e-07probabilidad_consumo3## [1] 2.431742e-08probabilidad_consumo4## [1] 3.249909e-06probabilidad_consumo5## [1] 0.001423056probabilidad_consumo6## [1] 0.02525045probabilidad_consumo7## [1] 0.08066608probabilidad_consumo8## [1] 0.1364339probabilidad_consumo9## [1] 0.1433312probabilidad_consumo10## [1] 0.1446728probabilidad_consumo11## [1] 0.1378342probabilidad_consumo12## [1] 0.07187003probabilidad_consumo13## [1] 0.02786572probabilidad_consumo14## [1] 0.005508971probabilidad_consumo15## [1] 1.443351e-05probabilidad_consumo16## [1] 6.424849e-07probabilidad_consumo17## [1] 0.003023241probabilidad_consumo18## [1] 0.0386389probabilidad_consumo19## [1] 0.04613378probabilidad_consumo20## [1] 0.03533945probabilidad_consumo21## [1] 0.03236875probabilidad_consumo22## [1] 0.03808713probabilidad_consumo23## [1] 0.048114probabilidad_consumo24## [1] 0.04166023