COVID-19 日本国内の潜在的な陽性者数を推定する試み

日本国内の潜在的な陽性者数を推定することは有益ですが、簡単ではありません。PCR検査がランダムになっていないことが推定を難しくしています。有症状者が検査されやすいというselection biasがあるからです。この記事ではいくつか仮定を置いて潜在的な陽性者数を推定したいと思います。

仮定

全国民のうち潜在的に陽性になっている割合

この割合は年代によらず一定と仮定します。ここでは と書きます（posはpositiveの略）。例えば0.0001なら日本人約1億2千万人中、おおよそ12000人が潜在的に陽性になっている計算です。

なお、国民の年代別人口の値はこのページの令和2年3月報 （令和元年10月確定値，令和2年3月概算値） （PDF：301KB） の「2019年10月1日現在（確定値）」の総人口 男女計の値を使用しました。

陽性者中の有症状者の割合

若年層で無症状が多いなど、年代で異なることが知られています。そこで、全員が検査されて診断されているダイアモンドプリンセスのデータをモデルに組み込みます。ここでは と書きます。は年代の添字（1: 0～9歳、 2: 10～19歳、…、 10: 90～99歳）を表します。

有症状者が厚生省の報道発表資料上に観測される割合

この割合は年代によらず一定と仮定します。ここでは と書きます（obsはobservation, symはsymptomaticの略です）。そこそこ高い値であると考えられます。*1

なお、厚生省の報道発表資料はこのページの一番下の方の新型コロナウイルス感染症の国内発生動向（2020年4月3日掲載分）の2ページ目のスライドの右のグラフの値を使用しました。4/2 18:00時点の値になります。

無症状者が厚生省の報道発表資料上に観測される割合

この割合も年代によらず一定と仮定します。ここでは と書きます（obsはobservation, asyはasymptomaticの略です）。有症状者の濃厚接触者など、数が限られており、低い値だと考えられます。

統計モデル

使用した統計モデルは以下になります。

ここで頭文字が大文字の変数はデータが与えられている変数です。各変数の説明は以下です。

• は年代別の国民の人数
• は年代別の陽性者数（推定したいもの）
• は年代別の有症状者数
• は年代別の無症状者数
• は厚生省の報道発表資料の年代別の有症状者数
• は厚生省の報道発表資料の年代別の無症状者数
• は年代別のダイアモンドプリンセスの有症状者数
• は年代別のダイアモンドプリンセスの陽性者数

各事前分布の説明は以下です。

• (1.8)式はが低い値であることへの事前分布
• (1.9)式はが高い値であることへの事前分布
• (1.10)式はが低い値であることへの事前分布
• (1.11)式はが年代に対してある程度滑らかであることへの事前分布

(1.8)式は楽観的、すなわち推定される陽性者数を減らす方向に働いています。色々変えて推定してみると良いと思います。

推定結果と考察

推定総陽性者数の分布は以下になります。分布の中央値でも7000人ぐらいは存在しそうという推定結果になりました。

年代別の推定総陽性者数は以下になります。黒点は中央値、灰帯は95%ベイズ信頼区間です。厚生省の報道発表資料と比べると、無症状かつ陽性者がかなり存在していると予想されます。特に10代以下の潜在的な陽性者はもっと多いと推測されています。

コード

Stanは整数値を推定するのが少し手間であるため、今回はJAGSを用いました。統計モデルを表すファイルは以下です。

JAGSコードでは、二項分布（dbinom）の引数は確率を先に書く点、正規分布（dnorm）の引数は標準偏差ではなく精度であることに注意です。あとは統計モデルのままなので説明を省きます。

推定を実行して図を描くRコードは以下です。ここでは悲観的に「症状確認中」を「有症状」に入れてますが、「無症状」に入れて実行してみるのも良いと思います。「無症状」としてカウントすると、推定総陽性者数の分布の中央値は6800人ぐらいでそこまで変わりませんでした。

その他の推定結果は以下です。

*1:ちなみにこの割合が年代で異なるモデルでも結果に大きな差はありませんでした。

『わけがわかる機械学習』中谷秀洋（著）の書評

僕が中谷さんと初めて会ったのはみどりぼんの読書会で、初めて話したのは岩波DSの打ち合わせだったと思います。今でもそんなに親しくはないと思います。しかし、中谷さんのブログは10年ぐらい前から読んでいました。自然言語処理を中心とする機械学習に関連する理論（の解釈）・論文レビュー・数値実験の記事が多く、他のブログでは見られない独特かつ理解の深い内容で、毎日勉強させてもらっていました。今でも何度も読むべきブログです。その中谷さんが機械学習についてまるごと一冊書いたものが本書になります。もともと買うつもりでしたが、献本いただいたので簡単にご紹介いたします。

わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する

目次は以下になります。

• 0章： はじめに
• 1章： 機械学習ことはじめ
• 2章： 確率
• 3章： 連続確率と正規分布
• 4章： 線形回帰
• 5章： ベイズ確率
• 6章： ベイズ線形回帰
• 7章： 分類問題
• 8章： 最適化
• 9章： モデル選択
• 10章： おわりに
• 付録A： 本書で用いる数学

特長

一読して最初の感想は前半（2章～3章, 5章）が「確率および分布の初中級者向けの完ぺきな入門」、後半（4章～7章）が「わけがわかるPRML」あるいは「声に出して読みたいPRML」です。

残りの部分は、深層学習を含む機械学習全般（ところどころ深層学習にも言及しています）で必須の内容（モデルとは、交差検証、など）を扱っています。

理論の本にありがちな「天下りで定義や式が与えられて、わけわからんけど式変形していく」ではなく、「経験から確率は〇〇という条件を満たしてほしいよね。仮に〇〇として議論を進めるとどうだろう？」という感じに、ボトムアップで話が展開されていくので納得しながら読みやすいです。

5章は4章を拡張するために導入され、8章は7章の問題を解くために導入されます。章のタイトルだけ見ると一見関係性が不明ですが、流れるように議論をすすめるためにこの配置になっています。

数式に関しては0.1節に

本書は「数式がなくてもわかる本」ではありません。

と書いてあるように大学学部レベルの数式が出てきます。しかし、遠慮なく式が出てくるのは4章以降ですし、式変形の内容や考え方が非常に丁寧に書いてあるので、じっくり読めば理解できて力になることは間違いありません。

以降では部分的に補足します。

2～3章

0, 1, 9, 10章は合計して25ページぐらいなのに、2章と3章は合計して65ページもあって力の入れようが分かります。一見普通の章タイトルですが、内容はめちゃくちゃいいです。特に例題をまじえて「同時分布・周辺分布・条件付き確率分布」を説明しているあたりがいいです。

僕のブログの以前の記事で

色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布とヒストグラム、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換、検定、散布図と箱ひげ図、回帰、相関あたりをRやPythonなどを使いながらシンプルに説明していく本があるといいと思うのですが、なかなかバランスのとれたいい本がありません。初歩の初歩しか説明してない、グラフが少ない、検定にページを割きすぎ、分厚い、ちょっと難しいなどの不満点があります。立ち読みして自分にあった本を選ぶのがいいと思います。

と書きましたが、「確率、確率変数、確率密度、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換」に関して僕が想定する内容と比べて200%ぐらいの完成度です。これからは「そのあたりがわかんなかったら、この本の2,3章読んでおいてね」と言うことにします。

4～7章

4.1節 最小二乗法の最後の段落（p.88）で

ここまで、自然な流れに従って「正解」を出したように感じるかもしれません。しかしこの短い話の中で大きな仮定を4つも使っています。仮定1：関数の形が1次式。 仮定2：二乗和誤差が最小＝良い答え。仮定3：誤差独立の仮定。仮定4…（略）

というように前提条件をきっちり説明したうえで、「どこが拡張しやすいか」「拡張したらどうなるか」を続けて論じていきます。一貫してこのような理屈の説明が出てくるのが本書の特徴と思います。

7章は中谷さんの専門の自然言語処理の例題で話がすすむためか、モデルの考え方や苦労などもちらほらあって、勉強になるのはもちろんとても楽しく読める章になっています。

4～7章は、PRMLの3～4章を読もうとしたけど挫折した人、読んだけどイマイチ分かってない人に特にすすめたいです。

Enjoy!

NIMBLEでノンパラベイズを試したメモ

NIMBLEというRのライブラリがあります。BUGS言語風の文法でC++にコンパイルされるタイプの確率的プログラミング言語です。実装されている推定のアルゴリズムはここに書いてあります。MCMCの他にも以下のようなアルゴリズムがデフォルトで実装されており、実行速度もかなり速いです。

• particle filter（Sequential Monte Carlo, SMC）
• 空間モデルのCARモデル（Gaussian Markov Random Fieldsと等価）
• ノンパラベイズのChinese Restaurant Process（CRP）とStick-breaking Process（SBP）

ここでは、ノンパラベイズのCRPとSBPを試してみます。参考にしたのは公式ドキュメントです。

データは以下の以前の記事と同じで、モデルも似たようなものです。

statmodeling.hatenablog.com

はじめにCRP版、次にSBP版を紹介します。

Chinese Restaurant Process（CRP）

library(nimble)

code <- nimbleCode({
  group[1:N] ~ dCRP(alpha, size=N)
  alpha ~ dgamma(1, 1)
  for(c in 1:C) {
    mu_mix[c] ~ dnorm(20, 100)
    s2_mix[c] ~ dinvgamma(5, b)
  }
  b ~ dgamma(0.01, 0.01)
  for(n in 1:N) {
    Y[n] ~ dnorm(mu_mix[group[n]], var=s2_mix[group[n]])
  }
})

set.seed(1)
load('data-nonpara-bayes.RData')
constants <- list(N = data$n, C = data$C)
data <- list(Y = data$velocity)
inits <- list(mu_mix = seq(0, 20, len = constants$C),
              s2_mix = rep(4, constants$C),
              group = sample(1:10, size = constants$N, replace = TRUE),
              b = 1, alpha = 1)

fit <- nimbleMCMC(code = code, constants = constants,
                  data = data, inits = inits,
                  nchains = 4, niter = 10000, nburnin = 2000,
                  summary = TRUE, samplesAsCodaMCMC = TRUE,
                  monitors = c('group', 'mu_mix', 's2_mix' , 'alpha'))

library(coda)
pdf(file='fit-traceplot-CRP.pdf')
traceplot(fit$samples)
dev.off()

• 4行目の group[1:N] ~ dCRP(alpha, size=N) でCRPに従うサンプリングができます。専用のアルゴリズムになっているとのことで、BUGSやStanに比べると動作がだいぶ高速です。
• 6行目の C はクラスター数で、Nよりはだいぶ小さく、ただし、クラスターとクラスターがつぶれて合体しないように少し大きめにとる必要があります（nimbleを実行する際にその旨のwarningが表示されます）。データ（data-nonpara-bayes.RData）に含まれていますが、ここでは20に設定しています。
• 20行目　Stanと異なり、全てのパラメータの初期値をもれなくきちんと与えておかないと、尤度が無限やNAになってRごと落ちたりしますので注意。
• 28行目 samplesAsCodaMCMC = TRUEを設定しておくと{coda}ライブラリで扱いやすい形のサンプリング結果となり、31～34行目のようにしてtraceplotを出力することができます。

その他にも、コンパイルとサンプリングの実行を別にできたりしますが、ドキュメントが若干分かりにくいです。

Stick-breaking Process（SBP）

library(nimble)

code <- nimbleCode({
  for(c in 1:(C-1)) {
    v[c] ~ dbeta(1, alpha)
  }
  alpha ~ dgamma(1, 1)
  w[1:C] <- stick_breaking(v[1:(C-1)])
  for(c in 1:C) {
    mu_mix[c] ~ dnorm(20, 100)
    s2_mix[c] ~ dinvgamma(5, b)
  }
  b ~ dgamma(0.01, 0.01)
  for(n in 1:N) {
    group[n] ~ dcat(w[1:C])
    Y[n] ~ dnorm(mu_mix[group[n]], var=s2_mix[group[n]])
  }
})

(...CRP版と同じ...)

inits <- list(mu_mix = seq(0, 20, len = constants$C),
              s2_mix = rep(4, constants$C),
              group = sample(1:10, size = constants$N, replace = TRUE),
              v = rbeta(constants$C, 1, 1),
              b = 1, alpha = 1)

(...CRP版と同じ...)

CRP版とほとんど同じです。

どちらの場合においても、サンプリング結果からクラスター数や密度を算出する関数を自作する必要があります。またモデルのコードにバグが含まれていると、コンパイルエラーではなくRごと落ちたりするのは若いライブラリのご愛敬です。多少不便でもMCMCの性能はかなり良いと思うので積極的に使いたいと思います。

Enjoy!