Este analisis se realizo con datos de 93 dias, desde que se detecto la primera infeccion por Covid19. Los modelos usan x=dia y=numero de infectados. El calculo se hizo para estimar el numero de posibles infectados para el dia 30 de mayo del 2020. Los datos se obtuvieron de: https://www.unionguanajuato.mx/articulo/2020/03/22/cultura/casos-de-coronavirus-en-mexico-por-estado-estadisticas-covid-19
**Nota*.Este es un analisis propio del autor, no es oficial, solo son estimaciones.
El modelo: \(y=\dfrac{a\mathrm{e}^{cx}}{\mathrm{e}^{b-x}+1}\) Su primera derivada: \(\frac{dy}{dx}= \dfrac{a\mathrm{e}^{cx-x+b}}{\left(\mathrm{e}^{b-x}+1\right)^2}+\dfrac{ac\mathrm{e}^{cx}}{\mathrm{e}^{b-x}+1}\) Simplificada se reescribe: \(\frac{dy}{dx}= \dfrac{a\left(c\mathrm{e}^x+\mathrm{e}^bc+\mathrm{e}^b\right)\mathrm{e}^{cx+x}}{\left(\mathrm{e}^x+\mathrm{e}^b\right)^2}\)
# Modelo 1 ----------------------------------------------------------
n=length(x)
plot(x, y, col="lightblue", pch=19, xlab = "Dia numero", ylab = "Numero de infectados", cex=2, main="Prediccion del numero de infectados con el modelo uno")
abline(v = 31.3, col = "black", lwd = 1, lty = 2)
abline(v = 38.3, col = "black", lwd = 1, lty = 2)
Fit1<- nlsLM(y ~ a0/(1 + exp((a1 - log(x))/a2)),
data = datos,
start = list(a0 = 80000, a1 = 4, a2 = 1),
algorithm = "port",
control = list(maxiter = 500))
#text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=1)
#text(15, 1500, "Estimados para el dia 05 de abril del 2020 = 2041", col = "red")
summary(Fit1)
##
## Formula: y ~ a0/(1 + exp((a1 - log(x))/a2))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a0 5.929e+05 4.815e+04 12.31 <2e-16 ***
## a1 4.950e+00 2.428e-02 203.87 <2e-16 ***
## a2 2.320e-01 1.527e-03 151.95 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 260.7 on 90 degrees of freedom
##
## Algorithm "port", convergence message: Relative error in the sum of squares is at most `ftol'.
#confint(Fit1)
# Add in predicted values into the plot
lines(x, fitted.values(Fit1), lwd=3, col="red")
AIC(Fit1)
## [1] 1303.663
BIC(Fit1)
## [1] 1313.794
#Distribucion de residuales
#plot(fitted(Fit1), residuals(Fit1), xlab="Valores estimados", ylab = "Residuals",
#abline(a=0, b=0), col="red", pch=16)
#Analisis de autocorrelacion de los residuales
#acf(residuales_Fit1<-(y-fitted(Fit1)))
#Osbervados vs estimados
#plot(fitted(Fit1), y, abline(a=0, b=1), col="red", xlab="Valores estimados", ylab ="Valores observados", pch=16)
#Los coeficientes de regresion
coef(summary(Fit1))
## Estimate Std. Error t value Pr(>|t|)
## a0 5.929251e+05 4.815294e+04 12.31337 5.272824e-21
## a1 4.949551e+00 2.427829e-02 203.86738 9.581088e-122
## a2 2.320081e-01 1.526856e-03 151.95148 2.728211e-110
df<-as.data.frame(coef(summary(Fit1)))
#Predicciones cuando x sea 26
x0<-n+1
df[1, 1]
## [1] 592925.1
df[2, 1]
## [1] 4.949551
df[3, 1]
## [1] 0.2320081
#Numero de posibles infectados estimados por el modelo uno.
Prediccion_m1<- (df[1, 1]/(1 + exp((df[2, 1] - log(x0))/df[3, 1])))
Eq1<-round(Prediccion_m1)
this_day <- today()+1
hoy <- Sys.Date();
mnn<-(hoy+1);
mnnn<-format(mnn, format="%d %B %Y")
text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=2)
text(7, (Eq1-5350), labels = Eq1, col = "black", cex=1.8)
text(17, (Eq1-11800), "Posibles infectados para el dia", col = "black")
text(45, (Eq1-11800), labels = mnnn, col = "black")
#Bondad de ajuste del modelo
#La R2
df.residual(Fit1)
## [1] 90
CMErr_Fit1<-sum((y-fitted(Fit1))^2)/(df.residual(Fit1))
CMTot_Fit1<-(sum((y-mean(y))^2))/((length(y)-1))
R2_Fit1<-1-(CMErr_Fit1)/(CMTot_Fit1)
#El Error
Sxy_Fit1<-sqrt(CMErr_Fit1)
#El Coeficiente de variacion
CV_Fit1<-(Sxy_Fit1)/(mean(y))*100
#Intervalos de confianza (95 %) Simulacio Monte Carlo modelo
CONC_Fit1 <- seq(n+1, n+1, by = 1)
Pred_Fit1 <- predict(Fit1, newdata = data.frame(x = CONC_Fit1), nsim = 100000)
IC_95_Fit1 <- predictNLS(Fit1, newdata = data.frame(x = CONC_Fit1), nsim = 100000)
## predictNLS: Propagating predictor value #1...
El modelo: \(y=\dfrac{a}{\mathrm{e}^{cx+b}+1}\) Su primera derivada: \(\frac{dy}{dx}= -\dfrac{ac\mathrm{e}^{cx+b}}{\left(\mathrm{e}^{cx+b}+1\right)^2}\)
# Modelo 2 ----------------------------------------------------------
plot(x, y, col="pink", pch=19, xlab = "Dia numero", ylab = "Numero de infectados", cex=2,
main="Prediccion del numero de infectados con el modelo dos")
abline(v = 31.3, col = "black", lwd = 1, lty = 2)
abline(v = 38.3, col = "black", lwd = 1, lty = 2)
Fit2 <- nlsLM(y ~ b0/(1 + exp(b1 + b2*x)),
start=list(b0=12000, b1=8, b2=-0.5), trace=FALSE,
algorithm = "port")
#text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=1)
#text(15, 1500, "Estimados para el dia 05 de abril del 2020 = 1986", col = "black")
summary(Fit2)
##
## Formula: y ~ b0/(1 + exp(b1 + b2 * x))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 1.475e+05 4.654e+03 31.69 <2e-16 ***
## b1 6.817e+00 4.661e-02 146.25 <2e-16 ***
## b2 -7.599e-02 1.100e-03 -69.06 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 740.7 on 90 degrees of freedom
##
## Algorithm "port", convergence message: Relative error in the sum of squares is at most `ftol'.
#confint(Fit2)
# Add in predicted values into the plot
lines(x, fitted.values(Fit2), lwd=3, col="black")
#AIC(Fit2)
#BIC(Fit2)
#Distribucion de residuales
#plot(fitted(Fit2), residuals(Fit2), xlab="Valores estimados", ylab = "Residuals",
#abline(a=0, b=0), col="black", pch=16)
#Analisis de autocorrelacion de los residuales
#acf(residuales_Fit2<-(y-fitted(Fit2)))
#Osbervados vs estimados
#plot(fitted(Fit2), y, abline(a=0, b=1), col="black", xlab="Valores estimados", ylab ="Valores observados", pch=16)
#Los coeficientes de regresion
coef(summary(Fit2))
## Estimate Std. Error t value Pr(>|t|)
## b0 1.474657e+05 4.653519e+03 31.68907 1.328727e-50
## b1 6.816862e+00 4.661250e-02 146.24538 8.432536e-109
## b2 -7.599455e-02 1.100391e-03 -69.06140 9.374874e-80
df2<-as.data.frame(coef(summary(Fit2)))
#Predicciones cuando x sea 26
x1<-n+1
df2[1, 1]
## [1] 147465.7
df2[2, 1]
## [1] 6.816862
df2[3, 1]
## [1] -0.07599455
#Numero de posibles infectados estimados por el modelo dos.
Prediccion_m2 <- df2[1, 1]/(1 + exp(df2[2, 1] + df2[3, 1]*x1))
Eq1<-round(Prediccion_m2)
this_day <- today()+1
hoy <- Sys.Date();
mnn<-(hoy+1);
mnnn<-format(mnn, format="%d %B %Y")
text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=2)
text(7, (Eq1-5350), labels = Eq1, col = "black", cex=1.8)
text(17, (Eq1-11800), "Posibles infectados para el dia", col = "black")
text(45, (Eq1-11800), labels = mnnn, col = "black")
#Bondad de ajuste del modelo
#La R2
df.residual(Fit2)
## [1] 90
CMErr_Fit2<-sum((y-fitted(Fit2))^2)/(df.residual(Fit2))
CMTot_Fit2<-(sum((y-mean(y))^2))/((length(y)-1))
R2_Fit2<-1-(CMErr_Fit2)/(CMTot_Fit2)
#El Error
Sxy_Fit2<-sqrt(CMErr_Fit2)
#El Coeficiente de variacion
CV_Fit2<-(Sxy_Fit2)/(mean(y))*100
#Intervalos de confianza (95 %) Simulacio Monte Carlo modelo
CONC_Fit2 <- seq(n+1, n+1, by = 1)
Pred_Fit2 <- predict(Fit2, newdata = data.frame(x = CONC_Fit2), nsim = 100000)
IC_95_Fit2 <- predictNLS(Fit2, newdata = data.frame(x = CONC_Fit2), nsim = 100000)
## predictNLS: Propagating predictor value #1...
El modelo: \(y=\mathrm{e}^{\frac{b}{x}+a}\) Su primera derivada: \(\frac{dy}{dx}= -\dfrac{b\mathrm{e}^{\frac{b}{x}+a}}{x^2}\)
# Modelo 3 ----------------------------------------------------------
plot(x, y, col="lightskyblue", pch=19, xlab = "Dia numero", ylab = "Numero de infectados", cex=2, main="Prediccion del numero de infectados con el modelo tres")
abline(v = 0, col = "black", lwd = 1, lty = 2)
abline(v = 31.3, col = "black", lwd = 1, lty = 2)
abline(v = 38.3, col = "black", lwd = 1, lty = 2)
abline(v = 41.3, col = "black", lwd = 1, lty = 2)
abline(v = 44.3, col = "black", lwd = 1, lty = 2)
axis(1, seq(0,n+1,2))
#axis(2, seq(0,5000,500))
Fit3 <- nls(y ~ exp(c0+c1/x),
start=list(c0=1, c1=20), trace=FALSE)
summary(Fit3)
##
## Formula: y ~ exp(c0 + c1/x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## c0 14.5515 0.0307 474.0 <2e-16 ***
## c1 -301.0527 2.5379 -118.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 962.8 on 91 degrees of freedom
##
## Number of iterations to convergence: 22
## Achieved convergence tolerance: 8.396e-07
#confint(Fit3)
# Add in predicted values into the plot
lines(x, fitted.values(Fit3), lwd=3, col="black")
#AIC(Fit3)
#BIC(Fit3)
#Distribucion de residuales
#plot(fitted(Fit3), residuals(Fit3), xlab="Valores estimados", ylab = "Residuals",
# abline(a=0, b=0), col="lightskyblue", pch=19, cex=2)
#Analisis de autocorrelacion de los residuales
#acf(residuales_Fit3<-(y-fitted(Fit3)))
#Osbervados vs estimados
#plot(fitted(Fit3), y, abline(a=0, b=1), xlab="Valores estimados", ylab ="Valores observados", col="lightskyblue", pch=19, cex=2)
#Los coeficientes de regresion
coef(summary(Fit3))
## Estimate Std. Error t value Pr(>|t|)
## c0 14.55154 0.03070161 473.9666 3.602863e-156
## c1 -301.05265 2.53785720 -118.6247 1.519336e-101
df3<-as.data.frame(coef(summary(Fit3))); df3
## Estimate Std. Error t value Pr(>|t|)
## c0 14.55154 0.03070161 473.9666 3.602863e-156
## c1 -301.05265 2.53785720 -118.6247 1.519336e-101
#Predicciones cuando x sea
x2<-n+1
df3[1, 1]
## [1] 14.55154
df3[2, 1]
## [1] -301.0527
#Numero de posibles infectados estimados por el modelo tres.
Prediccion_m3<- exp(df3[1, 1]+df3[2, 1]/x2)
b<-df3[1, 1]
Eq1<-round(Prediccion_m3)
this_day <- today()+1
hoy <- Sys.Date();
mnn<-(hoy+1);
mnnn<-format(mnn, format="%d %B %Y")
text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=2)
text(7, (Eq1-5350), labels = Eq1, col = "black", cex=1.8)
text(17, (Eq1-11800), "Posibles infectados para el dia", col = "black")
text(45, (Eq1-11800), labels = mnnn, col = "black")
#arrows(0, Eq1-100, 31.3, Eq1-100, col = "red")
arrows(0, 100, 31.3, 100, col = "#CD1076")
arrows(31.3, 100, 38.3, 100, col = "#CD1076")
arrows(38.3, 100, 41.3, 100, col = "#CD1076")
arrows(38.3, 100, 44.3, 100, col = "#CD1076")
text(34, 180, "mil", col = "#CD1076", cex=0.7)
text(27, 180, "mil", col = "#CD1076", cex=0.7)
text(39, 180, "mil", col = "#CD1076", cex=0.7)
#Bondad de ajuste del modelo
#La R2
df.residual(Fit3)
## [1] 91
CMErr_Fit3<-sum((y-fitted(Fit3))^2)/(df.residual(Fit3))
CMTot_Fit3<-(sum((y-mean(y))^2))/((length(y)-1))
R2_Fit3<-1-(CMErr_Fit3)/(CMTot_Fit3)
#El Error
Sxy_Fit3<-sqrt(CMErr_Fit3)
#El Coeficiente de variacion
CV_Fit3<-(Sxy_Fit3)/(mean(y))*100
I<--((df3[2, 1])*(exp(((df3[2, 1])/(x))+(df3[1, 1]))))/(x^2)
df_I <- data.frame(x, I)
df_II <-round(df_I,0)
plot(x, I, xlab="Dia numero", ylab="Incremento medio infectados por dia", col="lightskyblue", pch=19, cex=2)
lines(x, I, lwd=3, col="black")
text(I~x, labels=I, data=df_II, cex=0.7, font=1, pos=2)
axis(1, seq(0,n+1,2))
#abline(v = 0, col = "black", lwd = 1, lty = 2)
abline(v = 31.3, col = "black", lwd = 1, lty = 2)
abline(v = 38.3, col = "black", lwd = 1, lty = 2)
## Intervalos de confianza (95 %) Simulacion Monte Carlo modelo 3
CONC_Fit3 <- seq(n+1, n+20, by = 1)
Pred_Fit3 <- predict(Fit3, newdata = data.frame(x = CONC_Fit3), nsim = 100000)
IC_95_Fit3 <- predictNLS(Fit3, newdata = data.frame(x = CONC_Fit3), nsim = 100000);
## predictNLS: Propagating predictor value #1...
## predictNLS: Propagating predictor value #2...
## predictNLS: Propagating predictor value #3...
## predictNLS: Propagating predictor value #4...
## predictNLS: Propagating predictor value #5...
## predictNLS: Propagating predictor value #6...
## predictNLS: Propagating predictor value #7...
## predictNLS: Propagating predictor value #8...
## predictNLS: Propagating predictor value #9...
## predictNLS: Propagating predictor value #10...
## predictNLS: Propagating predictor value #11...
## predictNLS: Propagating predictor value #12...
## predictNLS: Propagating predictor value #13...
## predictNLS: Propagating predictor value #14...
## predictNLS: Propagating predictor value #15...
## predictNLS: Propagating predictor value #16...
## predictNLS: Propagating predictor value #17...
## predictNLS: Propagating predictor value #18...
## predictNLS: Propagating predictor value #19...
## predictNLS: Propagating predictor value #20...
#IC_95_Fit3$summary
H<-data.frame(IC_95_Fit3$summary)
#Extrayendi
IC_inf_m3<-IC_95_Fit3$summary[1, 5]
IC_sup_m3<-IC_95_Fit3$summary[1, 6]
Dia<-seq(Sys.Date()+1, length.out =20, by="1 day")
#data.frame(fechas)
#select(IC_95_Fit3$summary, Prop.Mean.1)
#H[, 5:6, drop=FALSE]
Prediccion<-(data.frame(H[, c(1,5,6),] ))
Pred_li<-df[Prediccion[1,5]]
#cbind(Dia, Prediccion)
as.datatable(formattable(cbind(Dia, Prediccion)))
#formattable(cbind(Dia, Prediccion))
El modelo: \(y=\dfrac{a\mathrm{e}^{cx}}{\mathrm{e}^{b-x}+1}\) Su primera derivada: \(\frac{dy}{dx}= \dfrac{a\mathrm{e}^\frac{b-\ln\left(x\right)}{c}}{cx\left(\mathrm{e}^\frac{b-\ln\left(x\right)}{c}+1\right)^2}\) Simplificada se reescribe: \(\frac{dy}{dx}=\dfrac{a\left(c\mathrm{e}^x+\mathrm{e}^bc+\mathrm{e}^b\right)\mathrm{e}^{cx+x}}{\left(\mathrm{e}^x+\mathrm{e}^b\right)^2}\)
# Modelo 4 ----------------------------------------------------------
plot(x, y, col="darkolivegreen1", pch=19, xlab = "Dia numero", ylab = "Numero de infectados", cex=2, main="Prediccion del numero de infectados con el modelo cuatro")
abline(v = 31.3, col = "black", lwd = 1, lty = 2)
abline(v = 38.3, col = "black", lwd = 1, lty = 2)
Fit4 <- nlsLM(y ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x),
start=list(p1=410, p2=20, p4=.03), trace=FALSE)
## Warning in nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower, : lmdif: info = -1. Number of iterations has reached `maxiter' == 50.
#text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=1)
#text(15, 1500, "Estimados para el dia 05 de abril del 2020 = 2227", col = "blue")
summary(Fit4)
##
## Formula: y ~ p1/(1 + exp(p2 - x)) * exp(p4 * x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## p1 6.982e+02 3.813e+01 18.31 <2e-16 ***
## p2 4.359e+01 5.747e-01 75.85 <2e-16 ***
## p4 5.208e-02 6.448e-04 80.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1573 on 90 degrees of freedom
##
## Number of iterations till stop: 50
## Achieved convergence tolerance: 1.49e-08
## Reason stopped: Number of iterations has reached `maxiter' == 50.
#confint(Fit3)
# Add in predicted values into the plot
lines(x, fitted.values(Fit4), lwd=3, col="blue")
#AIC(Fit4)
#BIC(Fit4)
#Distribucion de residuales
#plot(fitted(Fit4), residuals(Fit4), xlab="Valores estimados", ylab = "Residuals",
#abline(a=0, b=0), col="blue", pch=16)
#Analisis de autocorrelacion de los residuales
#acf(residuales_Fit4<-(y-fitted(Fit4)))
#Osbervados vs estimados
#plot(fitted(Fit4), y, abline(a=0, b=1), col="black", xlab="Valores estimados", ylab ="Valores observados", pch=16)
#Los coeficientes de regresion
coef(summary(Fit4))
## Estimate Std. Error t value Pr(>|t|)
## p1 698.17419845 3.812513e+01 18.31270 4.178967e-32
## p2 43.59344101 5.747071e-01 75.85332 2.326912e-83
## p4 0.05207766 6.448242e-04 80.76257 8.927059e-86
df4<-as.data.frame(coef(summary(Fit4)))
#Predicciones cuando x sea
x3<-n+1
df4[1, 1]
## [1] 698.1742
df4[2, 1]
## [1] 43.59344
df4[3, 1]
## [1] 0.05207766
#Numero de posibles infectados estimados por el modelo
Prediccion_m4<- df4[1, 1] / (1 + exp(df4[2, 1] - x3)) * exp(df4[3, 1] * x3)
Eq1<-round(Prediccion_m4)
this_day <- today()+1
hoy <- Sys.Date();
mnn<-(hoy+1);
mnnn<-format(mnn, format="%d %B %Y")
text(y~x, labels=y, data=datos, cex=0.7, font=1, pos=2)
text(7, (Eq1-5850), labels = Eq1, col = "black", cex=1.8)
text(17, (Eq1-11300), "Posibles infectados para el dia", col = "black")
text(45, (Eq1-11300), labels = mnnn, col = "black")
#Bondad de ajuste del modelo
#La R2
df.residual(Fit4)
## [1] 90
CMErr_Fit4<-sum((y-fitted(Fit4))^2)/(df.residual(Fit4))
CMTot_Fit4<-(sum((y-mean(y))^2))/((length(y)-1))
R2_Fit4<-1-(CMErr_Fit4)/(CMTot_Fit4)
#El Error
Sxy_Fit4<-sqrt(CMErr_Fit4)
#El Coeficiente de variacion
CV_Fit4<-(Sxy_Fit4)/(mean(y))*100
## Intervalos de confianza (95 %) Simulacio Monte Carlo modelo 3
CONC_Fit4 <- seq(n+1, n+1, by = 1)
Pred_Fit4 <- predict(Fit4, newdata = data.frame(x = CONC_Fit4), nsim = 100000)
IC_95_Fit4 <- predictNLS(Fit4, newdata = data.frame(x = CONC_Fit4), nsim = 100000)
## predictNLS: Propagating predictor value #1...
IC_95_Fit4$summary
## Prop.Mean.1 Prop.Mean.2 Prop.sd.1 Prop.sd.2 Prop.2.5% Prop.97.5% Sim.Mean
## 1 93318.04 93182.22 786.6016 810.3336 91572.35 94792.09 93182.17
## Sim.sd Sim.Median Sim.MAD Sim.2.5% Sim.97.5%
## 1 810.9868 93235.69 775.9953 91441.33 94630.27
IC_inf_m4<-IC_95_Fit4$summary[1, 5]
IC_sup_m4<-IC_95_Fit4$summary[1, 6]
media<-(round(Prediccion_m1)+round(Prediccion_m2)+round(Prediccion_m3)++round(Prediccion_m4))/4
Curva=read.csv("G:/UAAAN/MATERIAS/2019/REGRESION APLICADA/NLS/Curva.csv")
attach(Curva)
head(Curva)
## Dia_z z
## 1 27/02/2020 1
## 2 28/02/2020 2
## 3 29/02/2020 3
## 4 01/03/2020 4
## 5 02/03/2020 5
## 6 03/03/2020 6
Prediccion_z<- exp(df3[1, 1]+df3[2, 1]/z)
plot(z, Prediccion_z, main="Veamos primero la prediccion desde el dia uno hasta\n 200 dias (13-sep-2020) usando el modeo tres", xlab="Dias desde el primer contagio", ylab="Numero de infectados", col="lightskyblue", pch=19)
abline(v = n, col = "black", lwd = 1, lty = 2)
text(n, (55000), labels = "hoy", col = "black", cex=1.2) #
text(n, (30000), labels = hoy, col = "#CD1076", cex=1.8)
Incremento_z<--((df3[2, 1])*(exp(((df3[2, 1])/(z))+(df3[1, 1]))))/(z^2)
plot(z,Incremento_z, main="Fecha en que se predice se alcance el pico maximo\n de infectados", xlab="Dias desde el primer contagio", ylab="Incremento medio de infectados por dia", col="lightskyblue", pch=19)
abline(v = n, col = "black", lwd = 1, lty = 2)
text(n, 50, labels = hoy, col = "#CD1076", cex=1.8)
max_z<-round(max(Incremento_z))
text(62, (40), "hoy", col = "black", cex=1.2)
arrows(144, 0, 144, 3213, col = "#CD1076")
n2<- 149 #N Fecha del maximo ICA
hoy <- Sys.Date()
n2-n #Dias que faltan
## [1] 56
#arrows(n, 300, n2, 300, col = "lightskyblue")
#"Dias que faltan "
text(n/2+75, 1500, "Faltan ... días", col = "black", cex=1)
text(n/2+80, 1500, labels = n2-n, col = "#1E90FF", cex=1.8)
Que_dia<-(hoy+days(n2-n))
aaa<-format(Que_dia, format="%d %B %Y")
text(120, 2500, labels = aaa, col = "#CD1076", cex=1.8)
#el_dia<-(hoy+(labels=max_z)); el_dia
El modelo \(y=a\mathrm{e}^{-b\mathrm{e}^{-kx}}\) su derivada \(dy/dx=abk\mathrm{e}^{-b\mathrm{e}^{-kx}-kx}\)
datos_2=read.csv("Coah_29-05-20.csv")
attach(datos_2)
head(datos_2)
## Fecha t N y2
## 1 01/04/2020 1 64 4
## 2 02/04/2020 2 68 3
## 3 03/04/2020 3 71 7
## 4 04/04/2020 4 78 10
## 5 05/04/2020 5 88 13
## 6 06/04/2020 6 101 12
n2=length(t)
plot(t, N, col="#98FB98", pch=19, xlab = "Dia numero", ylab = "Numero de infectados", cex=1.5, main="Prediccion del numero de infectados Coahuila")
Fit5 <- nlsLM(N ~ g1*exp(-g2*exp(-g3*t)),
start=list(g1=410, g2=20, g3=.03), trace=FALSE)
text(N~t, labels=N, data=datos_2, cex=0.7, font=1, pos=2)
#text(15, 1500, "Estimados para el dia 05 de abril del 2020 = 2041", col = "red")
summary(Fit5)
##
## Formula: N ~ g1 * exp(-g2 * exp(-g3 * t))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## g1 5.110e+03 8.374e+02 6.102 1.04e-07 ***
## g2 4.117e+00 1.265e-01 32.537 < 2e-16 ***
## g3 1.669e-02 1.250e-03 13.352 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.6 on 56 degrees of freedom
##
## Number of iterations to convergence: 25
## Achieved convergence tolerance: 1.49e-08
# Add in predicted values into the plot
lines(t, fitted.values(Fit5), lwd=3, col="#FF4500")
#Los coeficientes de regresion
coef(summary(Fit5))
## Estimate Std. Error t value Pr(>|t|)
## g1 5.109710e+03 837.4455466 6.101543 1.038644e-07
## g2 4.117110e+00 0.1265348 32.537364 4.631939e-38
## g3 1.669434e-02 0.0012503 13.352271 4.779878e-19
df5<-as.data.frame(coef(summary(Fit5)))
#Predicciones cuando x sea
x4<-n2+1
df5[1, 1]
## [1] 5109.71
df5[2, 1]
## [1] 4.11711
df5[3, 1]
## [1] 0.01669434
#Numero de posibles infectados estimados por el modelo para el dia
Prediccion_m5<-df5[1, 1]*exp(-df5[2, 1]*exp(-df5[3, 1]*x4))
zz<-1:100
Prediccion_zz<- df5[1, 1]*exp(-df5[2, 1]*exp(-df5[3, 1]*zz))
plot(zz, Prediccion_zz, main="Veamos primero la prediccion desde el dia uno (01 de abril del 2020) hasta\n 100 dias (09-jul-2020)", xlab="Dias desde el primer contagio", ylab="Numero de infectados", col="#98FB98", pch=19)
abline(v = n2, col = "black", lwd = 1, lty = 2)
text(n2, (500), labels = "hoy", col = "black", cex=1.2) #
text(n2, (1000), labels = hoy, col = "#CD1076", cex=1.8)
Incremento_zz<-((df5[1, 1])*(df5[2, 1])*(df5[3, 1]))*exp(((-df5[2, 1])*exp(-(df5[3, 1])*zz))-(df5[3, 1])*(zz))
plot(zz,Incremento_zz, main="Fecha en que se predice se alcance el pico maximo\n de infectados", xlab="Dias desde el primer contagio", ylab="Incremento medio de infectados por dia", col="#98FB98", pch=19)
#
abline(v = n2, col = "black", lwd = 1, lty = 2) #
text(n2, 15, labels = hoy, col = "#CD1076", cex=1.8) #
max_zz<-round(max(Incremento_zz)) #
text(n2, (20), "hoy", col = "black", cex=1.2) #
arrows(84, 0, 84, 30, col = "#CD1076") #
n3<- 84 #N Fecha del maximo ICA #
hoy <- Sys.Date() #
n3-n2 #Dias que faltan
## [1] 25
#arrows(n, 300, n2, 300, col = "lightskyblue")
#"Dias que faltan "
text(n2/2, 30, "Faltan ... días", col = "black", cex=1)
text(n2/2+2, 30, labels = n3-n2, col = "#98FB98", cex=1.8)
Que_dia<-(hoy+days(n3-n2))
aaa<-format(Que_dia, format="%d %B %Y")
text(85, 25, labels = aaa, col = "#98FB98", cex=1.8)
Modelo | I.C. inferior (x 1000) | Prediccion (x 1000) | I.C. superior (x 1000) |
---|---|---|---|
Modelo uno | 87.703 | ||
Modelo dos | 85.668 | ||
Modelo tres | 84.042 | 84.868 | 85.696 |
Modelo cuatro | 91.572 | 93.318 | 94.792 |
En promedio se esperan 87.889 infectados para el dia 30 mayo 2020.
El modelo tres predice que para el dia 05 de junio podriamos tener 100, 000 infectados (ver la Simulacion Monte Carlo de I.C. 95 %).
IMPORTANTE. El pico maximo de infectados se esta prolongando a como se habia estimado inicialmente