接触確認アプリ(COCOA)の実機デバッグ環境を整えようとするができなかった話

最初に書いておきますが、iPhoneでの接触確認アプリの実機環境は整えることができません。もともと暴露通知(Exposure Notification)が、ひとつの国でひとつのアプリでしか使えないという制限のためという理由もあるのですが、Exposure Notification を有効にしたまま実機(iPhone)にインストールすることができません。

目的としては、

  • 暴露通知(Exposure Notification)の動作を把握する
  • COCOA の不具合っぽいものを確認しておく(直すことはできないので、確認のみ)

とするので、内部的なところには踏み込まないようにします。

ソースコードをダウンロードする

cocoa の元ネタのソースコードを github からダウンロードします。

現在の cocoa は ver.1.1.2 なのですが、ここにあるソースコードは 1.1.1 のままになっています。このため、実際の cocoa とは違うので「動作チェック」という訳にはいかないのですが、Exposure Notification の仕組みを知るには良い材料でしょう。

プログラムが Xamarin.Forms を使っているので、Visual Studio 2019 と Visual Studio for Mac を使います。

Visual Studio 2019 のほうは、Windows 上で Android アプリのビルドができます。iPhone アプリの開発は Windows から mac を通しても可能なのですが、mac のほうに Visual Studio for Mac を入れたほうが便利です。

ソースコードを開く

Android と iOS のアプリは、Covid19Radar.sln を Visual Studio で開きます。

  • Covid19Radar.sln
  • Covid19Radar.Functions.sln

Covid19Radar.Functions.sln のほうは、サーバー側のプログラムです。きちんと設定してやれば、ローカルの Azure エミュレータと Cosmos DB を使ってサーバーサイドも構築することも可能なのですが、これはまだ試していません。

  • Covid19Radar
  • Covid19Radar.Android
  • Covid19Radar.iOS

という3つのプロジェクトがあります。Xamarin.Forms で Android/iOS のアプリを作るとき、共通コードとなる Covid19Radar プロジェクトがあって、それぞれ機種特有のコードを Covid19Radar.Android や Covid19Radar.iOS に書くことになっています。

肝心の Exposure Notification のコードは、内部的に Google(Android) と Apple(iOS) では別々のコードになるところなのですが、XamarinComponents/XPlat/ExposureNotification at master · xamarin/XamarinComponents という形で、Android と iOS の動作をひとつにまとめて扱えるようになっています。

NuGet の場合は、以下からダウンロードができます。

setting.json を書き替える

Exposure Notification の仕組みで優れているところは、

  • 自分が陽性と解ったときだけ、サーバーに TEK(Temporary Exposure Key)を送る
  • 大勢の人は、定期的に陽性者の TEK をダウンロードして、自分のスマホ内で照合する

ところです。個人データと思われる TEK の情報を、「陽性情報の登録」のときのみサーバーにアップロードします。HER-SYS 番号との連係はさておき、自発的に「陽性情報の登録」をしたときだけしか、サーバーにアクセスしにいかないので、個人データの漏洩リスクが減ります。

もうひとつ、定期的にサーバーから TEK をダウンロードしますが、この中にある RPI(Rolling Proximity Identifier)だけでは、個人を特定しにくいです。実は、全くできないわけではなくて、別途 Beacon を使って接触通知アプリを使って広域的に探っていけばできないこともないのですが、現在のところ難しい、というところです。

Covid19Radar プロジェクトの中に、setting.json というファイルがあります。

これを開くといくつかの設定があります。実際の cocoa では、APP_VERSION などが設定されているはずです。

{
  "appVersion": "APP_VERSION",
  "apiSecret": "API_SECRET",
  "apiUrlBase": "https://API_URL_BASE/api",
  "supportedRegions": "440",
  "blobStorageContainerName": "c19r",
  "androidSafetyNetApiKey": "ANDROID_SAFETYNETKEY",
  "cdnUrlBase": "https://CDN_URL_BASE/",
  "licenseUrl": "https://covid19radarjpnprod.z11.web.core.windows.net/license.html",
  "appStoreUrl": "https://itunes.apple.com/jp/app/id1516764458?mt=8",
  "googlePlayUrl": "https://play.google.com/store/apps/details?id=jp.go.mhlw.covid19radar",
  "supportEmail": "SUPPORT_EMAIL"
}

apiUrlBase は、「陽性情報の登録」をしたときの呼び出し先で、Azure Functions が使われています。このアドレスと秘密キーは本番の cocoa だけのもなので、このままにしておきます。逆に言えば、、間違って「陽性情報の登録」を押しても、勝手にサーバーに接続できないということですね。

陽性者の TEK をダウンロードするための URL は、cdnUrlBase となるので、ここは書き替えます。

  "cdnUrlBase": "https://covid19radar-jpn-prod.azureedge.net/",

デバッグページ(DebugPage)を表示させる

もうひとつ、動作確認用のデバッグページを表示できるようにします。

ViewModels/MenuPageViewModel.cs を開いて、以下のコードを追加しておきます。

#if DEBUG
            MenuItems.Add(new MainMenuModel()
            {
                Icon = "\uf0c0",
                PageName = nameof(DebugPage),
                Title = nameof(DebugPage)
            });
#endif

Debug_Mock と Debug の違い

プロジェクトをビルドする前に、Debug_Mock と Debug を確認しておきます。

  • Debug は、Exposure Notification API を使ったデバッグモード
  • Debug_Mock は、Exposure Notification API を使わずにエミュレート

の違いがあります。

Nearby.EXPOSURE_NOTIFICATION_API エラーが出る

EN api を使ったまま、実機 Android にインストール(Visual Studio からのデバッグ実行)はできるのですが、何かのタイミングで、次のようなエラーが出ます。

**Android.Gms.Common.Apis.ApiException:** '17: API: Nearby.EXPOSURE_NOTIFICATION_API is not available on this device. Connection failed with: ConnectionResult{statusCode=UNKNOWN_ERROR_CODE(39507), resolution=null, message=null}'

前後のログ出力はこんな感じ。

08-22 09:26:56.703 I/MonoDroid(23458): UNHANDLED EXCEPTION:
08-22 09:26:56.713 I/MonoDroid(23458): Android.Gms.Common.Apis.ApiException: 17: API: Nearby.EXPOSURE_NOTIFICATION_API is not available on this device. Connection failed with: ConnectionResult{statusCode=UNKNOWN_ERROR_CODE(39507), resolution=null, message=null}
08-22 09:26:56.713 I/MonoDroid(23458):   at Android.Gms.Nearby.ExposureNotification.IExposureNotificationClient.IsEnabledAsync () [0x00067] in <a915b6bc332b428d88e3e93dc94335a6>:0 
08-22 09:26:56.713 I/MonoDroid(23458):   at Xamarin.ExposureNotifications.ExposureNotification.PlatformGetStatusAsync () [0x00074] in <e42d902759884cb2b9bcb6b7e1a5859c>:0 
08-22 09:26:56.713 I/MonoDroid(23458):   at Covid19Radar.Services.ExposureNotificationService.UpdateStatusMessageAsync () [0x00025] in D:\git\cocoa\Covid19Radar\Covid19Radar\Covid19Radar\Services\ExposureNotificationService.cs:88 
08-22 09:26:56.713 I/MonoDroid(23458):   at Covid19Radar.Services.ExposureNotificationService.OnUserDataChanged (System.Object sender, Covid19Radar.Model.UserDataModel userData) [0x00057] in D:\git\cocoa\Covid19Radar\Covid19Radar\Covid19Radar\Services\ExposureNotificationService.cs:73 
08-22 09:26:56.713 I/MonoDroid(23458):   at System.Runtime.CompilerServices.AsyncMethodBuilderCore+<>c.<ThrowAsync>b__7_0 (System.Object state) [0x00000] in /Users/builder/jenkins/workspace/archive-mono/2020-02/android/release/mcs/class/referencesource/mscorlib/system/runtime/compilerservices/AsyncMethodBuilder.cs:1021 
08-22 09:26:56.714 I/MonoDroid(23458):   at Android.App.SyncContext+<>c__DisplayClass2_0.<Post>b__0 () [0x00000] in <eaa205f580954a64824b74a79fa87c62>:0 
08-22 09:26:56.714 I/MonoDroid(23458):   at Java.Lang.Thread+RunnableImplementor.Run () [0x00008] in <eaa205f580954a64824b74a79fa87c62>:0 
08-22 09:26:56.714 I/MonoDroid(23458):   at Java.Lang.IRunnableInvoker.n_Run (System.IntPtr jnienv, System.IntPtr native__this) [0x00008] in <eaa205f580954a64824b74a79fa87c62>:0 
08-22 09:26:56.714 I/MonoDroid(23458):   at (wrapper dynamic-method) Android.Runtime.DynamicMethodNameCounter.1(intptr,intptr)
08-22 09:26:56.714 I/MonoDroid(23458):   --- End of managed Android.Gms.Common.Apis.ApiException stack trace ---
08-22 09:26:56.714 I/MonoDroid(23458): com.google.android.gms.common.api.ApiException: 17: API: Nearby.EXPOSURE_NOTIFICATION_API is not available on this device. Connection failed with: ConnectionResult{statusCode=UNKNOWN_ERROR_CODE(39507), resolution=null, message=null}
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.internal.ApiExceptionUtil.fromStatus(com.google.android.gms:play-services-base@@17.1.0:4)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.api.internal.ApiExceptionMapper.getException(com.google.android.gms:play-services-base@@17.1.0:2)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.api.internal.zaf.zaa(com.google.android.gms:play-services-base@@17.1.0:15)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.api.internal.GoogleApiManager$zaa.zac(com.google.android.gms:play-services-base@@17.1.0:175)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.api.internal.GoogleApiManager$zaa.onConnectionFailed(com.google.android.gms:play-services-base@@17.1.0:95)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.internal.zag.onConnectionFailed(com.google.android.gms:play-services-base@@17.1.0:2)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.internal.BaseGmsClient$zzf.zza(com.google.android.gms:play-services-basement@@17.2.1:6)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.internal.BaseGmsClient$zza.zza(com.google.android.gms:play-services-basement@@17.2.1:25)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.internal.BaseGmsClient$zzc.zzo(com.google.android.gms:play-services-basement@@17.2.1:11)
08-22 09:26:56.714 I/MonoDroid(23458): 	at com.google.android.gms.common.internal.BaseGmsClient$zzb.handleMessage(com.google.android.gms:play-services-basement@@17.2.1:49)
08-22 09:26:56.714 I/MonoDroid(23458): 	at android.os.Handler.dispatchMessage(Handler.java:107)

要は、EN api を利用するときには、アプリが正式に Google Play に登録しておく必要があり、アプリが EN api を使うときに Google Play Service が呼び出しのチェックをしているそうです。

このあたりのデバッグのしにくさは、ドイツ版の接触確認アプリ(corona-warn-app)や Google の exposure-notifications-android でも話題になっています。

ちなみに、corona-warn-app のほうでは、Frida を使った Google Play Service の偽装の方法が提案されているので、これだと上手くいくかもしれません(ちょっと試してない)。

Google Play 開発者サービスでガードが掛かる

つまり、Android 実機の場合は、

  • EN api を有効したままインストールはできるが、実行時に EN api を呼び出したときに Google Play Service でガードがかかる。

という訳で、UI 周りのチェックをしたいときは Debug_Mock で動かし、EN api の挙動をを調べる場合は今のところ手がないというパターンです。

iOS の場合は、インストール時にガードが掛かる

iPhone 版の場合は、Entitlements.plist に EN api の記述があります。

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
<plist version="1.0">
<dict>
	<key>com.apple.developer.exposure-notification</key>
	<true/>
	<key>aps-environment</key>
	<string>development</string>
</dict>
</plist>

この「com.apple.developer.exposure-notification」の値が true になることで、iPhone で EN api を利用できるようになります。

ただし、ビルドして実機に iPhone にインストールするときにプロビジョニングでこの値をチェックするのですが、はねられてしまいます。この値を有効にするには、あらかじめ申請が必要で
Exposure Notification Entitlement Request – Contact Us – Apple Developer で送ります。これはい各国の保健省(日本では厚生労働省)の認可が必要なので、日本の場合は cocoa の開発サイドのみになります。

カテゴリー: 開発 | 接触確認アプリ(COCOA)の実機デバッグ環境を整えようとするができなかった話 はコメントを受け付けていません

PCR検査の偽陽性問題を具体値で考える

PCR検査の拡充に対して偽陽性の問題がたびたび出てくる。8月現在で陽性率の高さもあって、再びPCR検査を増やすあるいは増やさないの意見が飛び交うわけだが、具体的に既にPCR検査を行っている国の場合、どの程度の「偽陽性」を許容しているのかを計算してみることにしよう。

感度と特異度

まずは用語を明確にしておく。

検査/東京大学 保健センター

ベイズだと感度も特異度も同等に扱うのだけど、疫学では別々に扱う。

  • 感度が、検査によって真の陽性者が陽性と判明する確率
  • 特異度が、検査によって真の陰性者が陰性と判明する確率

を表す。東大保健センターの図をそのまま借用すると以下のようになる。

ここで問題になるのは、

  • 検査によって真の陽性者ではあるが陰性として判別される「偽陰性」
  • 検査によって真の陰性者ではあるが陽性として判断される「偽陽性」

の2種類がある。この数値が増加することによる問題はそれぞれ別になる。

  • 偽陰性が増加すると、市中に陰性と思っている陽性患者を放つことになる → 二次感染の増加
  • 偽陽性が増加すると、陽性と思っている陰性者(正常者)が病院に収容することになり病床を圧迫する → 医療崩壊の可能性

「陽性的中率」は 70/(70+9) = 検出した真の陽性者数 / 検査での陽性者数 となる。陽性的中率が悪ければ、罹患していない人を患者として病院に収容することになる。これも病床が埋まり、医療崩壊の原因となる。

新型コロナウィルスのPCR検査では、どの程度の感度と特異度なのかが問題なのであるが、医療機関の解説を見ると、以下のようになる。

  • 感度 70%
  • 特異度 99%

特異度が極めて高いのは、偽陽性を多くして余計な治療を行わないように注意するためだ。通常の医療であっても、薬の投与には副作用があり健康体に害が多い。ゆえに、無闇に薬を出さないようにするため(正常体を間違って副作用で病気にさせないため)に、特異度を極めて高く=偽陽性を極めて低くなるように慎重に調査する。

感度が特異度よりも低いのは、病気の疑いがある場合も患者を拾うためだ(と思う)。PCR検査の場合は、ウィルスの数を敷居とするため、無闇に高い閾値にする≒感度を上げると本来のウィルスで病気になる人を見落としてしまう恐れがある。それよりも、まずは閾値を少し下げて≒感度を下げて病気の疑いがあることを示す。それから再検査という手順になる。よくある人間ドックの精密検査がその方式で、人間ドックでは大まかな病気疑いを血液検査やバリウム検査などで異常値を取得する。このときに、ちょっと疑いありの状態では精密検査として後日再調査をする。
ときには、精密検査をしても何も出ずに健康体であったことがわかる。しかし、何分の一の確率で実際に病気であることがわかる。このとき、最初の人間ドックの検査が無駄であったという訳ではなく、精密検査をするためのよい閾値となっていることがわかる。
この方式を「スクリーニング」と呼ぶ。

スクリーニングの手法は、医療関係だけでなくテスト手法や製造工程の以上検出などでよく使われる。製造工程が異常であること(不良品を出さないという意味で)を素早く検知したいが、無闇に工程を止めるのも時間と資金の無駄であるのでできるだけ工程を止めたくないというジレンマがある。検出器のスピードや精度にも関わってくる。この場合、不良品検出のための閾値を低めに設定しておき、不良品疑いのものをもう一度検品するという方法をとる。最初の不良品検出はざっと行うスクリーニングになる。苦情受付の電話窓口と裏に控えている部長の関係もスクリーニングになる。

このように不良であること、この場合は陽性として検出して患者として治療にあたりたいという目的に沿う場合、二段階の検査が有効に働く。

日本のPCR検査の場合、CTスキャンによる検査や医者の問診が良いスクリーニングになっている。CTスキャンによる肺炎疑いや、専門医よる病状の早く、あるいは濃厚接触者というスクリーニングを経て、PCR検査という手順を踏む。当然、PCR検査自体もできるだけ陽性患者を拾う(隔離という意味でも)たい意図があるので、感度はあまり上げていないのだと思われる。
この感度が90%だったり、特異度を99%近くに上げることが妥当かどうかはよくわからない。ただ、ウィルスの数を閾値としているのだから、ある程度の幅を持たせないと意味がない。少しのウィルスでも病気になるのか、それとも大量のウィルスでないと病気にならないのかという違いがある。単純に言えば、1個だけの新型コロナウィルスが肺に入っても病気にはならないだろう。1個だけなのに急激に増殖するとは考えにくい。ならば10個ならばどうか、100個ならばどうか、という形で閾値を決めることになる。そこには確率的に罹患するかどうかの境目があるので、明確な数があるわけではない。

これにより、感度と特異度は、凡その値として以下を使ってもよいのだろうと思う。

  • 感度 70%
  • 特異度 99%

このスクリーニングが結構重要な役割をしていることが解る。ベイズで言えばこれが「事前条件」となる。スクリーニングがあるなしでは、この事前条件が異なるため、その後の確率の見方が異なる。

一方で陽性率の計算として

  • 陽性率 = 陽性者数 / PCR検査数

がよく使われる。陽性率が高いと危険でPCR検査の数が足りていない、というのがPCR検査を多くする場合の理由の一つとなるのだが、

  • スクリーニングの精度が良く、ほどよく感度を上げている場合には、陽性率は高くなる
  • スクリーニングの精度が悪く、患者疑いの人を大量に組み上げている場合には、陽性率は低くなる

という特徴がある。例としては、前者がCTスキャン等でスクリーニングを行った場合、後者が自己診断によるPCR検査の希望、という形になる。
ただし、「スクリーニングの精度」にも問題があり、例えば「確実に、新型コロナウィルスに罹った患者しか通さないスクリーニング」を施した場合でも、陽性率は高くなる。4月頃に言われていた、37度が4日間継続というのはそれに近いところがある。
故に、陽性率の高い低いだけでは、PCR検査の数が妥当であるかどうかは解らないということが、解るのである。

ニューヨーク市の偽陽性を推測する

さて、ほどんどCTスキャン検査をしていない(と思う)状態で、PCR検査を最大化している(と思われる)ニューヨーク市の状態で偽陽性の数を計算してみよう。


ニューヨーク州の人口は1945万人で、抗体検査によりニューヨーク市は12%程度は既に感染済みであることがわかっている。本来は市のほうで計算したかったのだが、検査済み検査数が州のほうにしかないので「ニューヨーク州」で計算する。

8/1 時点

  • 感染者数 754 人
  • 実施済み検査数 82,737 人

ほんとうは実際に検査した中から感染者数を割り出さないといけないのだが、今回は概算なのでこの数値を使う。

同様に東京都の8/1では

  • 感染者数 472 人
  • 検査数 4324 人

となっている。

これを先の表に埋め込む

  • 検査で陽性の計に「感染者数」を入れる
  • 検査で陰性の計に「検査数 – 感染者数」を入れる。
  • 合計値が「検査数」になる。

このとき、真の罹患している/罹患していない人数を計算したい。
ここで、初期値として、検査が100%正しいと仮定して、

  • 罹患している人 = 感染者数
  • 罹患していない人 = 検査数 – 感染者数

を入れる。当然のことながら、合計値は同じになる。

ここで、検査の仮定である「感度」と「特異度」を入れる。

  • 感度 70%
  • 特異度 99%

罹患している人の数から、感度の割合で、検査の陽性/陰性を割り振る。
罹患していない人の数から、特異度の割合で、検査の陽性/陰性を割り振る。

検査の陽性の合計を、「推測」の列に入れる。
同じように検査の陰性を「推測」の列に入れる。

ここで、罹患している/していない人数が、正しいければ、計の列と推測の列は一致するはずなのだが、当然のことながら一致しない。これは、感度と特異度があるから、当然のことなのだ。検査は実際の真の値とはズレる。

これを違いとして列をつくっておく。

さて、初期値である罹患している人数が違っていることが分かったので、これを真の値に近づけてみよう。
本来ならば、最小二乗法などを使いながら近づけるのだけど、ぽちぽちと「違い」を見ながら値を近づけていくだけでも可能だ。

東京都の場合は、最初の472から620まで増やしていくと「違い」の列がほぼゼロになることがわかる。つまり、この検査では、

  • 偽陰性は186人ほど取りこぼし、
  • 偽陽性で37人ほど間違って隔離している

ことがわかる。

同じことをニューヨーク州で実験してみよう。
東京都の場合と違い、罹患している患者(オレンジのセル)を少しずつ下げていくと違いが減少していくことがわかる。この違いが「0」になれば、真の罹患している患者数がわかるはずだ。
しかし、ニューヨーク州の場合は、罹患している人数を 10 人にしても、違いが0にならない。
これはどういうことだろうか?

  • PCR検査の感度や特異度の仮定が間違っているため、計算があわない。
  • PCR検査が多すぎるため、偽陽性多すぎて、正しく検査できていない。

のどちらかだろう。

このシミュレーションでは、ニューヨーク州でのPCR検査では、罹患している人を正しく検出できていないことが判明する(感度や特異度の仮定が正しければ)。つまり、感度が70%、特異度が99%程度では、真の罹患している患者が0人であっても、検査で陽性となる患者が700人以上で出てしまうということだ。

特異度を 99.9% と仮定する

実際の特異度はもっと高くて、ほぼ100%に近いという記事もある。

RIETI – PCR検査体制の拡充と偽陽性の問題

では、特異度を10倍精度良くして、99.9%と99.99%で考えてみよう。

この位の特異度になると、PCR検査は有効に働く。

  • 特異度が 99.9% の場合は、真の罹患者が 970 人程度
  • 特異度が 99.99% の場合は、真の罹患者が 1070 人程度

と推測される。つまり、偽陽性が極端に少ない状態=特異度の精度が高い状態の場合には、大量なPCR検査は有効に働くということだ。

所感

実際の感度がどの位なのかは解らないが、特異度が極めて高い(偽陽性が限りなく低い)状態でないと、大量PCR検査は有効に働かない。
また、感度の関係から、検査をしても一定量の真の感染患者を取りこぼしてしまう。取りこぼしをどうするのかという問題が残ってしまう。

対策としては、「検査で陰性であっても、偽陰性を疑い一定期間隔離する」ことになるのだが、この表に出てくる「検査で陰性」かつ「罹患している」患者というのは具体的に知ることはできない。「検査で陰性」の人数に含まれてしまっている。このため、偽陰性を疑い、すべての人を隔離してしまうとそれだけで病床が埋まってしまい医療崩壊を起してしまう。難しいところだ。

このため、偽陰性自体を減らす方法として、

  • 検査の総数自体を減らして、偽陽性となる人数を減らす

ことになるので、「検査の総数」を減らすために事前のスクリーニング(事前条件となるCTスキャンや医師による判断、濃厚接触者の判断)が必要になると思われる。このあたりのシミュレーションは別途行いたい。

参考文献

カテゴリー: 開発 | PCR検査の偽陽性問題を具体値で考える はコメントを受け付けていません

隠れ SIR モデルによる感染予測の解説

日頃ツイッターで更新しているこの図の説明がまとまっていなかったので、記録しておきます。

過去の値から感染者を推測する

発想元は西浦教授のこの論文 感染症流行の予測:感染症数理モデルにおける定量的課題 からです。現在分かっている新規感染者の数から、過去に遡って感染者数と感染率を推測し、再び現在と未来の感染者数を予測するという方法です。正確な計算では、過去の感染率(あるいはRt値)は分布を持つので、未来の推測も幅を持った分布になるわけですが、そこまで計算するのは大変なので、簡略的にExcelで可能なように書き直します。

計算方法

  1. 新規感染者dI_nを7日間(潜伏+発症期間)遡って、隠れ新規dI’_n-7と定義する
  2. 隠れ患者I’nを続く7日間の合計値と定義する。
  3. 隠れ新規dI’n/隠れ患者I’nを隠れ感染率rと定義する。
  4. 隠れ実効再生産数Rt’を隠れ患者と続く7日間の隠れ新規の合計値の比と定義する。

なお、潜伏+発症期間を7日間とする理由は、現在の新型コロナウィルスでは感染期間が6.7日間と発表されていることが理由です。

特別な行動変容(都市封鎖や自粛要請など)がない限り、感染率や実効再生産数は変わらないと考えられる。この仮定から、過去7日間の隠れ感染率rをその前日の感染率とほぼ同じと推測できる。

これにより、過去7日間の隠れ新規dI’nが推測できる。
隠れ新規dI’nは7日後に新規患者数dInとして現れる関係があるため、過去7日間の隠れ新規n’から未来の日にあたる新規患者数nが推測できる。

未来においては、隠れ患者I’nの値を、新規患者dInと隠れ新規dI’nの値から計算でき推測が可能である。

数式にすると何やらややこしいのですが、要するに、

  • 潜伏期間+発症期間を7日間とする
  • この期間で、二次感染が発生する
  • 現在の新規患者は7日前に何処かで感染している。
  • どこかで感染した後に、自覚症状が出て新規感染として隔離される前に、どこかで二次感染を発生させる。

というように、感染した直後から新規感染者としてカウントされる前に「滞留」する現象に注目します。実際、SIR や SEIR モデルでは回復の期間に二次感染を発生させるというモデルです。

Excel で計算する

これを Excel で計算できるようにします。感染率の計算が循環関数にならないように1日ずらしているのがミソです。循環関数になるとイテレーションが発生して、計算がややこしくなってしまうため、こうして1日ずらしています。

実際にはそのままプロットすると、土日の検査がないときに凹んでしまうので、7日間の移動平均でプロットします。

7日前よりも後にある感染率rは、推測値になるため手動で設定しています。先に解説した通り、特別な行動変容の変化がない限り、実効再生産数や感染率は大きく変わらないので、直前の値の平均などを使います。

実験として、都市閉鎖を行った場合を想定し、感染率rを 0.100 程度に変化させることで予想グラフが作成できます。

現在の感染率のままの場合

都市封鎖をして 7/13 より感染率を 0.100 に変えた場合

オレンジ色の新規患者のラインが変わることがわかります。

東京都以外でも計算してみる

日次データがあれば、他の都市の様子を観察できます。

ロサンゼルスの場合

テキサス州のダラスの場合

ドイツのベルリンの場合

灰色のラインが隠れ患者数です。テキサス州のダラスでは単調増加となり全く抑え込めていません。ロサンゼルスでは、7/1頃に少し抑え込むことに成功しましたが、再び増加中です。

ドイツのベルリンでは、直近では増加していますが、隠れ患者の数が上下に揺れています。これが抑え込めている状態です。直近のデータが6/30までしかないので、最近のものは不明です。

この推測計算の欠点

この概算では計算できないものを上げておきます。

  • 院内感染数を考慮しない。
  • 外部からの流入を考慮しない。
  • 未来予測では、新規感染者を一定率で取り除けることと仮定している。
  • 自然回復が考慮されていない。

もともと分布のある推測値をそのままプロットしてしまうので、推測される新規患者数には範囲があります。ただし、将来において同じ感染率とは限らない(自粛や PCR検査強化、何らかの検査の偏りなど)ので、もともと実際の新規患者数とは一致しません。なので、ひとまずの増加傾向と対処を行ったときの遅延(滞留する患者による二次感染のため)を確認するためという点では使えると思います。

参考データ

計算式が入った Excel は以下でダウンロードができます。

各国のデータは以下から取得しています。

カテゴリー: 開発 | 隠れ SIR モデルによる感染予測の解説 はコメントを受け付けていません

Vue.js で子コンポーネントから親に値を送る(SVG編)

この手の記事は結構あるのだけど、毎度忘れてしまうので備忘録的に。

やりたいこと

こんな風にSVGにピンを置きます。Vue.js を使って、SVG画像のピンの位置を親の App.vue に伝えます。

App.vue

親コンポーネント(App.vue)のデータで持つのは map_x と map_y で、data() の戻り値で保持します。これは、そのまま App.vue 内で {{map_x}} で参照が可能です。

方眼紙とピンで表示は Hello.vue で指定できるようにします。
Hello.vue 側は max_x と max_y プロパティ、変更したときのイベントを updatePos で返せるようにしてあります。

<template>
  <div id="app">
    <div>map_x: {{map_x}}</div>
    <div>map_y: {{map_y}}</div>
    <Hello 
     :map_x="map_x" 
     :map_y="map_y" 
     @updatePos="updatePos" />
  </div>
</template>

<script>
import Hello from './components/Hello.vue'

export default {
  name: 'App',
  components: {
    Hello,
  },
  data() {
    return {
      map_x: 100,
      map_y: 100,
    };
  },
  methods: {
    updatePos(x,y) {
      this.map_x = x ;
      this.map_y = y ;
    }
  }
}
</script>

Hello.vue

親の App.vue から max_x と map_y を受け取るのがプロパティで props に記述します。ピンはマウスのクリック(v-on:mousedown)で動かせるようにしてあるので、別途データとして pin_x と pin_y を用意します。
プロパティの max_x と max_y は親から引き渡される読み取り専用のパラメータなので別途代入しています。

実際は、pin_x と pin_y は、ピンの先端に座標が合うように少し補正してあります。

<template>
  <div>
    <h1>Hello SVG</h1>
    <div>
      <svg
           ref="pic"
           viewBox="0 0 600 400"
           v-on:mousedown="mousedown()"
      >
            <image 
              href="@/assets/images/hogan.jpg" />
            <image
               :x="pin_x"
               :y="pin_y"
               width="20px"
               height="58px"
               href="@/assets/images/marker.svg"
            />
      </svg>
    </div>
  </div>
</template>

<script>
export default {
  name: 'Hello',
  props: {
    map_x: Number,
    map_y: Number,
  },
  data() {
    return {
      pin_x: this.map_x - 10 ,
      pin_y: this.map_y - 40 ,
    }
  },
  methods: {
    mousedown() {
        const vbox = this.$refs.pic.viewBox ;
        const pic_w = vbox.baseVal.width ;
        const pic_h = vbox.baseVal.height ;
        var r = this.$refs.pic.getBoundingClientRect();
        var x = Math.round(Math.round(event.clientX-r.left) * pic_w/r.width) ;
        var y = Math.round(Math.round(event.clientY-r.top)  * pic_h/r.height);
        this.pin_x = x - 10 ;
        this.pin_y = y - 40 ;
        this.$emit('updatePos', x, y );
    },
  },
}
</script>

マウスをクリックしたときの mousedown 関数では、クリックした座標を補正するためにちょっとややこしいことになっています。

クリックしたときの座標を返すために、this.$emit で updatePos イベントを呼び出します。

Vue 親子コンポーネントの関係は MVVM パターンではない

つまりは、

  • 親から子コンポーネントに渡すときは max_x, max_y のようにプロパティ経由で渡す
  • 子から親コンポーネントに渡すときは updatePos イベントのように this.$emit を通して渡す

のようになるので、MVVM パターンのように対称ではありません。Vue.js が MVVM パターンなのは、ひとつのコンポーネント内での仮想 DOM とのやり取りの話であって、親子コンポーネントは oneway のプロパティ経由なのだ、と。

という訳で、勘違いして迷っていたのでメモ書き。

カテゴリー: 開発 | Vue.js で子コンポーネントから親に値を送る(SVG編) はコメントを受け付けていません

各国の接触確認アプリを探索してみよう

プログラムを書く時に、どのような目的でそのアプリあるいはシステムが使われるのか?を定義しないとあらぬところに労力を掛け過ぎてしまいがちです。労力も時間も有限なのだから(時には無限の場合もあるけど)、注力するべきところに注力して仕上げをしたいところです。

規模の大きいITプロジェクトの場合、なかなか全容を把握するのは難しく、新人の頃は各工程(要件定義や設計工程)の全容が見ないままプロジェクトが終わることも度々あります。
今回の「接触確認アプリ」ですが、実はひと通りITプロジェクトの工程が見えてくる「好例」です。

  • 要件定義&概要設計
  • 必要な API
  • クライアント(Android, iOSアプリ)
  • サーバーサイド

という組み合わせでがあり、しかも各国のソースコードが github に公開されています。実装の違いも含めて、比較できることは珍しいことなので、ざっとでも情報を見ておくと今後に役立つでしょう。という訳で、いくつか見た範囲内での実例をリンクを記述しておきましょう。

接触確認アプリの目的

「接書確認アプリ」の目的は中国などを含めると様々な理由が含められますが、日本版のものは、5/9の有志気概者会議で明確に定義されています。

第1回 接触確認アプリに関する有識者検討会合 開催 | 政府CIOポータル

スマートフォンを活用して、

  1. 自らの行動変容を確認できる、
  2. 自分が感染者と分かったときに、プライバシー保護と本人同意を前提に、濃厚接触者に通知し、濃厚接触者自ら国の新型コロナウイルス感染者等把握・管理支援システム(仮称)に登録できるようにすることで、健康観察への円滑な移行等も期待できる。

これは「定義」なので、これ以上でもこれ以下でもありません。実際にはこの目的から提案依頼書(RFP)や要件定義書に落として「契約」を結ぶわけですが、今回は政府筋なので必要はないでしょう。発注先とは結ばれているとは思うのですが、ここでは問いません。

接触確認(暴露通知)の仕様

接触確認アプリの仕様は有識者会議と厚生労働省が作っています。

第1回 接触確認アプリに関する有識者検討会合 開催 | 政府CIOポータル
第2回 接触確認アプリに関する有識者検討会合 開催 | 政府CIOポータル
接触確認アプリに関する仕様書等の公表 | 政府CIOポータル

「接触確認アプリ及び関連システム仕様書」を見ると、詳細設計まで書いてあるのはさておき、官庁系のものとしては丁寧に書かれているので、これをもとに接触確認アプリとサーバーサイドを作ります。

このシーケンスと内部構造は、Google/Apple の仕様書にも書かれているものなので、そこは一致してますね。陽性患者の登録部分は日本独自なので、ここが作成時の肝となるところです。設計書や実装サイドとしては、この「仕様書」に従ってシステムを構築します。

いくつか興味深い「非機能要件」が定義されてます。

運用に関しては72時間以内(3日間)の復旧程度なので、即時復旧は意味しませんが、それなりに運用体制が必要です。あと、データは「クラウドを利用」しリージョンが「日本国内あること」の記述があります。この辺りも踏まえて受注と運用を検討します。

接触確認アプリの仕様

接触確認アプリは当初は GPS を利用した行動追跡が必須になっていたのですが、Google と Apple が共同で仕様を合わせた Exposure Notification(暴露通知)を使うことで、Bluetooth/BLE を利用したものに置き得られています。

アップルとグーグルが新型コロナ暴露通知アプリのサンプルコードやUI、詳細ポリシーを公開 | TechCrunch Japan

仕様の詳細は、Google/Apple の各サイトを見て貰うことにして、いくつか限定条件があります。

  • 暴露通知アプリは1国1アプリとする
  • 暴露通知とGPSを共有してはいけない
  • 暴露通知を他のアプリにバンドルしてはいけない

1国1アプリの理由は提供元が各国の保険省となるためです。複数出したり、保険省(日本では厚生労働省)以外が出すと審査で落ちるでしょう。
暴露通知とGPSとを共有しないのは、もともとGPSを使わないためなので、それが理由です。逆に言えば、組み合わせて詳細な行動追跡アプリを作ることも可能(昔の彼ログのような)なのですが、これに関しては駄目ですね。
暴露通知は、主に単機能で扱います。おそらくできるだけ許可を少なくする目的と、変な形で企業がバンドルできないようにするためでしょう。実質、厚生労働省以外は出せないので、ここは厚生労働省の要件定義次第になります。1国1アプリとのことなので、自治体が独自に使うことできません。大阪で予定してクーポン付きとかも駄目ですね。

Exposure Notifications: Helping fight COVID-19 – Google

ExposureNotification | Apple Developer Documentation

Googleのほうには、Android アプリの(Kotlin)とサーバーサイド(Java)の例があります。Apple のほうは API リファレンスだけなのですが、各国それなりに Swift で作ったコードがあるので参考にしやすいでしょう。

google/exposure-notifications-android: Exposure Notifications Android Reference Design

試してみるならば、Android のコードを使うのが便利です。Android Studio を使うと簡単にビルドができます。

google/exposure-notifications-server: Exposure Notification Reference Server | Covid-19 Exposure Notifications

サーバーサイドは試してはいませんが docker で動くようです。多分、ローカル環境で Android アプリとサーバーの動作環境を作れると思います。

ドイツの場合

ドイツの接触確認アプリが GitHub で公開されています。

Corona-Warn-App

いくつかのプロジェクトがありますが、Android, iOS, Server を見ておけばよいでしょう。

ドイツの corona-warn-app はドキュメントが綺麗に整備されていて、こんな風にシーケンス図やユースケースが残されています。

cwa-documentation/solution_architecture.md at master · corona-warn-app/cwa-documentation

この手のものが残されていると、ソースコードを修正するときに何が目的なのかが分かりやすいですね。

ちなみに、サーバーサイドは Java + PostgreSQL なので、日本の要件定義に合致しません。

フランスの場合

フランス版は GitLib で公開されています。

StopCovid sources · GitLab

ちょっと、興味深いのは「ROBERT」の文字が入っているところです。

これ、ロベルト・コッホの「ROBERT」なんですね。

ロベルト・コッホ研究所 – Wikipedia
RKI – Startseite

シンガポールの場合

シンガポールですが、コードは公開されてません(多分)。
シンガポール版は、実はハードウェアで実装されていて、不必要になったら捨てることができます。

なぜ新型コロナの追跡はアプリではなくハードウェアで行うべきなのか? – GIGAZINE

子供や高齢者にも利用してもらうことを考えると、こんな形でハードを配布するのがいいんですけどね。今回は、要件定義的に「スマホアプリ」となっているので割愛します。

日本の場合(Covid19Radar)

日本版のベースは Github で公開されています。
基本は Xamrin.Forms(C#) を利用していて、サーバーでは Azure Functions+CosomsDB になります。

Covid-19Radar/Covid19Radar: Open Source / Internationalization/ iOS Android Cross Platform Contact Tracing App by exposure notification framework Xamarin App and Server Side Code

settings.json でアップロード先のURL(API_URL_BASE)やAPI_SECRETが書き替えになっているので、これをローカル環境で動かすようにすれば、ローカルでの実験が可能です。

{
  "appVersion": "APP_VERSION",
  "apiSecret": "API_SECRET",
  "apiUrlBase": "https://API_URL_BASE/api",
  "supportedRegions": "440",
  "blobStorageContainerName": "c19r",
  "androidSafetyNetApiKey": "ANDROID_SAFETYNETKEY",
  "cdnUrlBase": "https://CDN_URL_BASE/",
  "licenseUrl": "https://covid19radarjpnprod.z11.web.core.windows.net/license.html",
  "appStoreUrl": "https://itunes.apple.com/jp/app/id1516764458?mt=8",
  "googlePlayUrl": "https://play.google.com/store/apps/details?id=jp.go.mhlw.covid19radar",
  "supportEmail": "SUPPORT_EMAIL"
}

サーバーサイドの Azure Functions と CosmosDB はローカル版もあるので、ローカル環境を構築して運用試験をすることも可能です。

実際には陽性者番号との突き合わせがサーバーサイドで発生するので、そこは自前で補ってやらないといけないのですが、複数の Android 機を使って一律通知のテスト位はできるかなと思います。

あと電力消費の問題も、複数台の Android に入れて1週間位実験すれば実情がわかるんじゃないでしょうか。

日本の場合(まもりあいJapan)

1国1アプリということで、開発は止まってしまったのですが、まもりあいJapan版も GitHub で公開されています。

まもりあいJapan

もともと Google/Appleの暴露通知ができる前からスタートしているのでGPS利用だった訳ですが、Exposure Notifications版に直して、通知先を Azure Functions に直せばローカルでは動かせるようになるかなと思っています。ちょっと試してみたいところですね。
サーバーサイドは nodejs + Firebase が使われています。

おまけ

接触確認アプリは BLE の近接距離を測っている(実際は電波強度)だけなので、もともと BLE を使うと実装できる話なのです。

接触確認アプリが有効になっているかを調べるアプリです。数十センチ程度まで近づくと素早く点滅します。M5 Atom 用です。https://twitter.com/ksasao/status/1274385507565178885 参照。

だから、本当はアップロードはフリーWiFiなどを利用すれば、ハードウェアでもいいはずなんですよね。とりあえず、皆様がBTを常時有効にし始めるので、データ調査が捗ってしまうという面もあったりなかったり。

すれちがいフレームワークのためのBLEを用いた近接検知機構の実装と評価

BLEはその名の通り低電力で動くのでボタン電池で1年間とか動いたりします。それを利用した近接検知はBLEが出た当初から研究されているので、色々調べると面白いですよ。

カテゴリー: 開発 | 各国の接触確認アプリを探索してみよう はコメントを受け付けていません

プログラマにもわかる SEIR モデルシミュレーション

最初に「SEIRモデル」を知ったのは2月末の頃です。まだダイヤモンド・プリンセス号で新型コロナウィルスが流行っていた頃で東京都にも他の都市もヨーロッパもアメリカもそれほど問題にはしていませんでした。

中国の武漢の状態から「感染力がインフルエンザ並みであるが、クラスター(当時はスーパースプレッダーという用語が使われていました)があり、条件が揃うと十数倍の感染力になってしまう」ことが指摘されていました。

さて、私自身は一介のプログラマなので疫学の専門ではありません。しかし、SEIRモデルが示す通り、数理としての疫学があり専門家会議の西浦教授が「数理」を通して現象を解いているところに共感を覚えます。同時に、コンピュータシミュレーションや統計学の分野として、疫学の「SEIRモデル」を解くことができます。

タイトルに「プログラマにもわかる」という言葉を入れましたが、敢えて言えば「プログラマだから分かる SEIRモデル入門」というところです。SEIRモデル自体は大学の1年ぐらいで学べるような疫学の入門レベルのものです。実際のシミュレーションには数々の現実のパラメータが含まれることになりますが、統計学の主要因分析のように主な要因(原因)特定して、現象を概略を把握しておくとは重要なことです。

というわけで、ちょっとしたプログラム(JS)を利用しながら、SEIRモデルの動作を実験してみましょうというのがこの記事の主旨です。

以下、書くときの効率を重んじて口調が教科書風になりますがご容赦を。

SEIRモデルとは何か?

SEIRモデルというのは、今回のような疫病の場合に患者の状態を4つに分けたモデルである。

  • 未感染者(S:Susceptible)
  • 潜伏期間中(E:Exposed)
  • 発症中(I:Infectious)
  • 免疫獲得者(R:Recoverd)

SEIRモデルは、4つの状態に分かれて、それぞれが直線的に遷移する。

UMLで使われるステートチャート(状態遷移図)と同じだ。4つの状態があって、遷移は、以下の3つの線だけになる。

  • 未感染者から潜伏期間中
  • 潜伏期間中から発症中
  • 発症中から免疫獲得者

未感染者は初期状態で、健康な状態。インフルエンザで言えば、今年の流行に罹っていない人だ。
この未感染者が誰からインフルエンザのウィルスを貰うと、まず潜伏期間中になる。潜伏期間中の場合は、誰にうつす訳ではない期間のことを言う。
潜伏期間が過ぎると、発症期間中に移行して誰かを感染させる状態になる。インフルエンザの場合は、体が動かなくなるので家に必然的に家に籠ることが多い。
発症期間が過ぎると何らかの免疫ができて、免疫獲得者となり回復する。縁起が悪いが、発症期間中に死亡しても免疫獲得者の枠に入れる。免疫を獲得すると誰かにウィルスをうつすこともないし、同時に誰かからウィルスをうつされることもない。

現実には、未感染者から一気に免疫獲得者になる道筋があって、これがワクチンの接種、つまりインフルエンザの予防接種だ。予防接種をすることによって、潜伏期間と発症期間を一気にすっ飛ばして(疫学的には素早く移動するということになるが)、免疫獲得者の枠に入る。

世の中で言われる「集団免疫」という手法は、この「未感染者から一気に免疫獲得者」になることを示している。インフルエンザなどのワクチンが開発された状態(インフルエンザは毎年変わるので、毎年予防接種を受けないといけない)の場合には、途中の発症期間中の死亡率がきわめて低い(実は予防接種でもゼロではない)ので、この「集団免疫」という手段を取る。麻疹とか風疹も同じだ。
しかし、新型コロナウィルスのように決定的なワクチンがない状態で集団免疫の手法を取ると、発症期間に死亡する確率が高くなるため手段として妥当ではない。しかしながら、国内死亡率というのは、単純に新型コロナウィルス由来だけのものではなく、経済的な自殺者も含めるものなので、そのバランスを鑑みて国ごとに方策がを決めることになる。かなりシビアな話だ。

遷移割合を入れる

UMLのステートチャート(状態遷移図)の場合は、100%次の状態に遷移するのだが、疫学の場合は遷移する「割合」が入ってくる。例えば、3000人の未感染者に1人の感染者が含まれていたら、新たに潜伏期間中に遷移する人数はどのくらいだろう?という具合だ。

これを学術的には「摂動論」という言う。解析学が微分/積分を使って数式を解析的に解くのだが、摂動論は状態の差分を逐一計算する。今で言えば、コンピューターシミュレーションのはしりのようなものだ。機械学習で使われる畳み込み計算も摂動論に一種である。一気に解析的に解いてしまうのでなく、状態を1手1手(この場合は1日1日)進めながら最終的な値に近づけていく計算手法である。

数理計算の場合、この手のコンピューターの計算力を使う方法が良く使われる。完全に予測をするためには、「遷移する確率分布」を求める(いわゆるベイズ)を使うことになるが、ここでは単純化させて「遷移する割合」を使う。
ここで、割合をつかっても良い理由がいくつかあって、

  • 都市閉鎖、自粛などで途中の「実効再生産数 Rt」が変化する
  • クラスターの発生などで、遷移確率は正規分布ではない(偏りがある、条件による)

このような理由があって、確率を正確に計算しても意味がない。後述するがダイヤモンド・プリンセス号のような「理想的な状態」の場合には、この SEIR モデルがよくマッチする。逆に言えば、数々の現実の条件が入ったときには、SEIR モデルはマッチしないと言える。

遷移する割合は差分として扱えるので次のように dS/dt としてΔ(差分)として表せる。

それぞれの状態をS,E,I,Rとして表して、n-1番目の状態からn番目の状態に遷移する、というのが次の式の意味になる。n番目というのは、このような1手1手のステップを使うときに良く利用される帰納法的な手法だ。

最初のSとRの式が他と異なるのは、S(未感染者)の場合は入り(プラス)がない。R(免疫獲得者)の場合は出(マイナス)がないことを示している。
つまり、Sは単調減少、Rは単調増加になる。途中のEとIは、増加と現象が前後の状態によって異なってくることが分かる。

このままだと、dS/dt 自体の計算が少しやりづらいので、dS/dt を中心にして式を書き直してみる。

それぞれの差分(dS/dt)には、何らかの割合(パラメータ)が加えられている。これを示すために、

  • S(未感染者)には、β と N(総人数)という割合
  • E(潜伏期間)には、α という割合
  • I(発症期間)には、γ という割合

を掛けることにする。N というのは全体の総人数(S+E+I+R)で全て足した数になる。

SEIR モデルとしては、

  • β : 感染率
  • α : 潜伏率
  • γ : 回復率

という呼び方をする。これは丁度、Compartmental models in epidemiology – Wikipedia の SEIR model の式と同じになる。英語版の SEIR モデルの場合、自然死亡率である μ が含まれているが、ここで扱う SEIR モデルは短期間となるので、μ = 0 とする。

潜伏率(α)や回復率(γ)というのをどのように定義するのかというと、次のように疫学の潜伏期間(lp)と発症期間(ip)によって表せる。

潜伏期間というのは、誰からウィルスを貰って発病(症状が出る)までの期間である。通常の場合、潜伏期間中にはウィルスを外に出さない、つまり誰にもうつさない状態となる。インフルエンザで言えば、2,3日ということになる。
発症期間というのは、症状が出てから回復して免疫を獲得するまでの期間である。モデル的には発症期間中の死亡もここに含まれるので症状が重い場合には、実は発症期間は短く換算される。インフルエンザで言えば1週間ということになる。小学生がインフルエンザに罹って学校に復帰できるまでの期間が1週間なので、これにあたる。

さて、肝心の感染率(β)はどのように計算すればよいだろうか。
今回の新型コロナウィルスで一気に有名なった「基本再生産数 R0」や「実効再生産数 Rt」がこの感染率に関係してくる。
「再生産数」というのは、疫学だけの用語ではなくて「ひとつの個体がどれだけの別の個体をうみだすか」の割合(数)になる。例えば、女性が一生のうちに子どもを生むという再生産数は、2.3程度が理想的と言える。2.0よりちょっと多いのは、途中で病気や事故などで死亡するからだ。現在の少子化により1.6位になっているので、ぐんぐん人口が減ることが分かる。

私自身の専門になるけど、炉物理の分野でも再生産数という意味で「断面積 σ」が使われる。ひとつの中性子が U235 に当たって、いくつの中性子を生み出すかという計算に断面積を使う。この断面積は、核種固有ではなく燃料棒や制御棒などを含めた値でもあり、最近のコンピューターシミュレーションでは、燃料棒、制御棒、一次冷却水を別々の断面積で扱い、再生産数を計算し臨界を保つことになる。ここは余談だが。

  • 「基本」再生産数は環境に依存しない、ウィルス固有の値
  • 「実効」再生産数は環境に依存する割合

という違いがある。「基本」(basic)と「実効」(effective)の違いは他の分野でも出てくる。
さて、ここでは基本再生産数 R0 を使ってみよう。先の英語の SEIR model の中に次のような式がある。

基本再生産数 R0 は、感染率、潜伏率、回復率、死亡率(μ)によって、計算されていることがわかる。ここで、潜伏率(α)と回復率(γ)は既知であり、死亡率(μ)は 0 とするので、次のように式が書き替えられる。

つまり、感染率は、基本生産数に回復率を掛けたものになる。回復率(γ)は、発症期間(ip)の逆数となるので、次の形になる。

ところで、基本再生産数 R0 の定義を思い出してみよう。ひとつの個体がその生存期間中に他の個体を出す割合である。疫学の教科書には、ひとりの感染者が他の感染者(二次感染者)を生み出す数、となっているが、つまりウィルスが生きている間(回復期間に入る前、つまり発症期間I)にどれだけの患者を作るか?という数になる。さきに書いた、少子化の問題になれば発症期間はイコール女性の一生の間という期間になる。

R0 は発症期間に依存していることになる。発症期間に依存している R0 を発症期間(ip)で割るというのが感染率(β)の意味になる。
つまり、感染率というのは単位あたり(この場合は1日あたり)に感染者を増やす割合を示している。
この部分、単位日あたりの感染率と基本再生産数を別に扱っているのは不思議な感じがするが、ともかく同じことを示していることが分かるだろう。多分、「倍加期間」というように、核種の半減期(指数関数)にあわせるための工夫だと思う。

何故、指数関数になるのか

SEIR モデルなど実測の経緯を示すために指数関数(片対数表)が良く使われる。急激に患者数が上がってしまうために指数関数を使っているように見えるが、実はこの SEIR モデル自体が指数関数的に増加するようになっているために、グラフがそうなる。

モデルと現実とが逆転しているようにみえるが、「現実の主要因要素を抜き出して関したものがモデル」であるので、モデルは現実の一部をうまく表せるようになっている。逆に、現実のパラメータがモデルに合わなくなると、そのモデルは捨てなくてはいけない。新しいモデルを適用する必要がある。

これは、パターンランゲージなど活用してプログラムを組むときに似ている。業務モデルを作ってパターンに当てはめたり、パターンをうまく使えるプログラム言語を活用すると、システムは効率よく作れる。しかし、現実をうまく業務モデルに落とし込めなかったり、業務モデルが確定した後に現実のほうが変わってしまった(業務モデル自身によって改善された場合も含む)場合でも、業務モデルから乖離してしまう。
この場合、業務モデルを新しく作り直すのは常だ。
同じ現象は、SEIR モデルにも言える。

指数関数というと自然対数 e の乗数を使ったものが多いが(これは数学的に理由がある。e で指数関数を使うと微分しても、e に戻るので計算が楽になるのだ。ガウス分布がこれにあたる)、なにも疫病の現実が指数関数にマッチしいるのではない。指数関数になる SEIR モデルを使っているから、シミュレーション結果が指数関数的になるのだ。
この因果関係がよくわかっていないと、指数関数にフィッテイングするように現実の感染者数をマッチさせてみたり、S字カーブにフィッティングする、という謎理論が出てくる。いや、謎理論は別にいいのだが、現実に即していないので「科学的」とは言えないだろう。

SEIR モデルの場合、感染して潜伏期間に入る、潜伏してから発症する、発症してから免疫を獲得するというステップがある。これは現実のウィルスも同じである。新型コロナウィルスの場合、潜伏期間中に誰かに移すという報告もあるので発症期間が少しずれるが、これは「誰かに移す期間を発症期間とする」とう定義に合わせて、2日間ほど発症期間に加えてしまえばよい。

あらためて、未感染者 S の差分の式をみていこう。感染率(β)が、未感染者(S)に掛け合わさって感染していくことがわかる。総数(N)で割ってあるのは、これらが単位のない割合となるからだ。

  • 未感染者(S)から一定の割合で感染して潜伏期間中(E)の移る。
  • 潜伏期間から発症期間があるが、潜伏期間中に滞留する。
  • 同じ理由で、発症期間にも滞留する。

これにより、それぞれの「滞留」数が、患者数の累計が指数関数的に増加するという意味あいになる。これは解析的に解いてもいいのだが、現実のパラメーター(実効再生産数 Rt)が日よって異なってくるので解析的に解くのは妥当ではない。日単位でシミュレーション解析をしたほうがよい。

これをシミュレーションしたのが、下の図になる。

オレンジ色の潜伏期間数が、急激に上がっている(線形的ではない)ではことが分かるだろう。特に立ち上がりのところは、指数関数的にあがる。実際は e にフィットする訳ではないが指数関数的に爆発的に増加する(いわゆるオーバーシュート)と言われるのはこの現象である。

SEIRモデルを JavaScript でシミュレートする

統計学的に解析が必要となるので、R を使ったり Python で既存のライブラリを使ったシミュレーションが多いのだが、実は Javascript でも簡単に作れる。

特に、基本再生産数 R0 を色々変えてみたり、潜伏期間や発症期間を変えて現実にフィットするかをチェックするためには、試行錯誤が簡単にできる Javascript がよいだろう。Python で GUI を使って作ることも可能なのだが、Vue.js を使い c3.js のようなグラフツールを作れば、JS でも十分使える。

        seir_eq: function(v,t,alpha,beta,gamma,N) {
            S = v[0]
            E = v[1]
            I = v[2]
            R = v[3]
            ds = - beta * I / N * S             // dS/dt = -βI/N*S
            de = beta * I / N * S - alpha * E   // dE/dt = βI/N*S-αE
            di = alpha * E - gamma * I          // dI/dt = αE - γI
            dr = gamma * I                      // dR/dt = γI 

            return [ds,de,di,dr];
        },

肝は、先の差分の式を JS に直すだけだ。
ここでは、seir_eq という関数を使って差分の部分を計算している。S,E,I,R が配列になっているのは、参考にした python プログラムの名残りである。

これを日単位でステップ実行さると先のようなグラフができる。
シミュレーターは、SEIR モデル シミュレータ で実行ができるようにしてある。

コードは、moonmile/seir-model: SEIR model simulator で公開しているので、詳しくは index.html の中味を見て欲しい。それほど難しいことはやっていないのが分かるだろう。シンプルな SEIR モデル程度ならば、それほど難しくはない。

考察

試しに、以下の数値で計算を実行してみてほしい。

  • 未感染者(S)が3,000人
  • 初期の感染者が5人
  • 基本再生産数 R0 = 10.0

基本再生産数 R0 が 10.0 というのは「理想的なクラスター」の値である。
1か月後には、潜伏期間の患者が1,500人になることが分かる。この数値は、実際にダイヤモンド・プリンセス号の患者数にうまくフィットしている。閉鎖空間における新型コロナウィルスの感染率は非常に高いと言える。逆に言えば、閉鎖空間ではない場合の Rt は小さい。感染させない人が8割ぐらいいるという根拠である。

東京都の都市閉鎖についても、場所によっては同じ現象が起こると考えられる。一概に閉鎖がよい訳ではないのは、経済活動的な理由もあるが、人の行き来を少なくしてしまうために場所によっては「感染に理想的な空間」が生まれやすいという理由でもある。

また、新型コロナウィルスでは以下の恐れがある(まだ事実とは言えない)。

  • 無症状の保菌者がいる。
  • 無症状の保菌者が感染させている可能性が高い

これは今後のニューヨークでの抗体検査、日本での抗体検査によって明らかになっていくことだろう。

新規患者数の遅延について

以下は、隠れ患者を想定した SIR モデルのグラフである。私的に計算したものなので、学術的な正しさは脇に置いておく。

新型コロナウィルスでは、潜伏期間が長い且つ無症状患者が多いと感がられるので「隠れ患者」を試算してみた結果である。

  • 新規患者数を2週間前に遡って、隠れ新規患者数とする
  • 2週間前の隠れ患者数の累計を出し、その時点での新規患者を回復者として取り除く
  • 感染率を推測し、現時点に当てはめ、2週間後の患者数を予測する

という手法を取っている。このグラフは Excel で集計し、以下でダウンロードが可能である。

他にも新規患者数について様々な予測グラフがあるが、以下の点を指摘しておきたい。

  • 現実をモデルにフィッティングさせてはいけない(特にSカーブ)。常に、現実のほうが正しい。
  • 全国一律で予測してはいけない。「クラスター発生」である通り、漠然と発生しているのではなく各地に偏在する、分散がある。よって、最低限、閉鎖する都市単位(東京都、大阪府、札幌市など)で傾向をみる必要がある。

完全な予測は難しいし、実際のところは完全に予測するのは不可能だ。IT 屋なら分かるだろうけど、ウォーターフォール方式でプロジェクトが盲目的に進んでしまえば、そこに待っているのはデスマーチしかない。偶然にもプロジェクトが成功するかもしれないが、それは現実が変わらないという前提で綿密な計画を立てたときに限る。
ならば、その時々にハンドリングをするアジャイル方式で取らねばならない。専門者会議の意見が週単位で変わってしまうのも仕方がない。長期的なスローガンは今回の場合は悪手に陥りやすい。手間暇はかかるが、状況を見据えてその場その場で状況に合わせて短期的な計画を立て、実行して結果を見直すというサイクルが必要だ。

時として、長期的な安心を求めるかもしれないが、それは無駄だ。しかし、暗中模索というわけではない。今回紹介する SEIR モデルのように、ある程度の先行きの予測は「シミュレーション」により可能である。シミュレーション結果がズレていれば、清くパラメータを変えて再びシミュレーションをすればよい。隠れSIRモデルを使って、独自に日々更新しているのはそれの意味もある。

と言うわけで、プログラマならばちょっと手を動かして JS で組んでみてください。Python でも組めるし、これならば C# でもいけるでしょう。Excel VBA でも可能です。是非ためしてみてください。

参考文献

カテゴリー: 開発 | プログラマにもわかる SEIR モデルシミュレーション はコメントを受け付けていません

Vue.js で SJIS コードの CSV をダウンロードさせる

結構探してベストプラクティスっぽいものが無かったので記録として。

サーバーサイドでデータを作る

CSV 形式でデータをダウンロードさせる方法はいくつかあるのですが、Javascript のコードが UTF8 になっているので、これ以外の形式でダウンロードさせるときには結構工夫が必要です。
特に、古めの Excel で直接開こうとするとき SJIS コードが必要になり、これが結構厄介です。

方法としては、

  • サーバーサイドの Laravel 側で SJIS 形式でデータを作って、ブラウザで直接ダウンロード
  • クライアントサイドの Javascript で SJIS 形式データを作って、ダウロードさせる

前者のほうが簡単に見えますが、ブラウザ側で Vue.js を使ったりして、WEB API 経由でデータを取ろうとする途端に厄介毎に巻き込まれます。
ならば、後者のほうで UTF8 で WEB API を呼び出して、Vue.js 側で SJIS 変換すればよい、と思いつつもピッタリのサイトが見つかりませんでした。

結論から言えば、encoding.js を使って UTF8 から SJIS に変換できます。

Vue.js でダウンロードさせる

手順は2つになります。

  • Web API 経由で、JSON 形式でデータを取得
  • 取得したデータを encoding.js でコンバートしてダウンロード

Web API 経由でデータを取得したほうが、Vue.js などの SPA では必須になります。API の呼び出しは axios を使っています。

axios.post(url, { params: params })
.then(function(response){
    var head = "...";
    var filename = "sample.csv";
    this.downloadCsv(head, response.data, filename );
}

適当なパラメーターをサーバーサイドへ POST します。
取得できるデータは JSON 形式としておきます。

CDN で encoding.js をインクルードして、downloadCsv 関数と downloadCsvBlob 関数を定義します。
downloadCsvBlob 関数は、Edge と Chrome の両方を動作させるための共通関数です。この関数は、【JavaScript】Axiosを使ってCSVをダウンロードしたい | ゆとって生きたい。 から持ってきています。

<script src="https://cdnjs.cloudflare.com/ajax/libs/encoding-japanese/1.0.30/encoding.js"></script>

function downloadCsv( head, items, filename ) {
    var csv = "";
    items.forEach(function(it) {
        csv += Object.keys(it).map(function(item){
        return it[item]
        }).join(',') + "\n";
    }.bind(csv))
    csv = head + "\n" + csv ;
    csv = new Encoding.stringToCode(csv)
    csv = Encoding.convert( csv, 'SJIS');
    csv = new Uint8Array( csv );
    this.downloadCsvBlob( csv, filename );
}

function downloadCsvBlob(blob, fileName) {
    if (window.navigator.msSaveOrOpenBlob) {
        // for IE,Edge
        var blob = new Blob([blob], { "type" : "text/csv" });    
        window.navigator.msSaveOrOpenBlob(blob, fileName);
    } else {
        // for chrome, firefox
        const url = URL.createObjectURL(new Blob([blob], {type: 'text/csv'}));
        const linkEl = document.createElement('a');
        linkEl.href = url;
        linkEl.setAttribute('download', fileName);
        document.body.appendChild(linkEl);
        linkEl.click();
        URL.revokeObjectURL(url);
        linkEl.parentNode.removeChild(linkEl);
    }
}

コード変換の肝は、以下の部分ですね。
いくつか試してみたのですが、この順番でないとうまく動かないようです。変数 csv を省略しようとする動かなくなったりするのは、内部データの使い方の問題だと思いますが、ひとまず、以下の形で UTF8 から SJIS に変換ができるようになります。

    csv = new Encoding.stringToCode(csv)
    csv = Encoding.convert( csv, 'SJIS');
    csv = new Uint8Array( csv );

これでブラウザ上のボタンをクリックすると SJIS 形式のデータをダウンロードさせることができます。当然のことながら、Web API を叩かないで、ブラウザ上だけで計算してデータをダウンロードさせることも可能です。

参考先

カテゴリー: 開発 | Vue.js で SJIS コードの CSV をダウンロードさせる はコメントを受け付けていません

学童でプログラミング教室を開いた3年間の話

少しプライベートな話になってしまうが記録として残しておこう。
ひょっとすると個人が特定できてしまうかもしれないが、特定せずに「マックで女子高生が話していた」程度に捉えてほしい。

プログラミング教室の事始め

板橋区のとある学童でプログラミング教室を隔週で開催している。板橋区は学童を「あいキッズ」という形で民間に委託しており、それを受けた民間の会社やNPO法人が委託した予算内でいろいろな試みをしている。私が関わっている「放課後NPOアフタースクール」もそのひとつで、あいキッズの中で、いくつかの学習企画を立てている。

その昔、学童といえば共働きの鍵っ子や、母子家庭の子が預けられる子というイメージが強かったのだが、最近はそうでもない。共働きの家庭がクラスの半分に及び、学校が終わってから親が家に変えるまでの数時間の間、小学生は学童ですごすことはごく普通な日常となっている。だいたい、学年で約半分の子は学童に通っていたりする。

学童は保育園と同じ様に保育の場ではあるもの、「あいキッズ」が民間に委託されている特性から、

  • 学童で何かを学習させる
  • 学童で何かを運動させる
  • 学童で何かを食べさせる

の用途もあいキッズに含まれている。学習や運動は、職員が行う折り紙教室やレゴ作り、サッカー、鬼ごっこなど、毎日細かなイベントが組まれている。お誕生日会や宿題の面倒などもあいキッズの役割である。

そして「放課後NPOアフタースクール」の場合は、職員が行うイベントとは他に、「市民先生」と呼ばれる板橋区から通えるような距離の民間の先生を呼んできて、1時間程児童に何かを教えるという学習スタイルが組まれている。

その中でプログラミング教室の打診を受けたのは、3年間のことで、まだ小学校でプログラム研修が本格化する前だった。しかし、既に Scratch などの子供向けのプログラミング教材は整いつつあった。最初は週1という話であったのだが、フリーランスにとって週1で時間を取られるのは難しく、隔週にしてもらった。これは金額的に安いというのも最初の理由である。

プログラミング教室の準備

子供向けプログラミング教室と言えば、月1万円の月謝を取るところも少なくはない。教材に至っては、LEGO などを使って5,6万円するものもある。一般向け(いわば富裕層向け)にならば、その値段でも殺到するかもしれないが、学童でのプログラミング教室では到底無理である。同時に、あいキッズに委託している金額もそう多くはないので、プログラミング教室だけにお金を掛けるわけにはいかない。

よって、

  • できるだけ、低予算(ノートパソコンのみ)で行う
  • 印刷教材は、あいキッズ側で受け持って貰う
  • 毎回のパソコンのセッティングはあいキッズ側で行って貰う

という条件で始めた。大抵のプログラミング教室では据え置きのPCかノートPCが使われる。専用の場所が確保されるため、PC や教材を据え置くことができる。しかし、今回のあいキッズの場合は、共同の教室での開催となるので、毎回ノートPCを設置/片付けが必要になる。これを毎回、私ひとりで行うのは大変なので、あいキッズの職員にやってもらう(結果、あいキッズのアルバイトがやることなった)。

同時に、児童分の印刷教材を家でプリントアウトすると大変なので、これもあいキッズ側で行って貰う。印刷教材も低予算にするために、毎回数枚のカラープリントで済ませている。実は、印刷教材はなくてもよかったのだが、プログラミング教室を始めるときに

  • 参加する児童から、いくらかの月謝(300円程度)を取る

ことに決めた。これは、プログラミングが人気になりつつあったと、パソコンを使わざるを得ないので、人数を制限したかったからだ。月謝はあまり意味はない。学童の場合には、月1,000円でもキツイ家庭があるので、300円程度と設定している。
このため、300円分の何かのお返しをしようというわけで、印刷教材をつくることにした。実際の時給は度外視である。

パソコンを毎回移動しないといけないために、ノートPCを購入することにした。
新品では難しいので、中古の Windows 7 を探し貰うことにした。ノート PC の選定や購入に関しては、本来ならば私がやるところなのだが、それだと「お金」が掛かってしまうので、これもあいキッズ側でやってもらうことにした。
ノートPC は機種はばらばらだが8台ぐらいまで用意した。これを2人交代で使うとして、16人までの教室として設定した。

この「費用」にいちいち細かく言及するのは、このプログラミング教室を「ボランティア」でやる気はなかったということだ。いや、実質ボランティアかもしれないが、完全な無償の場合には長く続かない(実際、市民先生で残っているのは私ひとりだけだそうだ)。いくらかの報酬を貰い、それに見合うだけの仕事というスタイルとることに決めているし、実際よかったと思っている。

金額的には月1万円(源泉徴収、振込手数料などが抜かれるので、実際は9,500円位になる)で請けおっている。
仕事としては到底割に合わないことが理解できるだろう。

初年度の意気込み

ノートPC も準備してもらい、ある程度の印刷教材を用意して、4月からスタートした。
教えるプログラム言語は Scratch と決めておき、いずれは Python か Arduino で電子工作、という妄想を抱いていた。

そう、妄想は第1回目から打ち砕かれてしまうのだ、

  • 教室にいる人数が予定よりも多い。どうやら水増しで 20 人位になっていた。
  • 白板みたいなものがないので、解説ができない、話を効聞かない
  • みな一斉に同じことをやる、ということが出来ない、待てない

一応、小学3年より上という制限を付けたのだが、いやいや無理、社会人の社員研修とは全く違う。これは自分の子を考えて分かったことで、20 人位あつまると教室は騒がしくなんともできない状態になってしまう。

3、4回ぐらいまでは自前な教材に沿ってぽちぽちやっていたのだが、諦めた。これは好き勝手やって貰うほうがよい。教えるという教師的な目線ではなくて、プログラミング教室に居るおっちゃんのイメージに切り替えた。

毎回、印刷教材を5分だけ解説する。その後は自由に Scratch で遊ぶ。何か聞かれたら答える。私自身も暇なときは Scratch で遊ぶ。そんな適当な雰囲気にすることにした。

学習障害っぽい子がいる

学習障害の場合、心理学的言えば、多動性、注意欠陥、アスペルガーと色々な症候群があるのだが、ちょっと発達障害っぽい子がいた。実際、板橋区では小学校の中に発達障害のクラスを別途もうけているところが多い。「5組」と言ったりする。

その昔、今で言えば発達障害っぽい子がクラスに1人か2人はいたものだ。授業中に無闇に騒いでみたり、椅子に座っていられなかったり、大声を出したり、逆に緘黙症だったりと。今でいえば、軽い発達障害なんだけど、みな同じクラスにてちょっと変な子だった。それでも友達だった。

そうい意味で、学童でもちょっと変な子が混じっていて、その子も友達の輪に入っている。
仮に「H君」と呼ぼう。正確な症例はあえて聞かなかったが、5組な子なそうで、あいキッズの職員も手を焼いているとのこと。H君は、当時小4だった。

  • ひとりでノートPCを占有できないと騒ぎ立てる
  • ペアで組んでいる子に順番を渡さない
  • プログラムが上手くいかないと拗ねる
  • 拗ねると、ロッカーとか入って中味をぐちゃぐちゃにする

もう、まったくまともに扱ったら胃が痛くなる感じである。
ちなみに、あいキッズの職員に「プログラミング教室をやっている間は、H君がいないから楽なんですよ」と言われてしまいました。その分、プログラミング教室は大変なんですけどね。まあ、いいんですが。

H君の対応に対して、あるとき気付いた、

  • うちの小1の子と同じじゃないか。

なるほど、自分の息子(当時、小1)に Scratch を教えようとしたら全然覚えなくて、癇癪をおこして、遊び出して、云うことはきかない。なんだ、小1の子と同じだ、ことに気付いた。

それ以来、H君対応は私にとって楽になって、小1だと思うと全然苦にならなくなった。発達障害というのはそういうものだし、後にH君は落ち着くことなる。小5の後半位から落ち着き始めて、小6だと全然大丈夫になっていった。今回の新型コロナウィルス騒ぎで、卒業のお祝いを言えなかったのが残念なのだが、彼に対しては少しだけスキンシップを行うにした。頭をなぜたり、何かやっているときに肩を叩いてあげたり、逆に好きなことをやっているときはほおっておいたり。他の子のガードもあり、彼だけに私のノートPCを使っても貰っていた。子どもに特性にもよるけれども、H君にはちょうどそれがよかったらしい。いや、私はそう思う。

学生ボランティアの存在

プログラミング教室を始めるときに、ひとつ注文があって「誰かボランティアを付けてください」ということにした。これは、16人(結果的に20人位だが)を1人で相手をするのは無理だと思ったからだ。実際無理だ。本格的なプログラミング教室ならば、児童ひとりのノートPCを一台ずつ付けて、8人位が限界だ。

最初は、交代であいキッズの職員が担当していたのだが、パソコンに不慣れというか、どうも「教育」に不慣れなところがあって、なにかと疲れた。

3か月ぐらい経った頃だったが、大東文化大の学生がボランティア(アルバイト?)で常駐することになった。後から分かったことだが、教育学科で、その後は大学院に行った学生である。今年の3月に大学院を修士で卒業している。
この学生が、子ども達に無茶な注文(あれこれ暴れる、あれこれ喧嘩する、あれこれ走り回る)に対応してくれたので非常に助かった。彼女がいなかったらプログラミング教室は廻らなかっただろう。今年から居ないので、どうなるかわからないけど、ノウハウ的に3年間溜まったのでなんとかなると思う、まだ再開していないけど。

学童は半分しつけもやるので、「子どもを叱る」というのが頻発する。大人が真剣に子どもを叱ると、子どものほうも真剣に泣いたり反省したりする。へらへらとあしらうと、へらへらと返してくる。
私は子どもたちを「教育」する立場でははない(教育免許も持ってないし、そういう学習も受け手は無い)ので、其処のは踏み込まないこにした。あくまでプログラマ的な立場や、心理学的な立場(ユング心理学とか)での立場を守って、踏み込み過ぎないようにしている。

Scratch を自由に扱う

最初は、Pyton や Arduino と思っていたのだが、横並びの授業はできないことがわかったので、Scratch を自由に扱うことにしている。パソコンを使っているので、

  • キーボードを自分で打つ
  • Youtube は見ない

というルールにしている。Youtube を見ると回線がパンクするという理由もある。
Scratch には色々なゲームがアップされているので、これは自由に検索して使ってよいことにしている。一応、毎回課題やプリントは配っているのだが、無視するのは自由(苦笑)。課題をこなすのも自由である。
女の子のほうが比較的まじめに課題シートをやるのだが、必ずしもそうではない。ずっと絵を書いている子もいれば、やっぱり飽きて来なくなってしまった子もいる。
別の Scratch 教室に通っている子は、どんどん課題をこなしたり自由に作品を作ったりする。作品を作っている横で3人で交代で Scratch のゲームをしている子たちもいる。いろいろ雑多だ。

キーボードを自分で打つというルールは、キーボードに慣れるためだ。他に、ゲームをやりたいときに検索ボックスに自分でローマ字で打つことにしている。横でアルファベットを教えてあげるけど、キーボードを打つのは子ど自身だ。間違っても根気よくアルファベットを言う。

何か友達がゲームをしていて、あのゲームがしたい、と私(先生)に聞いてきたときは、友達に直接聞いて、ということにしている。学童の場合、学年はばらばらだけど皆顔見知りなのでそれなりに「友達」だったりする。なので、聞けば教えてくれる。最初は友達じゃなくても教えてくれる。これは重要だ。

プログラミング教室内でも、いくつかのブームがあって、

  • ファミスタ野球ブーム
  • ペーパーマインクラフトブーム

なんてのがある。試しに2台ほど、本物の?マインクラフトを入れたのだが、そのうちにペーパーマインクラフトのほうに戻ってしまうのは不思議なところだ。

いまでは「支える」ことに集中している。

来なくなった子も何人かいる

3年間やっているが、途中で来なくなった子もいる。プログラミング教室は他の部活動?(サッカーとか)に比べると抜ける子はそう多くないそうなのだが、それでもいる。

  • タイピングソフトを所望していたが、そのうち来なくなった
  • 絵を書いていたが、来なくなった
  • 塾が忙しくなって来なくなった

私はパソコンが好きでやっているのだけど、向かない/面白くない人もいるのも分かる。なんで、プログラミング自体が向かない、ゲーム自体がたいして面白くないという子や人もいる。そういう場合は、サッカーとか書道とか別の活動をやればよいのだ。

小学校で本格的にプログラミング授業がはじまると、このあたりが問題になってくると思う。プログラミングだけが特別とは言わないし、逆にプログラミングをやらなくてもいいと思う場合も多い。プログラミングをやらないから「論理思考ができない」とは思わない。ただ、いくつかのプログラミング手法が、論理思考の手助けになることは確かなことだ。

あと、小6になると中学受験で来なくなった子が多くなる。その子は往々にして Scratch ができる子なので、もったいないといえばもったいないけど、限られた時間は限られた時間になるので、そこは中学受験に割り振られれうのだろう。

今後の方針

4月から4年目となるのだが、この新型コロナウィルスのおかげでまだスタートしていない。
一般的なプログラミング教室とは異なり、学童で行うプログラミング教室は、学童そして「子どもを安全に預かる」ことが最優先になる。なので、無理矢理教えようとしないのが長く続くコツとなる。教えたりするのはおまけだったりする。

が、さすがに3年目のペーパーマインクラフトブームは長く続きすぎたので、もうちょっと教材を考えようかなと思っている。

ちなみ、ペーパーマインクラフトは、マイクラ remix-2-2 on Scratch で遊べる。いやー、よくできているので、長々と時間が潰せます(苦笑)。お試しあれ。

カテゴリー: 開発 | 2件のコメント

隠れSIRモデルを試作して、東京都の感染者を予測する

実地な値を使って検証するのは避けていたのだけど、現在の状況を見て解説することにします。

隠れ感染者が一定量いる

現在、東京都の新規患者が100人/日を超えています。指数関数的増加(差分が1.0以上となるので複利的に増加する)が理解できていれば、次の増加が100人ではなく、それ以上にとなることは容易に推測ができます。

ただし、単純に「指数関数的増加」とはいえ、いったいどのくらいの増加数なのか?を予測しておくことは大切なことです。

新型コロナウィルスの厄介な点は、

  • 感染率がインフルエンザと同程度に低い(麻疹などよりもかなり低い)
  • 潜伏期間が2週間と長い
  • 感染しても無症状な人が一定量(8割程度)いるらしい

ところです。このため、麻疹のように一気に拡散&重症化する場合は、発病後すぐに家や病院にこもるために、あまり外側に拡散されることがないとありません。潜伏期間が長いと、いつ罹患したのかがやっかいです。クラスター対策班がクラスターの追跡ができなくなってきているのはこのためです。同時に、感染しても無症状であれば、感染したまま普通の生活を送ることになります。無症状であるために、病院に行くこともなく家で休むこともなく PCR 検査や検体検査を受けることしません。このため、本人の意図とは関係なく、ウィルスが蔓延してしまいます。

進化論的に言えば、これが新型コロナウィルスの「生存戦略」ということになります。

隠れ感染者を推測する

潜伏期間中であったり無症状な感染者を「隠れ感染者」と呼びましょう。この「隠れ感染者」は数値の上ではでてきません。毎日、定期的に発表される新規感染者とは別に扱う必要があります。

では、どうやって「隠れ感染者」を推測するのでしょうか?

感染の因果関係として、

  • 新規感染者が発生する
  • 新規感染者は潜伏期間(2週間)前に何処かで感染している
  • 新規感染者を感染させた元の感染者がいる

ことが条件になります。実際は海外から帰国した感染者もいるので誤差でるのですが、おおまかな東京都の感染者の状況はこの感染の因果関係にあてはめられます。

これを SIR モデルに当てはめてみます。

新規感染者 dIt/dt を発症した日(14日前)に「隠れ新規感染者」として設定します。
隠れ新規を累積して、隠れ累積患者数を計算します。このとき、その日の新規患者を取り除きます。これは新規患者としてカウントして隔離されるためです。
隠れ累計と隠れ新規の関係から、感染率が計算できます。感染率 β は、もとの SIR モデルの感染率や基本再生産数とは異なる定義となるため、このモデル独自のものです。もう少し式変形をして、基本再生産数と比較できるようにしていきます。

既に新規感染者が計測されている日から14日さかのぼったところの感染率βが算出できます。
これをもとにして、予測感染率β*を推測します。
一般に、何らかの環境の操作(外出の自粛など)が行われない限り、予測感染率は前日の感染率を踏襲するはずです。

この予測感染率 β* から、予測隠れ新規患者数を計算します。
予測隠れ新規の数は、14日後に予測新規患者数となり、新規の感染者数が予測できるというモデルです。

Excel で予測する

これを Excel 使って計算してみましょう。

東京都 新型コロナウイルス陽性患者発表詳細 – 東京都_新型コロナウイルス陽性患者発表詳細 – 東京都オープンデータカタログサイト

4/6 付けの東京都の患者データを使って予測してみます。

青いセルが予測の部分です。「患者累積」と「隠れ累積」が大きくずれているのは、新規患者をうまく拾えていないのと、発症までの潜伏期間が2週間と長いために潜在的な患者が滞留してしまっているのが理由です。
単純計算いけば、新規患者数100人×14日間 = 1400人以上は滞留していることになります。

2週間ほどさかのぼって、感染率の予測値を入れます

  • 0.160 は、3/22 までの大まかな平均値です。4/6 付けで新規患者数が83人と若干減っているので実効感染率が一時的に下がっていますが、いまのところ外れ値として処理しています。継続して下がっていれば、この感染率をさげます。
  • 4/4 からの 0.080 は土日で自粛が効果的に行われたとみなして下げてあります。

感染率の絶対値は、特に意味はありません。いずれ実効再生産数のような割り出しが必要になります。
ここでは、相対的な値となっています。

それぞれの値をプロットしたの下の図です。

予測患者数累積が非常にあがっていることが分かります。

対数軸にプロットしたとき、患者累積と隠れ累積の差が大きくならないように注意が必要です。

現時点でプロットしてみたところ

  • 4/11 頃に患者数の累積が2,000人を超える
  • 新規患者は500人のピークがある
  • 4/4 の自粛は、4/18 頃に効果が現れるので、それまでは同じペースで自粛が必要となる
  • 増え方は急ではあるが、きちんと自粛をすれば、患者数が1万人を超えるのは4月末以降の予想となる

問題は、新規患者数が500人/日程度発生するが、医療崩壊を起さずにこれを受け入れらるか?ということです。
このモデルで計算したときには、「隠れ累積患者数」が4,000人程度いることと2週間遅れで発病することをあわせると、かなり厳しい状態と言えます。ただし、軽症/重症者が混在するので、この比率が問題になります。

データ

Excel データは以下でダウンロードが可能です。

参考文献

カテゴリー: 開発 | 隠れSIRモデルを試作して、東京都の感染者を予測する はコメントを受け付けていません

感染症数理モデルでシミュレーションする(その2)

3週間ほど経ち状況が変わってきていますが、感染症数理モデルの話を引き続き

SEIR モデルを JS で書く

新型コロナウィルスの専門家委員会から「基本再生産数 R0」や「実効再生産数 R」という言葉でてきます。他にも専門用語がでてきますが、これらを一般的に解説するかどうかの前に「その専門分野ではどのような使い方をされているのか?」を正確に把握するする必要があります。まあ、それが「専門分野」なのですから、同じ言葉を使っていても専門分野ごとに多少の違いはあるし、逆に同じ用語を他の分野でも同じ様に使っている場合もあります。どちらにせよ、その分野での専門用語の把握が必要です。

そんな訳で、2週間ほど前になりますが、SEIR モデルを Javascript で書きました。元の Python からそのまま JS に書き替えて、そのあと R0 を取り込めるようにしています。ちょっと、R0 の意味が分からなくてコードが間違えていたのですが、先日直したので、ここにコードの解説をしておきます。

SEIR モデルの数式

実は、Wikipedia を見ると、日本語のものと英語のものと SEIR モデルの式が若干違います。感染率 β の部分で総数 N が出てくるのですが、日本語の場合は、β(N) のように総数に依存するのですが、英語版のほうは β は N から独立しています。
感染率が総数 N に依存するのは変な感じがするので、英語版のほうを使うことにしました。

The SEIR model

出生率 Λ と死亡率 μ は短期間なので、0 とみなせます。これを書き方えたのが下の式になります。

  • N: 全体の総数
  • S: 未感染者
  • E: 潜伏期間中
  • I: 発症者
  • R: 回復者/免疫を得た人

になります。イギリスで発表された全国民が感染して免疫を得るという方式は、この SEIR モデルの、R になるということです。SEIR モデルの場合、感染して潜伏→発症→回復という順番になるのですが、かならず回復するというモデルなので、今回のような新型コロナウィルスのように発症後に死亡するケースがある場合には、致死率が2%程度だとしても相当の死亡者が予想されます。

即効撤回されたのは、良い判断でした。

この式の中でのパラメータとして、以下の3つがあります。

  • α: 潜伏率
  • β: 感染率
  • γ: 回復率

これらは、潜伏期間 lp と発症期間 ip で書き替えることができます。

  • α = 1/lp : 潜伏帰還の逆数
  • γ = 1/ip : 発症期間の逆数

肝心の感染率 β の値ですが、基本再生産数 R0 との関係が以下のようになります。

死亡率 μ が 0 となるので、R0 は 感染率 β と 回復率 γ との比になります。
式の中では感染率 β を使うことになるので、次のように書き替えます。

この式から分かることは、

  • 感染率 β と回復率 γ が等しいときに R0 が 1.0 になる
  • 感染率 β が回復率 γ が大きいときに R0 は 1 を超えて、感染が拡大する
  • 感染率 β が回復率 γ が小さいときに R0 は 1 より小さくなり、感染は縮小していく

という現象になります。これが専門家会議の言う「基本再生産数 R0」の意味です。
なので、できるだけ R0 < 1.0 の状態を保ちながら、ピークを低く抑えましょうという対策の根拠になります。

SEIR モデルを JS に直す

この数式を Javascript に直すと次のようになります。
v のところが配列になっているのは Python コードの名残りです。

function seir_eq(v,t,alpha,beta,gamma,N) {

    S = v[0]
    E = v[1]
    I = v[2]
    R = v[3]
    ds = - beta * I / N * S             // dS/dt = -βI/N*S
    de = beta * I / N * S - alpha * E   // dE/dt = βI/N*S-αE
    di = alpha * E - gamma * I          // dI/dt = αE - γI
    dr = gamma * I                      // dR/dt = γI 

    return [ds,de,di,dr];
}

これを100日間繰り返すのが次の calc 関数です。
SEIR モデルに、試しに隔離率(T) を埋め込んでみて、発症者の一定割合を病院に隔離するというシミュレーションも入れておきます。
あと、S,E,I,R の数はマイナス値にはならないので繰り返し計算をするときに補正します。

function calc(state,alpha,beta,gamma) {
    var t_max = 100 ;
    var dt = 1 ;
    lst = []
    var N = Sinit + Einit + Iinit + Rinit
    console.log( state );

    for ( var i=0; i<t_max; i++ ) {

        var d =  seir_eq( state, i, alpha,beta,gamma, N )
        var Si = state[0]+d[0]
        var Ei = state[1]+d[1]
        var Ii = state[2]+d[2]
        var Ri = state[3]+d[3]
        // 感染者を発見して隔離する
        dx = Ii * T
        Ii = Ii - dx
        Ri = Ri + dx // 免疫者に加算

        // マイナス値を調節する
        if ( Si < 0 ) {
            Ei = Ei + Si; Si = 0;
        }
        if ( Ei < 0 ) {
            Ii = Ii + Ei; Ei = 0;
        }
        if ( Ii < 0 ) {
            Ri = Ri + Ii; Ii = 0;
        }

        state = [ Si, Ei, Ii, Ri ]
        // console.log( state );
        lst.push( state );
    }
}

グラフ化する

これに vue.js とグラフツールの c3.js を追加したが次の形です。

SEIR モデル シミュレータ

このグラフは基本再生産数 R0 が 10 のときのものです。10 というのは非常に多い値なのですが、いわゆるダイヤモンドプリンセス号のような閉鎖空間のクラスター(患者集団)の場合がこれに相当します。時間を追うごとに加速度的に感染者が広がります。

注目して欲しいのは、以下の3点です。

  • 発病者のピークは感染スタートよりも左にある
  • 潜伏期間中(黄色)のピークは、発病者(緑)のピークよりも左にある
  • 最終的にほとんどの人が回復者(赤)となる

このグラフでは潜伏期間2週間で考えているので、発病者(緑)のピークは感染スタートから2週間以降になります。感染スタート時にすべての人が感染する訳ではないので、幅を取って2週間から3週間というところです。これが、現在の花見や春休みのシーズがスタートとなったときに、2,3週間後に流行するのでは?という理由です。

当然のことながら、発病期間のピークの前に、潜伏期間のピークがあります。SEIR モデルは潜伏期間中には罹患させないモデルになっていますが、新型コロナウィルスは潜伏期間中も罹患させる可能性がある、または無症状の場合も感染さえる場合があるという「可能性」が指摘されています。
このため、潜伏期間中の人が活発に活動すると、その後に発症のピークが出るかもしれません。

SEIR モデルの場合は、最終的にほとんどの人が感染して回復します。というか、無限に計算を続ければ100%の人が感染し回復して終わるようになっています。
先に書きましたが、SEIR モデルの場合は発症→回復となるので、死亡者がいません。実際の新型コロナウィルスでは死亡するので、最終的な回復者数×死亡率が全体の死亡数になってしまいます。
このため、無防備な策は誤りです。

基本再生産数 R0 を 1 に近づける

基本再生産数 R0(あるいは実効再生産数 R)は、感染率 β と回復率 γ との比率になります。

R0 を 1.0 に近づける(1.0以下にする)というのはどういうことかというと、数式からは以下の2つの手段があります。

  • 感染率 β を下げる
  • 回復率 γ を上げる

感染率を低く押さえるというとは、いわゆる「誰かに罹患させる確率を減らす」ということです。移動を減らすとか、3つの条件が重なる場所を作らないあるいは行かない、というのがその手段です。
もうひとつは、回復率を上げるというのは、病床を増やす、体力を蓄えておいて軽症で済ませるなどがあります。

基本再生産数 R0 を 10,5,2 と変化させてみましょう。



発症者(緑)のピークが限りなく右にずれることがわかります。
つまり「ピークを平らにする」ということは、基本再生産数 R0 を限りなく1.0に近づける(あるいは1.0以下にする)ことであり、同時に感染者のピークを限りなく未来に追いやることができます。

そして、新型コロナウィルスに対しての有効なワクチンや予防接種ができるまでの時間稼ぎができます。

基本再生産数 R0 が 1.0 のときは?

実験として、基本再生産数 R0 が 1.0 のときはどうなるのかを見ていきましょう。

感染者数(緑)のグラフが非常に平らになりほとんど直線になります。
罹患した後に回復した人も全然増えていないように見えますが、実際は回復者が100日目で25人になります。

これから分かることは、

  • 感染がまったく広がらない訳ではない。わずかずつ感染は広がる
  • しかし、対処可能な感染者数となる

ことを示しています。確率的に感染者は発生するのですから、何処から感染するかは分かりません。しかし、感染して発病したとしても、十分に治療可能な病床が用意できる、ことを意味しています。

github

コードは github moonmile/seir-model: SEIR model simulator に公開しています。

おまけ 実効再生産数 R とは?

ここからは私的なメモです。

基本再生産数 R0 と実効再生産数 R の違いを記述しておきます。
専門家会議では、最初の NHK での解説で「基本再生産数 R0」という用語を使っていました。再生産数という言葉から一般的に「基本再生産数を 1.0 以下に抑えれば、なにかグラフが平たくなる」ことが理解できたと思います。

このグラフが出たのが 2/24 の NHK です。

この感染症の流行モデルの元ネタは当時の NHK の解説では何処にも出てきませんが、専門家会議に北大の西浦教授がいることと、西浦教授の2017の論文から、SEIR モデルがベースであることがわかります。

そのあと、色々な批判もあったのち、海外でもこの「平たくする」図がでてくるようになりました。

Flatten the curve | These guidelines are intended to help Flatten the Curve with the COVID19 outbreak, to help limit spread and reduce the load on hospitals and other healthcare.

ただし、これらの「ピークを平たくする」というイメージには、なんら数式的な根拠が示されていません。SEIR モデルを理解しているのか、それとも読者が SEIR モデルを理解できないと思っているのかは定かではないのですが、肝心の数理モデルはでてきません。twitter で #FlattenTheCurve というハッシュタグで拡散されつつありありますが、根拠は示されてません。

結果的に、発症時のピークを下げる=医療崩壊を防ぐことになるので、効果は同じだと思われるのですが、私が懸念するのは、

  • ピークを下げた時期の概算が示されていない
  • ピークを過ぎた後の発症数の裾野が、いつまで続くのか示されていない

のが問題だと思っています。

また、新型コロナウィルスの特徴として無症状な人が多いことから、新しい感染者を0人にすることはほぼ不可能です。天然痘のように全ての人にワクチンをうって撲滅しない限り 0 にはなりません。

同時に「クラスター対策」として言われるように、基本再生産数 R0 は状況によって大きく変わります。全体の平均として 1.0 以下におさえると、マクロの全体として発症者が押さえられるのは確かなのですが、個々の感染しやすい場所、あるいは感染しやすい季節や場所などのミクロの視点で見れば、基本再生産数は大きく変わり分布を持ちます。

このため「実効再生産数 R」として「実効」を付けて、基本再生産数 R0 になんらかのパラメータを掛けます。いわゆる事前条件を加味します。
どうやら疫学の場合の実効再生産数は、二次感染の平均を考慮するようです。このあたり、最近の専門家会議での解説をみると「実効再生産数 R」と区別されていることがわかります。

詳細は2020年3月19日付けの提言を読むとわかります。

「感染拡大地域では自粛検討を」専門家会議が提言【全文】

  • 基本再生産数(R0:すべての者が感受性を有する集団において1人の感染者が生み出した二次感染者数の平均値)
  • 実効再生産数(感染症の流行が進行中の集団のある時刻における、1人の感染者が生み出した二次感染者数の平均値)

基本再生産数は、一般的な都市での感染率 → 都市機能を保ったまま低く保つ必要がある(一般的なイベントの開催、通勤や通学の正常化など)
実効再生産数は、感染してしまった都市での感染率 → 緊急で下げる必要がある(北海道の緊急事態宣言、ロックダウン、都市封鎖)

のように区別します。

カテゴリー: 開発 | 感染症数理モデルでシミュレーションする(その2) はコメントを受け付けていません