[R-es] Pronósticos con modelos robusto de series de tiempo
Jose Betancourt B.
betanster en gmail.com
Vie Sep 6 00:29:27 CEST 2013
con el paquete EpiEstim
por favor ?como exactamente calculan las probabilidades para conformar
el segundo vector, aqui yo puse uno de SAR2003
#Análisis paso a paso
rm(list = ls())
#existe información de la incidencia de casos de una enfermedad
incidencia <-c (1,0,0,0,0,0,0,0,2,1,1,0,1,2,0,2,0,2,2,2,3, 1, 2, 0,2,
0, 1, 1, 2, 1, 1, 2, 3, 0, 1, 2, 0, 2, 1, 3, 1, 6, 3, 3,7, 4,
4, 3, 6, 4,5,2,4, 2,3, 1, 4, 4,8, 9, 13, 1, 0, 0, 0, 0, 0,
0,0, 0, 0,0, 0, 0, 0,0, 0, 0,
0, 0, 0, 28, 35, 8,
59, 13, 14, 1, 26, 15, 1,
1, 9, 41, 18, 56, 104,
54, 59, 50, 20, 2, 1, 1,
42, 20, 21, 24, 37, 34,
26, 50, 30, 17, 28, 14,
75, 62, 37, 19, 12, 11,
17, 27, 0, 34, 12, 0,
15, 0, 3, 8, 14, 11, 10,
11, 9, 5, 31, 11, 3,
2, 3, 0, 0, 0, 26,
0, 0, 0, 0, 16, 0,
8, 3, 2, 0, 0, 25, 2,
14, 5, 5, 2, 0, 5, 0,
2, 2, 3, 2, 0,
8, 2, 5, 3, 2, 1,
1)
#se explora la información
library (epicalc)
summ (incidencia)
hist(incidencia)
## se calcula el intervalo serial
MeanFluSI <- 9.418994 #puse los datos de la incidencia
SdFluSI <- 16.06
DicreteSIDistr <- vector()
for(i in 0:35)
{
DicreteSIDistr[i+1] <- DiscrSI(i, MeanFluSI, SdFluSI)
}
plot(0:35, DicreteSIDistr, type="h", lwd=10, lend=1, xlab="tiempo
(dias)", ylab="frequencia")
title(main="distribution Discreta del intervalo serial ")
# voy a poner unas probabilidades SARS 2003; no se como calcularlas
para este caso a partir de la incidencia y esa es mi duda en este
trabajo
SI.Distr<- c (0.000, 0.001, 0.012, 0.043, 0.078, 0.104, 0.117, 0.116,
0.108, 0.094, 0.078, 0.063, 0.049,
0.038, 0.028, 0.021, 0.015, 0.011, 0.008, 0.005, 0.004, 0.003,
0.002, 0.001, 0.001)
#Calculando R
EstimateR(incidencia, T.Start=90:140, T.End=95:145,
method="NonParametricSI",SI.Distr=SI.Distr, plot=TRUE,
leg.pos=xy.coords(1,7))
#uniendo os dos vectores
enfermedad <- list(incidencia, SI.Distr, .Names = c("Incidence",
"SI.Distr"))#como combinar los dos vectores
enfermedad<- list(incidence=a, SI.Distr=b)
enfermedad
#metodo Wallinga...
WT(incidencia, T.Start=90:140, T.End=95:145, method="NonParametricSI",
SI.Distr=SI.Distr, plot=TRUE, leg.pos=xy.coords(1,1.75), nSim=100)
WT(incidencia, T.Start=90:140, T.End=95:145, method="ParametricSI",
Mean.SI=9.418994, Std.SI=16.06, plot=TRUE, nSim=100)
#calcular infectividad general
lambda <- OverallInfectivity(incidencia, SI.Distr)
par(mfrow=c(2,1))
plot(incidence, type="s", xlab="tiempo (dias)", ylab="Incidencia")
title(main="Curva epidemica ")
plot(lambda, type="s", xlab="tiempo (dias)", ylab="Infectividad")
title(main="Infectividad general")
lambda
2013/9/5, Jose Betancourt B. <betanster en gmail.com>:
> este ha sido bastante bueno
> f52<-scan("f52.txt")
> sb<-f52[1:200]; sp<-f52[201:245];
> sb.ts <- ts(sb, start=c(2008, 1), frequency=52)
> sb.ts # is important to see the end of this time series.
> sp.ts<- ts(sp, start=c(2011, 45), frequency=52)
> plot(sb.ts,ylab="No. focos")
>
> # Main function
> predAll<-function(x,h,tsname){
> # INPUTS:
> # x: time serie
> # h: number of values to predict.
> # tsname: maximun number of element in the hidden layer (nhmax=20).
> #
> ##---- PACKAGES -----
> library(forecast)
> ##---- MODELS -------
> m1<-stl(x, s.window='periodic') #STL Decomposition
> m2<-auto.arima(x) # arima
> m3<-nnetar(x,lambda=0) #neural networks
> # Predictions
> f1<-forecast(m1, h=h);
> f2<-forecast(m2, h=h);
> f3<-forecast(m3, h=h);
> # Fit accuracy
> rmse1 <- mean((x-f1$fitted)^2)/var(x);
> rmse2 <- mean((x-f2$fitted)^2)/var(x);
> fitted<-f3$fitted;#fitted[1:52]<-x[1:52];
> rmse3 <-
>
> mean((x[(frequency(x)+1):length(x)]-fitted[(frequency(x)+1):length(x)])^2)/var(x[(frequency(x)+1):length(x)]);
> #lab<-data.frame("F","R2","R2a","RSS")
> mat<-cbind(c(rmse1),c(rmse2),c(rmse3))
> tab<-as.data.frame(mat)
> dimnames(tab)[[1]]<-c("RMSE")
> dimnames(tab)[[2]] <-c("STL","ARIMA","NNAR")
> print(tab)
> # Plot all predictions together
> par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
> plot(f1,xlab="Time in Weeks",ylab=tsname,main="All models")
> lines(f2$mean,col=1,lwd=2)
> lines(f3$mean,col=3,lwd=2)
> lines(x,col=2,lwd=2)
> legend.txt = c("STL","ARIMA","NNAR","Observed");
> legend("topright",inset=c(-0.32,0),legend.txt,col = c(4,1,3,2),
> lwd=c(2,2,2,2),lty = c(1,1,1,1,1), pch =
> c(NA,NA,NA,NA),title="Models",bty="n");
> return(list(stl=f1, arima=f2, nnar=f3))
> }
> result<-predAll(sb.ts,45,"Foci Number")
> lines(sp.ts,col=2,lwd=2)# to test these predictions
>
>
> ve al paquete AER
>
>
>
>
> El 05/09/13, Enrique RAMOS <ceramos0 en yahoo.com.mx> escribió:
>> Alguien me podría sugerir un paquete en R para generar pronóticos con
>> modelos robusto de series de tiempo.
>>
>> Saludos
>> Enrique RAMOS
>> [[alternative HTML version deleted]]
>>
>>
>
Más información sobre la lista de distribución R-help-es