(ISLR)2019-nCov新增人数预测模型
Title : TODO
Objective : TODO
Created by: chenjunyi
Created on: 2020/2/11
#Feb_11
#截止至昨日24:00 的感染人数
par(family=“STXihei”)
Today_confirmedNum <- c( #截至于Feb_11
44742,42744, 40224, 37162,
34594, 31197, 28060, 24363, 20471,
17238, 14411, 11821, 9720, 7736,
5997, 4535, 2761, 1985, 1297
)
#获得数据长度
date_len <- length(Today_confirmedNum) #数据长度
x <- 1:date_len
lastnum <- date_len
#截至前日24:00 的感染人数
Yesterday_confirmedNum <- c(Today_confirmedNum[2:date_len],830)
#截止至昨日24:00 的新增确诊人数序列
ConfirmNumIcr <- rev(Today_confirmedNum - Yesterday_confirmedNum)
days <- 30
#画出数据离散点 注：plot()必须加上type=""参数才能画出线
plot(x, ConfirmNumIcr,type = ‘b’,xlim = c(1,days+1),ylim =c(0,5000),xlab=‘1月23日起’, ylab=‘新增确诊人数’,cex = 1.5 )
points(x[lastnum], ConfirmNumIcr[lastnum],col=“black”,pch=19,cex=2)
#建立模型；已证明该模型p-value最小
model1 <-lm.fit3_2 <- lm(ConfirmNumIcr~I(x2)+I(x3))
#求出拟合模型
#为什么要用data.frame?
model_date <- data.frame(x=seq(1,days,1))
#predict()函数的意义：给定一个模型model1，给出横坐标向量，生成纵坐标向量，se=TRUE即指给出Standard ERROR
#Q:model_fitted_y也是有下个值的向量，如何lines()出来？
model_fitted_y <- predict(model1, newdata = model_date,se = TRUE)
#画出包括昨日数据的模型
lines(model_date
x
,
m
o
d
e
l
f
i
t
t
e
d
y
x, model_fitted_y
x,modelfittedyfit,col=‘red’,cex = 2,lwd = 2)
#points(model_date
x
[
d
a
y
s
o
f
m
o
n
t
h
]
,
m
o
d
e
l
f
i
t
t
e
d
y
x[days_of_month], model_fitted_y
x[daysofmonth],modelfittedyfit[days_of_month],col=“red”,cex=2,pch=19)
points(model_date
x
[
d
a
t
e
l
e
n
]
,
m
o
d
e
l
f
i
t
t
e
d
y
x[date_len], model_fitted_y
x[datelen],modelfittedyfit[date_len],col=“red”,cex=2,pch=19)
points(model_date
x
[
d
a
t
e
l
e
n
+
1
]
,
m
o
d
e
l
f
i
t
t
e
d
y
x[date_len+1], model_fitted_y
x[datelen+1],modelfittedyfit[date_len+1],col=“red”,cex=2,pch=19)
#画出去掉昨日数据的模型
ConfirmNumIcr.old <- ConfirmNumIcr[1:date_len-1]
x.old <- x[1:date_len-1]
model.old <- lm(ConfirmNumIcr.old~I(x.old2)+I(x.old3)) #得到的是拟合的函数
model_date.old <- data.frame(x.old=seq(1,days,1))
#给出预测变量x，得到在所拟合的函数对应的响应值
model_fitted_y.old <- predict(model.old, newdata = model_date.old,se = TRUE)
lines(model_date.old
x
,
m
o
d
e
l
f
i
t
t
e
d
y
.
o
l
d
x, model_fitted_y.old
x,modelfittedy.oldfit,col=‘green’,cex = 2,lwd = 2)
points(model_date.old
x
[
d
a
t
e
l
e
n
−
1
]
,
m
o
d
e
l
f
i
t
t
e
d
y
.
o
l
d
x[date_len-1], model_fitted_y.old
x[datelen−1],modelfittedy.oldfit[date_len-1],col=“green”,cex=2,pch=19)
points(model_date.old
x
[
d
a
t
e
l
e
n
]
,
m
o
d
e
l
f
i
t
t
e
d
y
.
o
l
d
x[date_len], model_fitted_y.old
x[datelen],modelfittedy.oldfit[date_len],col=“green”,cex=2,pch=19)
#备注信息
legend(“top”, legend = “空心点为实际数据，绿线为昨日模型，红线为今日模型”,text.col = “black”,text.width = 26)
text(date_len+7,model_fitted_y
f
i
t
[
d
a
t
e
l
e
n
+
1
]
,
"
预
测
明
日
(
2.13
)
新
增
点
"
)
t
e
x
t
(
d
a
t
e
l
e
n
+
7
,
m
o
d
e
l
f
i
t
t
e
d
y
.
o
l
d
fit[date_len+1],"预测明日(2.13)新增点") text(date_len+7,model_fitted_y.old
fit[datelen+1],"预测明日(2.13)新增点")text(datelen+7,modelfittedy.oldfit[date_len],“昨日模型预测今日新增点”)
diff <- ConfirmNumIcr[date_len] - model_fitted_y.old
f
i
t
[
d
a
t
e
l
e
n
]
d
i
f
f
p
e
r
<
−
d
i
f
f
/
C
o
n
f
i
r
m
N
u
m
I
c
r
[
d
a
t
e
l
e
n
]
p
r
i
n
t
(
s
p
r
i
n
t
f
(
′
昨
日
模
型
偏
差
:
p
r
i
n
t
(
s
p
r
i
n
t
f
(
′
今
日
模
型
预
测
明
日
新
增
数
据
为
:
fit[date_len] diff_per <- diff/ConfirmNumIcr[date_len] print(sprintf('昨日模型偏差: %4.2f%%', 100*diff_per)) print(sprintf('今日模型预测明日新增数据为: %4.2f', model_fitted_y
fit[datelen]diffper<−diff/ConfirmNumIcr[datelen]print(sprintf(′昨日模型偏差:print(sprintf(′今日模型预测明日新增数据为:fit[date_len+1]))
print(sprintf(‘算上偏差，今日模型预测明日(2.13)新增数据为: %4.2f’,(diff_per+1)*model_fitted_y$fit[date_len+1]))
points(model_date$x[date_len+1], 1275.72,col=“blue”,cex=2,pch=19)