File:Covid cases forecast from r0.svg
Original file (SVG file, nominally 900 × 450 pixels, file size: 90 KB)
Captions
Summary
[edit]DescriptionCovid cases forecast from r0.svg |
English: Covid cases from r0 in Finland. |
Date | |
Source | Own work |
Author | Merikanto |
- Calculate forecast of Covid-19 cases from r0 in Finland
- 31.8.2021
- v 0000.0000
new_in_skript=0
if (new_in_skript==1)
{
#install.packages("ggplot2", "plotly", repos ="https://ftp.acc.umu.se/mirror/CRAN/")
install.packages("svglite")
install.packages("ggplot2")
install.packages("rvest")
install.packages("readtext")
install.packages("stringi")
install.packages("datamart")
install.packages("XML")
install.packages("tidyr")
install.packages("stringr")
install.packages("stringi")
install.packages("tibble")
install.packages("readr")
install.packages("data.table")
install.packages("caTools")
install.packages("mgcv")
install.packages("repmis")
install.packages("lubridate")
install.packages("tidyverse")
install.packages("R0")
install.packages("EpiEstim")
install.packages("jsonlite")
install.packages("rjstat")
install.packages("dplyr")
}
library(ggplot2)
library(svglite)
library(rvest)
library(readtext)
library(stringi)
library(stringr)
library(datamart)
library(XML)
library(jsonlite)
library(tibble)
library(caTools)
library(mgcv)
library(repmis)
library(lubridate)
library(tidyverse)
library(dplyr)
library(tidyr)
library(readr)
library(data.table)
library(rjstat)
library (R0)
library(EpiEstim)
library(forecast)
library(prophet)
library(astsa)
library(MLmetrics)
library(tseries)
- choices
- 1 finnish wiki data, 2 aggregated cases data
- 3 solanpaa finnish data 4 thl cube json data
- 5 thl cube josn hospital data
load_data_from=3
- beginday1='28/02/2020'
- beginday1='1/1/2021'
beginday1='1/8/2020'
forecastendday1<-"2021/10/16"
- today=Sys.Date()-1
today=Sys.Date()-3
spanni=0.3
- yala=0.6
- yyla=3.2
yala=0.7
yyla=1.4
metodi="loess"
widthi=10
heighti=5
plottaa=1 ## must be 1
tulosta_svg=1 # plot to out svg 0, 1 of 2
tulosfilee1="./R0_Suomessa_2.svg"
beginday0=as.Date(beginday1)
beginday2=format(beginday0, "%Y/%m/%d")
today1=format(today, "%d/%m/%Y")
today2=format(today, "%Y/%m/%d")
print(today1)
datelimits1=c(beginday1, today1)
paivat1=seq(as.Date(beginday2), as.Date(today2), "days")
forecast_profet<-function(xy2, futuredays)
{
## calculate r0 w/r0 package
ds<-as.Date(xy2$Dates)
y<-as.integer(xy2$Cases)
nummeros<-1:length(ds)
lenu=length(ds)
df<-data.frame(ds,y)
m <- prophet(df)
future <- make_future_dataframe(m, periods = futuredays)
forecast <- predict(m, future)
#str(future)
#str(forecast)
futu_days=future$ds
futu_trendi=forecast$trend
futu_trendi_upper=forecast$trend_upper
futu_trendi_lower=forecast$trend_lower
futu_yhat=forecast$yhat
futu_yhat_upper=forecast$yhat_upper
futu_yhat_lower=forecast$yhat_lower
futu_weekly=forecast$weekly
futu_weekly_upper=forecast$weekly_upper
futu_weekly_lower=forecast$weekly_lower
xypaluu<-data.frame(as.Date(futu_days),futu_yhat)
# xypaluu<-data.frame(as.Date(futu_days),futu_weekly)
#plot(r0t1)
#print (r0t1)
#stop(-1)
names(xypaluu)<-c("paivat","r0")
return(xypaluu)
}
calculate_r0 <- function(time1, time2, val1, val2)
{
td=time2-time1
gr0<-log(val2/val1)
gr=gr0/td
td = log(2)/gr
tau<-5.0
k<-log(2.0)/td
r0<-exp(k*tau)
return(r0)
}
moving_average <- function(x, w, FUN, ...)
{
if (w < 1) {
stop("Window length: mustbe greater than 0")
}
output <- x
for (i in 1:length(x)) {
lower_bound <- i - w + 1
if (lower_bound < 1) {
output[i] <- NA_real_
## !!! assume NA 0
output[i] <- 0
} else {
output[i] <- FUN(x[lower_bound:i, ...])
}
}
return (output)
}
calculate_multiple_r0 <- function(daata1) {
lenu1<-length(daata1)
daata2<-1:lenu1
for (n in 2:lenu1){
valju1=daata1[n-1]
valju2=daata1[n]
timex1=0
timex2=1
r0<-calculate_r0(0, 1, valju1, valju2)
daata2[n]<-r0
#print (r0)
}
return(daata2)
}
load_data_from_finnish_wiki<-function()
{
url1="https://fi.wikipedia.org/wiki/Suomen_koronaviruspandemian_aikajana"
destfile1="./ward0.txt"
download.file(url1, destfile1)
texti000<-readtext(destfile1)
texti0<-texti000$text
etsittava1="1. huhtikuuta 2020 alkaen"
len1=nchar(texti0)
k1=regexpr(pattern=etsittava1, texti0)
k1b=len1-k1
texti1=strtail(texti0,k1b)
sink("out1.txt")
print (texti1)
sink()
etsittava2=""
k2=regexpr(pattern=etsittava2, texti1)
texti2=strhead(texti1,k2)
sample1<-minimal_html(texti2)
tabu1 <- html_table(sample1, fill=TRUE)1
colnames(tabu1) <- c("V1","V2", "V3","V4", "V5","V6", "V7","V8" )
- print(tabu1)
sairaalassa00<-tabu1$V4
sairaalassa=as.integer(sairaalassa00)
teholla00<-tabu1$V5
teholla=as.integer(teholla00)
uusiatapauksia00<-tabu1$V3
uusiatapauksia0<-gsub(" ", "", uusiatapauksia00)
uusia_tapauksia=as.integer(uusiatapauksia0)
uusiakuolleita00<-tabu1$V7
uusiakuolleita1=as.integer(uusiakuolleita00)
uusiakuolleita2<-uusiakuolleita1
uusiakuolleita2[uusiakuolleita2<0]<-0
uusia_kuolleita<-uusiakuolleita2
toipuneita00<-tabu1$V8
toipuneita01<-gsub(" ", "", toipuneita00)
toipuneita0<-gsub("[^0-9.-]", "", toipuneita01)
toipuneita=as.integer(toipuneita0)
tapauksia00<-tabu1$V2
tapauksia01<-gsub(" ", "", tapauksia00)
tapauksia0<-gsub("[^0-9.-]", "", tapauksia01)
tapauksia=as.integer(tapauksia0)
kuolleita00<-tabu1$V6
kuolleita=as.integer(kuolleita00)
pv0<-tabu1$V1
len1=length(pv0)
daates1 <- vector(mode="character", length=len1)
- print(pv0)
n=1
for(n in 1:len1)
{
it1<-pv0[n]
#print(it1)
qq1<-str_split(it1, "\\[")1
qq2<-qq1[1]
qq3<-gsub(" ", "", qq2, fixed = TRUE)
daates1[n]=qq3
}
daates2=as.Date(daates1, format="%d.%m.%Y")
print(daates2)
aktiivisia_tapauksia=tapauksia-kuolleita-toipuneita
- print (paivat1)
- print (teholla)
- print (sairaalassa)
- print (tapauksia)
- print (kuolleita)
- print (toipuneita)
- print (uusia_tapauksia)
- print (uusia_kuolleita)
- plot(paivat1,aktiivisia_tapauksia)
- xy<-data.frame(daates2, sairaalassa)
xy<-data.frame(daates2, uusia_tapauksia)
names(xy)<-c("Dates", "Cases")
xyz<-data.frame(daates2, sairaalassa, teholla)
dfout1<-data.frame(daates2, aktiivisia_tapauksia, uusia_tapauksia, sairaalassa, teholla, uusia_kuolleita )
names(dfout1)<-c("Pvm", "Aktiivisia_tapauksia","Uusia_tapauksia", "Sairaalassa", "Teholla", "Uusia_kuolleita")
write.csv2(dfout1, "./sairaalassa.csv",row.names=FALSE )
return(xy)
}
load_data_from_aggregated<-function()
{
- fetch the data
srkurl='https://datahub.io/core/covid-19/r/countries-aggregated.csv'
dfine000 <- read.csv(file=srkurl)
dfine<-as.data.frame(dfine000)
- head(dfine)
- class(dfine)
- print(dfine)
- tail(dfine)
- stop(-1)
dfinland <- dfine[ which(dfine$Country=='Finland'), ]
- head(dfinland)
- print(dfinland)
- quit(-2)
kols <- c("Date", "Confirmed","Recovered","Deaths")
tapaukset <- dfinland[kols]
- head(tapaukset)
len1=nrow(tapaukset)
- len1
len2=len1-1
len3=len2
confirmed<-tapaukset$Confirmed
deaths<-tapaukset$Deaths
dailycases <- vector()
dailycases <- c(dailycases, 0:(len2))
dailydeaths <- vector()
dailydeaths <- c(dailydeaths, 0:(len2))
m=0
dailycases[1]<-tapaukset$Confirmed[1]
- dailydeaths[1]<-tapaukset$Deaths[1]
dailydeaths[1]<-0
- confirmed
- deaths
m=1
for(n in 2:(len3+1)) {
a<-confirmed[n]
b<-confirmed[m]
#print (a)
#print (b)
cee<- (a-b)
#print(cee)
dailycases[n]=cee
m=m+1
}
mm=1
for(nn in 2:(len3+1)) {
aa<-deaths[nn]
bb<-deaths[mm]
#print ("_")
#print (aa)
#print (bb)
ceb=aa-bb
#if (ceb<0) ceb=0
#print(ceb)
dailydeaths[nn]=ceb
mm=mm+1
}
- deaths
- dailycases
- dailydeaths
dfout1<-dfinland
- print(nrow(dfinland))
- print(length(dailydeaths))
dfout1 <- cbind(dfout1, data.frame(dailycases))
dfout1 <- cbind(dfout1, data.frame(dailydeaths))
- head(dfout1)
dfout2<-within(dfout1, rm(Country))
names(dfout2) <- c('Date','Confirmed','Recovered','Deaths', 'DailyConfirmed','DailyDeaths')
- head(dfout2)
write.csv2(dfout2, "/Users/himot/akor1/finland_data1.csv");
daate1<-dfout2$Date
dailydeaths1<-dfout2$DailyDeaths
dailycases1<-dailycases
- daate1
- daate2<-gsub("2020-", "", daate1)
daate2<-daate1
leenu<-length(daate2)
- alkupvm<-50
alkupvm<-1
daate3<-daate2[alkupvm:leenu]
dailydeaths3<-dailydeaths1[alkupvm:leenu]
dailycases3<-dailycases1[alkupvm:leenu]
- daate3
- dailydeaths3
pv0<-tabu1$V1
len1=length(pv0)
daates1 <- vector(mode="character", length=len1)
- print(pv0)
n=1
for(n in 1:len1)
{
it1<-pv0[n]
#print(it1)
qq1<-str_split(it1, "\\[")1
qq2<-qq1[1]
qq3<-gsub(" ", "", qq2, fixed = TRUE)
daates1[n]=qq3
}
daates2=as.Date(daates1, format="%d.%m.%Y")
print(daates2)
# barplot(dailydeaths3, main="Koronaviruskuolemat päivittäin vuonna 2020",
# names.arg=daate3)
dataf1 <- data.frame("Date" = daates2, "Paivitt_kuolemat"=dailydeaths3)
- str(dataf1)
dataf2 <- data.frame("Date" = daates2, "Paivitt_tapaukset"=dailycases3)
- str(dataf2)
write.csv(dataf1, "/Users/himot/akor1/dailydeaths1.csv", row.names=T)
write.csv(dataf2, "/Users/himot/akor1/dailycases1.csv", row.names=T)
indf1 <- read.csv(file = '/Users/himot/akor1/dailycases1.csv')
#head(indf1)
cases1<-indf1$Paivitt_tapaukset
dates1<-indf1$Date
len1=length(cases1)
dates2<-as.Date(dates1)
paivat<-1:len1
xy<-data.frame(daates2, dailycases3)
return(xy)
}
download_solanpaa_finnish_data<-function()
{
solanpaa_fi="https://covid19.solanpaa.fi/data/fin_cases.json"
cache_file="solanpaa_fi.json"
download.file(solanpaa_fi, cache_file)
j1 <- fromJSON(cache_file)
## maybe errori
dates<-as.Date(j1$date)
dailycases<-j1$new_cases
dailydeaths<-j1$new_deaths
dataf1 <- data.frame("Date" = dates, "Paivitt_kuolemat"=dailydeaths)
dataf2 <- data.frame("Date" = dates, "Paivitt_tapaukset"=dailycases)
write.csv(dataf1, "./dailydeaths1.csv", row.names=T)
write.csv(dataf2, "./dailycases1.csv", row.names=T)
xy0<-data.frame(dates, dailycases)
names(xy0)<-c("Dates", "Cases")
xy<-na.omit(xy0)
return(xy)
}
calculate_r0_with_r0<-function(xy2)
{
## calculate r0 w/r0 package
dates<-as.Date(xy2$Dates)
cases<-as.integer(xy2$Cases)
cases[is.na(cases)] <- 1
cases[(cases<0)] <- cases*-1
cases[cases==0] <- 1
nummeros<-1:length(dates)
num<-cases
#names<-nummeros
names<-dates
lenu=length(dates)
bekini=as.Date(dates[1])
enti=as.Date(dates[lenu])
#print(bekini)
#print(enti)
#stop(-1)
#enti=lenu
#bekini=enti*0+1
#enti=as.integer(enti)
#bekini=as.integer(bekini)
df1 <- setNames(num, names)
mGT<-generation.time("gamma", c(3, 1.5))
#TD <- est.R0.TD(df1, mGT, begin=1, end=length(dates), nsim=200)
#TD <- est.R0.TD(df1, mGT, begin=bekini, end=enti, nsim=200)
TD <- est.R0.TD(df1, mGT, begin=bekini, end=enti, nsim=200)
TD.5D <- smooth.Rt(TD, 5)
paivat1<-TD.5D$epid$t
paivat2<-as.Date(paivat1)
r0t1<-TD.5D$R
conf1<-TD.5D$conf.int
xypaluu<-data.frame(paivat1,r0t1)
names(xypaluu)<-c("paivat","r0")
return(xypaluu)
}
calculate_r0_with_epiestim<-function(xy2)
{
## calculate r0 w/r0 package
dates<-as.Date(xy2$Dates)
cases<-as.integer(xy2$Cases)
nummeros<-1:length(dates)
num<-cases
#names<-nummeros
names<-dates
lenu=length(dates)
cases[is.na(cases)] <- 1
cases[(cases<0)] <- cases*-1
cases[cases==0] <- 1
incid<-cases
bekini=as.Date(dates[1])
enti=as.Date(dates[lenu])
config<-make_config( list(mean_si = 2.6,std_si = 1.5) )
res<-estimate_R(incid,method="parametric_si", config = config)
#plot(res)
resr<-res$R
str(resr)
meanr<-resr$Mean
medianr<-resr$Median
quantile95<-resr$Quantile.0.95
quantile05<-resr$Quantile.0.05
quantile75<-resr$Quantile.0.75
quantile25<-resr$Quantile.0.25
meanr
daydexes<-resr$t_start
daydexes
plot(daydexes, meanr)
dayss<-as.Date(dates[daydexes])
print (dayss)
#stop(-1)
#plot(dayss, meanr)
xypaluu<-data.frame(dayss,meanr)
names(xypaluu)<-c("paivat","r0")
return(xypaluu)
}
calculate_r0_with_simple_exponent_moving_average<-function(xy2, madays1, madays2)
{
## calculate r0 w/r0 package
dates<-as.Date(xy2$Dates)
cases<-as.integer(xy2$Cases)
nummeros<-1:length(dates)
num<-cases
#names<-nummeros
names<-dates
lenu=length(dates)
cases[is.na(cases)] <- 1
cases[(cases<0)] <- cases*-1
cases[cases==0] <- 1
# compute a MA(7)
ma1<-moving_average(cases,madays1,mean)
r0t1<-calculate_multiple_r0(ma1)
r0avg1<-moving_average(r0t1, madays2, mean)
xypaluu<-data.frame(dates,r0t1)
#plot(r0t1)
#print (r0t1)
#stop(-1)
names(xypaluu)<-c("paivat","r0")
return(xypaluu)
}
lataa_thl_tapaukset_kuolleet<-function()
{
url1<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=measure-492118&column=dateweek20200101-508804L"
cube1 <- fromJSONstat(url1, naming = "label", use_factors = F, silent = T)
res01 <- cube11
#res00
url2<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=measure-444833&column=dateweek20200101-508804L"
cube2 <- fromJSONstat(url2, naming = "label", use_factors = F, silent = T)
res02 <- cube21
#res02
#stop (-1)
paiva=as.Date(res01$dateweek20200101)
kuolleet=as.integer(res01$value)
tapaukset=as.integer(res02$value)
kuolin_prosentit=kuolleet/tapaukset
kuolin_prosentit=kuolin_prosentit*10000
kuolin_prosentit=as.integer(kuolin_prosentit)
kuolin_prosentit=as.double(kuolin_prosentit)
kuolin_prosentit=kuolin_prosentit/100.0
#print (paiva)
#print (kuolleet)
#stop(-1)
#print (tapaukset)
#print (kuolin_prosentit )
df1<-data.frame(paiva,tapaukset, kuolleet, kuolin_prosentit)
names(df1)<-c("Paiva", "Tapauksia", "Kuolleita", "Kuolinprosentti")
#write.csv2(df1, "./kuolleet_ikaryhmittain.csv", sep = ";" )
write.csv(df1, "./thl_tapaukset_kuolleet.csv")
xy0<-data.frame(paiva, tapaukset)
names(xy0)<-c("Dates", "Cases")
xy<-na.omit(xy0)
#return(df1)
}
nth_element <- function(vector, starting_position, n) {
vector[seq(starting_position, length(vector), n)]
}
get_thl_hospital_data<-function()
{
url_base1="https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19care/fact_epirapo_covid19care.json"
request1="?row=dateweek20200101-508804L&column=measure-547523.547516.547531.456732.&fo=1"
url1 <- paste0(url_base1, request1)
cube1 <- fromJSONstat(url1, naming = "label", use_factors = F, silent = T)
res1 <- cube11
days0<-as.Date(res1$dateweek20200101)
values1<-as.integer(res1$value)
dimx<-4
dimy<-length(values1)/dimx
values2<-values1
days1<-days0
dim(values2)<-c(dimx,dimy)
dim(days1)<-c(dimx,dimy)
days2<-nth_element(days1, 1, 4)
esh1<-values2[1,]
esh2<-values2[2,] ## esh
esh3<-values2[3,] ## perusterv
teho1<-values2[4,]
esh1[is.na(esh1)] <- 0
esh2[is.na(esh2)] <- 0
esh3[is.na(esh3)] <- 0
teho1[is.na(teho1)] <- 0
#print(esh1)
sairaalassa1<-esh1+esh2+esh3+teho1
#df1<-data.frame(days2, esh2, teho1)
df1<-data.frame(days2, sairaalassa1, teho1)
names(df1)<-c("Paiva", "Sairaalassa", "Teholla")
#write.csv2(df1, "./kuolleet_ikaryhmittain.csv", sep = ";" )
write.csv(df1, "./thl_sairaalassa_teholla.csv")
#xy0<-data.frame(paiva, tapaukset)
#names(xy0)<-c("Dates", "Cases")
#xy<-na.omit(xy0)
return(df1)
}
- main program
if(load_data_from==1)
{
xy<-load_data_from_finnish_wiki()
print (xy)
}
if(load_data_from==2)
{
xy<-load_data_from_aggregated()
}
if(load_data_from==3)
{
xy<-download_solanpaa_finnish_data()
}
if(load_data_from==4)
{
xy<-lataa_thl_tapaukset_kuolleet()
}
if(load_data_from==5)
{
xy<-get_thl_hospital_data()
names(xy)<-c("Dates","Cases")
}
print (xy)
- quit(-1)
#print (beginday1)
select_datelimit_begin=as.Date(beginday1,format="%d/%m/%Y")
select_datelimit_end=as.Date(today1, format="%d/%m/%Y")
#format(select_datelimit_begin, "%Y-%m-%d")
print(select_datelimit_begin)
print(select_datelimit_end)
#2020-12-16
xy2<-xy[xy$Dates >= select_datelimit_begin & xy$Dates <= select_datelimit_end,]
#xy2<-xy[xy$Dates >= select_datelimit_begin,]
print("xy2")
print(xy2)
#stop(-1)
cases1<-xy2$Cases
dates1<-xy2$Dates
len1=length(cases1)
dates2<-as.Date(dates1)
paivat<-1:len1
## test code
arrat0<-calculate_r0_with_simple_exponent_moving_average(xy2, 14,7)
#arrat1<-calculate_r0_with_r0(xy2)
#arrat2<-calculate_r0_with_epiestim(xy2)
#print("calcu ok")
#plot(arrat$paivat, arrat$r0)
# arrat<-arrat2
arrat<-arrat0
str(arrat)
head(arrat)
#plot(arrat$paivat, arrat$r0)
# stop(-1)
-
datelimits2=c(today1, as.Date(forecastendday1,"%Y/%m/%d"))
datelimits3=c(as.Date(beginday1, "%d/%m/%Y" ), as.Date(forecastendday1,"%Y/%m/%d"))
daysek1<-seq(today, as.Date(forecastendday1), "days")
forecastdaynum1<-length(daysek1)
lendaysek1<-length(daysek1)
daysek2<-seq(as.Date(beginday1),as.Date(forecastendday1), "days")
daysek3<-seq(as.Date(beginday1),as.Date(forecastendday1), "days")
daysek4<-seq(as.Date(today), as.Date(forecastendday1,"%Y/%m/%d"),"days" )
y=arrat$r0
print("R in")
y[is.na(y)] <- 1
y[!is.finite(y)] <- 1
basedatay=y
basedatadays=arrat$paivat
print(basedatay)
- stop(-1)
- print (length(y))
- print (forecastdaynum1)
sarima_forecast = sarima.for(y, n.ahead= forecastdaynum1,p=0,d=1,q=1,P=1,D=1,Q=0, S=6)
str(sarima_forecast)
y1=sarima_forecast$pred
y2=sarima_forecast$se
print("YYYYYYYYYY")
- print (y1)
- print (y2)
- plot(y1)
- stop(-1)
print (basedatadays)
print (length(basedatadays))
print (length(basedatay))
- as.Date(a$Number, origin = “1964-10-22”)
#basedata<-c( as.Date(basedatadays, origin="2021-06-01"), basedatay )
#basedata<-c( as.numeric(basedatadays), basedatay )
basedata<-data.frame(basedatadays, basedatay)
names(basedata)<-c("ds", "y")
print(basedata)
# stop(-1)
- farrat1<-forecast_profet(basedata, lendaysek1)
- stop(-1)
# profeta <- prophet(basedata,interval.width = 0.1, yearly.seasonality =0, weekly.seasonality = TRUE)
# profeta <- prophet(basedata,interval.width = 0.1, yearly.seasonality =0, daily.seasonality = TRUE)
# profeta <- prophet(basedata,interval.width = 1.0, yearly.seasonality =1,weekly.seasonality = TRUE, daily.seasonality = TRUE)
# profeta <- prophet(basedata, interval.width = 0.25, daily.seasonality = TRUE)
# profeta <- prophet(basedata, seasonality.mode = 'multiplicative')
profeta <- prophet(seasonality.mode = 'multiplicative')
#profeta <- add_seasonality(profeta, 'quarterly', period = 60, fourier.order = 8, mode = 'additive')
- profeta <- add_seasonality(profeta, 'quarterly', period = 90, fourier.order = 8, mode = 'additive')
- profeta <- add_seasonality(profeta, 'monthly', period = 60, fourier.order = 8, mode = 'additive')
- profeta <- add_seasonality(profeta, 'weekly', period = 30, fourier.order = 8, mode = 'additive')
## profeta <- add_seasonality(profeta, 'p0', period = 365, fourier.order = 8, mode = 'additive')
# profeta <- add_seasonality(profeta, 'p1', period = 180, fourier.order = 8, mode = 'additive')
- profeta <- add_seasonality(profeta, 'p2', period = 90, fourier.order = 8, mode = 'additive')
- profeta <- add_seasonality(profeta, 'p3', period = 30, fourier.order = 8, mode = 'additive')
profeta <- add_seasonality(profeta, 'p1', period = 120, fourier.order = 8, mode ='multiplicative')
profeta <- add_seasonality(profeta, 'p2', period = 60, fourier.order = 8, mode = 'multiplicative')
profeta <- add_seasonality(profeta, 'p3', period = 30, fourier.order = 8, mode = 'multiplicative')
#profeta <- add_regressor(profeta, 'regressor', mode = 'additive')
profeta<-fit.prophet(profeta,basedata)
future <- make_future_dataframe(profeta, periods = lendaysek1 )
forecast1 <- predict(profeta, future)
svg(filename="prophet_r10.svg", width=widthi, height=heighti, pointsize=12)
plot(profeta, forecast1, main="R0 Suomessa", xlab="pvm", ylab="R0",font.main=4, font.lab=4, font.sub=4, cex.main=1.25, cex.lab=3, cex.axis=3)
dev.off()
str(forecast1)
yhat1<-forecast1$yhat
ds<-forecast1$ds
rohve<-data.frame(ds, yhat1)
# stop(-1)
- forecasti
#md = rwf(y,h=forecastlen1,drift=T,level=c(90,95),fan=FALSE,lambda=NULL)
md = rwf(y,h=lendaysek1,drift=T,level=c(30,40),fan=FALSE,lambda=NULL)
str(md)
fitted1<-md$fitted
mean1<-md$mean
lower1<-md$lower[,1]
upper1<-md$upper[,1]
#plot(fitted1)
plot(mean1)
lines(lower1)
lines(upper1)
#stop(-1)
plot(md)
dif_data <- diff(y)
arima1<-auto.arima(dif_data)
foca1<-forecast(arima1)
str(foca1)
plot(foca1)
foca1_fitted1<-foca1$fitted
foca1_mean1<-foca1$mean
foca1_lower1<-foca1$lower[,1]
foca1_upper1<-foca1$upper[,1]
plot(foca1_mean1)
lines(foca1_lower1)
lines(foca1_upper1)
#mdarima1_fitted1<-focal_fitted1
mdarima1_mean1<-mean1
mdarima1_lower1<-lower1
mdarima1_upper1<-upper1
#fit1 <- stl(y, s.window="periodic",t.window=length(y) )
#fit1 <- stl(y )
#arima <- forecast(fit1,h=lendaysek1,method ='arima')
#fit = arima(y, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 7))
#fit = arima(y, c(0, 1,1), seasonal = list(order = c(0, 1,1), period = 7))
fit = arima(y, c(0, 1,1), seasonal = list(order = c(0, 1,1), period = lendaysek1))
pred <- predict(fit, n.ahead = lendaysek1)
mdarima1_pred<-pred$pred
mdarima1_se<-pred$pred
#mdarima1_pred<-pred$se
#mdarima_mean1<-arima$mean
bat1 <- tbats(y)
fc2 <- forecast(bat1, h=lendaysek1)
plot(fc2, ylab="Tapauksia")
sensor <- ts(y,frequency=7) # consider adding a start so you get nicer labelling on your chart.
fit <- auto.arima(sensor)
fcast <- forecast(fit, h=lendaysek1)
plot(fcast)
grid()
print(" Fcast 1 ...")
str(fcast)
mdarima_pred<-fcast$mean
#stop(-1)
#plot(fc)
str(pred)
#plot(pred)
print(mdarima1_pred)
#quit(-1)
# stop(-1)
- r0=1.2
r0=arrat$r0
tau=5.5
log(r0)
k=log(r0)/tau
- print(k)
print(arrat)
print
print("Koo ....")
- sel1=as.Date(today1, format="%d/%m/%Y")
sel2=as.Date(forecastendday1,format="%Y/%m/%d")
sel1=as.Date("2021/08/28", format="%Y/%m/%d")
print (sel1)
print (sel2)
errat<-rohve
- print (errat)
names(errat)<-c("paivat", "r0")
- print (errat)
- str(arrat)
- errat5<-arrat[errat$paivat >= sel1 & errat$paivat <= sel2,]
errat5<-errat[errat$paivat >= sel1,]
- arrat5 <- arrat[arrat$paivat > sel1, ]
- arrat5<-subset(arrat, paivat> sel1 & paivat < sel2 )
- arrat5<-filter(arrat, filter(paivat>sel1 & paivat<sel2))
- arrat5 <- subset(arrat, paivat > "2021-08-28" & paivat < "2021-09-15")
print("ERRAT5")
errat5
ennustepaivat=errat5$paivat
- print (xy2)
print (str(xy2))
xytoday<-xy2[xy2$Dates == sel1,]
print (xytoday)
nykymaara=xytoday$Cases
print (nykymaara)
nykypaiva=xytoday$Dates
- print (length(errat5))
r0=errat5$r0
k=log(r0)/tau
maara=nykymaara
edmaara=maara
print ("FORECSZT")
mamma<-c(nykymaara)
pappa<-c(nykypaiva)
for (i in 1:length(k))
{
kaija=exp(k[i])
maara=edmaara*kaija
paiva=ennustepaivat[i]
print(paiva)
print (maara)
if(i>1) {
print("*")
mamma=append(mamma,maara)
pappa=append(pappa,paiva)
}
edmaara=maara
}
print("Mamma")
print (mamma)
mennuste<-data.frame(pappa, mamma)
names(mennuste)<-c("paivat","maarat")
print (mennuste)
svg(filename="covid_cases_forecast_from_r0.svg", width=widthi, height=heighti, pointsize=12)
ggplot( mennuste, aes(x =paivat , y = maarat) ) +
ggtitle("Koronavirustapauksia - ennuste R0+Prophet") +
xlab("Pvm") + ylab("Uusia tapauksia")+
theme(title=element_text(size=15), axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=14,face="bold"))+
geom_point() +
geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.99, size=3)+
geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) +
geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5)
dev.off()
- r0
- k
stop(-1)
if(tulosta_svg==1)
{
#svg(filename=tulosfilee1, width=8, height=3, pointsize=12)
svg(filename=tulosfilee1, width=widthi, height=heighti, pointsize=12)
}
if(plottaa==1)
{
metodi="loess"
print ("ggplot")
#ggplot(arrat, aes(x =paivat , y = r0)) +ylim(0.6, 1.8)+xlim(as.Date(datelimits1, format="%d/%m/%Y") )+
ggplot(arrat, aes(x =paivat , y = r0)) +ylim(yala, yyla)+xlim(as.Date(datelimits1, format="%d/%m/%Y") )+
ggtitle("Arvioitu koronaviruksen perusuusiutumisluku R0") +
xlab("Kuukausi") + ylab("R0")+
theme(title=element_text(size=15), axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=14,face="bold"))+
geom_point() +
geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.99, size=3)+
geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) +
geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5) +
geom_hline(yintercept=1.0, linetype="dashed", color = "red", size=1)
}
if(tulosta_svg==1)
{
dev.off()
}
Licensing
[edit]- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 13:10, 31 August 2021 | 900 × 450 (90 KB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Width | 720pt |
---|---|
Height | 360pt |