データサイエンティスト上がりのDX参謀・起業家

データサイエンティスト上がりのDX参謀・起業家のブログ。データ分析や事業について。自身はアーティスト、経営者、事業家。

コロナ陽性者数の推移への正規関数の当てはめ

コロナの感染者数の推移に関して、感染拡大が上がって下がる変化の仕方を関数で表現できないか、またそれによってある程度拡大は収束の推測ができないかと思っていました。自分が知っている関数の中では、正規分布の確率関数が、拡大と収束の形に一番近いため、正規関数に当てはめて分析してみました。


※この取り組みは、私個人の見解で行っているものであり、また厳密な研究をしているわけではないので、必ずしも正解というわけではなく、あくまでも参考例として捉えてもらえればと思います。



結果的にいうと、正規関数のあてはまりは悪くなく、一定程度の参考になりそうかなというものです。東京の感染者数に対して、いま感染拡大している第6波に当てはめてみると、第6波が始まってから52日程度(2/24あたり)がピークで、7日平均が3.5万人程度という計算になっています(本記事一番下↓の図参考)。7日平均が約3.5万人なので、曜日ピークの水・木曜あたりは4〜5万人くらい?→パラメータを追加して分析したら47日後(2/19あたり)に約2.7万人という結果に


ただこの計算も、この1週間で結構ぶれているので、また数日や1週間程度したら大きく変わるかもしれません。第6波とみられる日から21日間のデータを使って計算したピーク人数は約6.6万人、28日間のデータを使ったピーク人数は約3.3万人と大きく減少しました。この1週間くらい、増加スピードが緩やかになっている影響かと思います。


今回行った計算は、”増えだしてからx日間のデータを使って正規関数を当てはめ、もっとも当てはまりがいいパラメータを選択する”、というものです。この方法を東京の第4波と第5波のデータに当てはめてみたところ、そこそこ当てはまりが良かったです。第4波のときは実際のピークが935人(過去7日間の平均)であるのに対し、推定値が961人〜1,022人。第5波のときは4,923人(同上)であるのに対し、4024人〜7,347人(35〜42日の間で拡大スピードが少し下がった)でした。


この方法が、全国のデータで通用するのか、また世界各国のデータで通用するのかは未確認です。下記、Rで行ったプログラムと、結果の一部です。データは、下記の東京都のオープンデータから取得したcsvを利用しています。


データ取得元:東京都 新型コロナウイルス陽性患者発表詳細 - 東京都_新型コロナウイルス陽性患者発表詳細(全期間一括版) - 東京都オープンデータカタログサイト



第4波を2021/3/18〜6/18としてデータを作成したときの結果。

・3/18から21日(3週間)のデータを使って計算した結果
f:id:isseing333:20220206142356p:plain

縦軸↑:7日平均感染者数
横軸→:開始日からの日数(このグラフは3/18が1、50の時点は5/6)
丸:観測された値
青丸:計算に使ったデータ
黒丸:その後実際に観測されたデータ
曲線:計算した正規関数


・3/18から56日(8週間)のデータを使って計算した結果
f:id:isseing333:20220206142652p:plain




第5波を2021/7/5〜10/6としてデータを作成したときの結果。

・21日間のデータを利用
f:id:isseing333:20220206142905p:plain


・28日間のデータを利用(この時点ではピークが大きく計算されている)
f:id:isseing333:20220206142937p:plain


・35日間のデータを利用(28〜35日間の拡大が少なかったためか、ピークも小さくなった)
f:id:isseing333:20220206143002p:plain


・70日間のデータを利用(若干当てはまりが悪く見える。正規関数でなく、歪正規関数のように、他の関数の方が当てはまりがいいのかもしれない)
f:id:isseing333:20220206143130p:plain




第6波を2022/1/4〜としてデータを作成したときの結果。最大で2/5までの33日間のデータを利用。

・21日間のデータを利用
f:id:isseing333:20220206143409p:plain


・28日間のデータを利用(感染者が急に増えたためか、ピークが約6.6万人と計算)
f:id:isseing333:20220206143439p:plain


・33日間のデータを利用(この1週間くらいは感染拡大が鈍化したためか、ピークが約3.5万人に下がる。2/24日あたりをピークに減少に向かう?拡大の鈍化が続けば、ピークの人数ももっと減るかも)
f:id:isseing333:20220206143458p:plain


↓追記:パラメータ追加して計算するとピークまでの日数とピーク人数が減少(2/19あたりがピーク、7日間平均が約2.7万人、水・木のピークだと3.5〜4万人くらい?)。

f:id:isseing333:20220206200316p:plain



上に書いたように、本記事は正確な研究に基づいたものではないですし、流行学的な感染モデルのようなものを仮定したりしているわけではありません。ただ観測された感染者数だけを使って計算しています。今後の個々人の行動や政策などによって大いに変わる可能性はありますが、ある程度の参考や目安になればと思っています。


Data <- read.csv(".../130001_tokyo_covid19_patients.csv", as.is=T)

table(Data$発症_年月日)

Kansensha <- table(Data$公表_年月日)

#  7日間移動平均を利用→正規関数のブレが少なくなった
Kansensha2 <- filter(Kansensha, rep(1, 7))/7


plot(Kansensha)
plot(Kansensh2)


Prm <- NULL

#  調査するパラメータの設定。当てはまりが良さそうな値の周辺を設定している
#  他の地域や時期を分析するときにはパラメータが大きく変わる可能性はある
peak <- c(100, 200, 300, 500, 1000, 3000, 5000, 10000, 30000, 50000)
mean <- c(45, 49, 50, 51, 52, 53, 55)
sigma <- c(0.05, 0.054, 0.055, 0.056, 0.057, 0.058, 0.059, 0.06)

for(a in 1:length(peak)){
  for(b in 1:length(mean)){
    for(c in 1:length(sigma)){
      Prm <- rbind(Prm, c(peak[a], mean[b], sigma[c]))
    }
  }
}

head(Prm)
dim(Prm)


Ryukou4 <- Kansensha[391:483]
Ryukou5 <- Kansensha[500:593]
Ryukou6 <- Kansensha[683:714]

R <- Ryukou3

Len <- length(R)
Len <- 40



Est <- function(R, Len){
  
  # base <- mean(R[1:7])
  base <- min(R[1:7])
  Ryukou <- R - base
  
  SDS <- NULL
  
  # i <- 1
  
  for(i in 1:nrow(Prm)){
    DS <- NULL
  
    for(x in 1:Len){
      peak <- Prm[i, 1]
      mean <- Prm[i, 2]
      sigma <- Prm[i, 3]
      
      # sigmaが小さい方が分布の幅が大きくなる結果になったので、式のプラスマイナスや分母が間違ってる可能性?
      # ただ、yを計算してプロットしたら当てはまりの良い正規分布になっているので、結果には問題なさそう
      DS <- c(DS, (Ryukou[x] - peak * 1/(2*pi*sigma^2)^0.5 * exp(- (x - mean)^2 / 2*sigma^2))^2)
    }
    
    SDS <- c(SDS, sum(DS))
    # print(i)
  }

  # SDS
  
  print(min(SDS))
  print(Prm[(SDS == min(SDS)), ])
  
  print(max(R))
  
  
  x <- 1:(length(Ryukou)+100)
  
  prm <- Prm[(SDS == min(SDS)), ]
  peak <- as.numeric(prm[1])
  mean <- as.numeric(prm[2])
  sigma <- as.numeric(prm[3])

  # sigma <- 0.1
  # peak <- 200
  
  y <- base + peak * 1/(2*pi*sigma^2)^0.5 * exp(- (x - mean)^2 / 2*sigma^2)
  print(max(y))
  
  
  plot(R, ylim=c(0, max(y*1.1)), xlim=c(0, max(x)))
  points(R[1:Len], col="blue")
  lines(y)
  
  y
}



Ryukou4.2 <- Kansensha2[391:483]
Ryukou5.2 <- Kansensha2[500:593]
Ryukou6.2 <- Kansensha2[683:715]



Est(Ryukou4, 21)
Est(Ryukou4, 35)
Est(Ryukou4, 42)

y4.2.21 <- Est(Ryukou4.2, 21)
y4.2.35 <- Est(Ryukou4.2, 35)
y4.2.42 <- Est(Ryukou4.2, 42)
y4.2.49 <- Est(Ryukou4.2, 49)
y4.2.56 <- Est(Ryukou4.2, 56)
y4.2.63 <- Est(Ryukou4.2, 63)
y4.2.70 <- Est(Ryukou4.2, 70)

y5.2.21 <- Est(Ryukou5.2, 21)
y5.2.35 <- Est(Ryukou5.2, 35)
y5.2.42 <- Est(Ryukou5.2, 42)
y5.2.42 <- Est(Ryukou5.2, 70)

y6.2.14 <- Est(Ryukou6.2, 14)
y6.2.21 <- Est(Ryukou6.2, 21)
y6.2.28 <- Est(Ryukou6.2, 28)
y6.2.33 <- Est(Ryukou6.2, 33)



パラメータ追加

peak <- c(100, 200, 300, 500, 1000, 3000, 4000, 5000, 6000, 7000, 10000, 30000, 50000)
mean <- c(44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 55, 60, 70)
sigma <- c(0.01, 0.05, 0.054, 0.055, 0.056, 0.057, 0.058, 0.059, 0.06, 1)

for(a in 1:length(peak)){
  for(b in 1:length(mean)){
    for(c in 1:length(sigma)){
      Prm <- rbind(Prm, c(peak[a], mean[b], sigma[c]))
    }
  }
}

head(Prm)
dim(Prm)