SIR 感染症モデルのシミュレーター

// この投稿は ライフゲーム Advent Calendar 2017 の 19 日目の記事です。

2017年11月11日の稲葉寿先生還暦記念祝賀研究集会「数理人口学・数理疫学・構造化個体群モデル」で講演してきました。
そのときに発表したシミュレーション ツールを紹介します。

数理生物学で「SIR 感染症モデル」という数理モデルがあり、

S(t): 未感染者の人口
I(t): 感染者の人口
R(t): 回復者の人口

としたとき、

S’ = θR – βSI
I’ = βSI – γI
R’ = γI – θR

のような微分方程式系で表されるものを指します。ここで、各定数は

β: 感染率
γ: 回復率
θ: 免疫喪失率

を表します (このように免疫喪失を考慮する場合は SIRS ともいう)。

このモデルをもとに、感染症の伝播を視覚的に表現するセルオートマトンを WPF で作成しました。
各ファイルは GitHub にあります。

EpidemicSimulator.exe を実行し、右下のトグルスイッチを押せばシミュレーションを開始できます。
感染率などのいくつかのパラメーターは実行時にリアルタイムに変更できます。

S と I が隣り合っているとき、一定の割合で感染が発生するようになっています (したがって、数式で表した状況とは厳密には異なります)。
また、I と R は一定の割合でそれぞれ R と S に移動します。

実行前に設定するパラメーター:

  • 高さ
  • SIR の初期人口比

実行中も設定できるパラメーター:

  • 感染率
  • 回復率
  • 免疫喪失率
  • Looping Map: マップの端でループするかどうか
  • ターンの時間間隔

Epidemic Simulator

 

実装方法については技術的に目新しいところはありませんが、特徴としては以下が挙げられます。

(1) シミュレーションの演算は UI スレッドではなく、バックグラウンド スレッドで実行
(2) 各フレームで画像データを生成して、Image コントロールで表示

(1) については、重い処理を UI スレッドで実行するとアプリケーションがフリーズしてしまうため、
各フレームで非同期的にデータのスナップショットを取得しています。
技術的な説明は、以前に

で書いた通りです。

バージョン情報
.NET Framework 4.5
ReactiveProperty 3.6.0
ToggleSwitch 1.1.2

参照
稲葉寿先生還暦記念祝賀研究集会「数理人口学・数理疫学・構造化個体群モデル」

広告
カテゴリー: ツール, 数学. タグ: . Leave a Comment »

3D における回転の表現と相互変換

以前に投稿した WPF で 3D オブジェクトを回転させるではオブジェクトの回転の状態を行列で表していましたが、
3 次元空間における回転を表現する方法は、次のように何通りか考えられます。
なお、原点を中心に回転させるものとし、回転角度は、回転軸を表すベクトルの方向に右ねじを押し込む場合を正とします。

  1. 行列
    • 3 次正方行列。とくに、回転を表すものは直交行列であり、P^{-1} = {}^t P が成り立つ
    • WPF では Matrix3D 構造体
  2. 四元数 (クォータニオン)
  3. オイラー角 (ロール、ピッチ、ヨー)
    • 任意の回転を 3 回の単純な回転の合成により表す。詳細の定義はロール・ピッチ・ヨー ― Kinectで学ぶ数学を参照
    • 座標系ごと回転させる
    • ロール、ピッチ、ヨーをそれぞれどの軸に対応させるか、どの順番で作用させるかで結果が変わってしまうため注意が必要
  4. 特定の 2 点の回転前後の座標
    • とくに、直交する 2 つのベクトルを用いるとよい

3D のプログラミングをしていると、API によってどの表現を利用するかが異なることがあります。
以下では、これらの表現を互いに変換する方法について考えます。

■ 行列 → 2 点の座標

任意のベクトルに対して行列を作用させれば回転後のベクトルが求められます。
WPF では演算子 * が用意されています。

■ 四元数 → 行列

四元数の成分から行列の各成分を計算できます (詳細は四元数による回転行列の表現を参照)。
WPF では Matrix3D.Rotate メソッドとして用意されています。

■ オイラー角 → 行列または四元数

ここでは、ヨーは y 軸、ピッチは x 軸、ロールは z 軸のまわりの回転を表し、この順に適用されるとします。
(元の座標系における) ヨー、ピッチ、ロールを表す回転行列をそれぞれ R_y, R_p, R_r とし、それぞれの回転角度を \theta_y, \theta_p, \theta_r とします。
すなわち、

R_y = \left( \begin{array}{ccc} \cos \theta_y & 0 & \sin \theta_y \\ 0 & 1 & 0 \\ -\sin \theta_y & 0 & \cos \theta_y \end{array} \right), R_p = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \theta_p & -\sin \theta_p \\ 0 & \sin \theta_p & \cos \theta_p \end{array} \right), R_r = \left( \begin{array}{ccc} \cos \theta_r & -\sin \theta_r & 0 \\ \sin \theta_r & \cos \theta_r & 0 \\ 0 & 0 & 1 \end{array} \right)

このとき、座標系ごとヨー、ピッチ、ロールの順に適用した回転は、元の座標系で R_y R_p R_r で表されます (適用の順序が逆になる)。
証明は次のようにできますが、実際の回転をイメージするとわかりやすいと思います。

証明

ベクトル {\bf x} にヨーを作用させると、
R_y {\bf x}
次に、これにピッチを作用させるには、いったん座標系を戻して R_p を作用させるから、
R_y R_p R_y^{-1} \cdot R_y {\bf x} = R_y R_p {\bf x}
同様に、これにロールを作用させると、
(R_y R_p) R_r (R_y R_p)^{-1} \cdot R_y R_p {\bf x} = R_y R_p R_r {\bf x}
(証明終)

WPF での実際の演算では四元数を使うとよいでしょう。

■ 2 点の座標 → オイラー角

ここでは、2 点として {\bf e}_z = (0, 0, 1), {\bf e}_y = (0, 1, 0) を選びます。
これらの回転後のベクトルがそれぞれ {\bf u}, {\bf t} で与えられたとすると、オイラー角は以下の手続きにより求められます。

まず、{\bf u} の x 要素と z 要素の比はピッチおよびロールの影響を受けないから、
\tan \theta_y = \dfrac{{\bf u}_x}{{\bf u}_z}
により、\theta_y, R_y が決まる。

ピッチおよびロールの決め方から、
\tan \theta_p = \dfrac{- (R_p {\bf e}_z)_y}{(R_p {\bf e}_z)_z}, \tan \theta_r = \dfrac{- (R_r {\bf e}_y)_x}{(R_r {\bf e}_y)_y}

また、前項と同様に考えて、{\bf u} = R_y R_p {\bf e}_z, {\bf t} = R_y R_p R_r {\bf e}_y であるから、
R_p {\bf e}_z = R_y^{-1} {\bf u}, R_r {\bf e}_y = R_p^{-1} R_y^{-1} {\bf t}

したがって、\theta_p, \theta_r も順に決まる。
なお、この決め方の場合、回転角度の範囲は、
- \pi < \theta_y \le \pi, - \dfrac{\pi}{2} \le \theta_p \le \dfrac{\pi}{2}, - \pi < \theta_r \le \pi
となる。
(終)

arctan を求めるには、Math.Atan2 メソッドを使うとよいでしょう。

以上により、行列 → 2 点の座標 → オイラー角 → 四元数 → 行列と変換する方法が与えられたので、
回転の表現を相互に変換できるようになります。

では、これらを実装してみます。
ベクトル、行列、四元数を扱うための System.Numerics.Vectors という、WPF のライブラリよりも高機能なライブラリもあるのですが、
今回は WPF のライブラリのみを利用して実装したいと思います。

全体のソースコードは RotationTest (GitHub) にあります。
これらのメソッドを呼び出すテストが付いています。

前回: WPF で 3D オブジェクトを回転させる
次回: Leap Motion で手の回転状態を取得する

バージョン情報
.NET Framework 4.5

参照

命題論理の複雑な問題を解く (C#)

// この投稿は 数学 Advent Calendar 2016 の 23 日目の記事です。

前回の命題論理を実装する (C#) では命題論理を扱うためのライブラリを作成しました。
今回は、このライブラリを利用して複雑な問題をいくつか解いてみます。

騎士と悪漢

まず、前回でも紹介した書籍 Raymond Smullyan「記号論理学: 一般化と記号化」から、問題を 2 つ引用します。


すべての住民は騎士または悪漢のいずれかであり、騎士はつねに本当のことを言い、悪漢はつねに噓をつく。

問題 1.3
住民 A, B がおり、A は「私たちは二人とも悪漢だ」と言った。
A, B はそれぞれ騎士か悪漢か?

問題 1.5
住民 A, B がおり、A は「私たちは同種 (二人とも騎士または二人とも悪漢のいずれか) だ」と言った。
A, B はそれぞれ騎士か悪漢か?


普通に考えて解くこともできますが、この状況を記号で表すことを考えます。
住民 X が騎士であるという命題を k_X とし、X が命題 p を主張したとすると、
k_Xp の真偽性は等しいということになり、k_X \equiv p が成り立ちます。
したがって、問題 1.3 では k_A \equiv (\lnot k_A \land \lnot k_B) が、問題 1.5 では k_A \equiv (k_A \equiv k_B) が成り立ちます。

前回作成したライブラリ Blaze を使って、次のように解くことができます。

実行すると、問題 1.3 では A が悪漢、B が騎士と確定します。
問題 1.5 では B が騎士と確定しますが、A は確定できません。

PropositionsConsole

 

数当てゲーム

さて次に、David Gale「Tracking the Automatic ANT: And Other Mathematical Explorations」という書籍から、
数当てゲームを紹介します。


  • プレイヤーは 2 人
  • それぞれのプレイヤーに数が割り当てられており、自身の数を知っているが、相手の数は知らない
  • 2 人の持つ数は、連続する 2 つの自然数 (つまり、差は 1)
  • ターン制で、相手の数を当てれば終了 (わかったら宣言し、わからなければ発声しない)
  • 先攻・後攻の区別はなく、いずれが宣言してもよい

どちらのプレイヤーが何ターン目に相手の数を当てることができるか?

人型ロボットと数当てゲームをする人のイラスト


自分の数が 9 だとすると、与えられた条件から相手の数は 8 または 10 の 2 通りに絞られます。
ここまではすぐにわかるのですが、いずれかを確定させるには相手の将来の行動あるいは実際の行動の原因を推測する必要があります。
相手がどう考えているかを考える、という日常生活においても重要な能力を問われているようでもあります。

小さい数値でいくつか試してみましょう。
2 人のプレイヤーには対称性がありますが、ここでは大小関係を固定して、プレイヤー A が n、プレイヤー B が n+1 を持つとします。
以下では、各プレイヤーの心理状況を書くことにします。

n=1 の場合、

  • A「B=2 しかありえない。」
  • B「A=1 または A=3。A=1 ならば、A は 1 ターン目に当てるはず。A=3 ならば、A は 1 ターン目に静観するはず。」

したがって、A が 1 ターン目で当てて終了します。
n=2 の場合、

  • A「B=1 または B=3。B=1 ならば、B は 1 ターン目に当てるはず。B=3 ならば、B は 1 ターン目に静観するはず。」
  • B「A=2 または A=4。A=2 ならば、A は 2 ターン目に当てるはず (B が 1 ターン目に静観するため)。
    A=4 ならば、A は 2 ターン目まで静観するはず。」

したがって、1 ターン目は両方静観し、A が 2 ターン目で当てて終了します。
このように、相手がどのように行動してくるかをあらかじめシミュレートしておけば、実際の行動に応じて相手の数を確定できることになります。

この方針でプログラミングしたものが NumberGuessConsole (GitHub) です (コードが長いため、ここには記載しません)。
このプログラムを A=8, B=9 で実行した結果です。A が 8 ターン目に当てました。
(「X @ n」は X が n ターン目に当てることを表します。)

NumberGuessConsole

各プレイヤーの脳内であらかじめシミュレーターを実行することで、値を確定させるための条件を得ます。
A の Knowledge には ((B = 7)≢(B = 9))∧((B = 7)⇒(B @ 7))∧((B = 9)⇒~(B @ 7)) と入っています。
B が 7 ターン目に静観したため、A は B=9 だとわかりました。

NumberGuessConsole

結果として、小さいほうの数 n を持つプレイヤーが n ターン目に相手の数 n+1 を当てることになります。
(プログラミングをしなくても、同様の方針により少し複雑な数学的帰納法で証明できます。)

前回: 命題論理を実装する (C#)

作成したライブラリ
Blaze (GitHub)
Blaze (NuGet Gallery)

作成したサンプル
PropositionsConsole (GitHub)
NumberGuessConsole (GitHub)

素材
不気味の谷現象のイラスト (いらすとや)

バージョン情報
.NET Framework 4.5

参照
Raymond Smullyan「記号論理学: 一般化と記号化」
David Gale「Tracking the Automatic ANT: And Other Mathematical Explorations」

カテゴリー: 数学. タグ: . 1 Comment »

命題論理を実装する (C#)

// この投稿は 数学 Advent Calendar 2016 の 16 日目の記事です。

命題といえば高校数学で「かつ」「または」「ならば」などを学習すると思いますが、
この投稿では、これらを一般化、形式化した命題論理をプログラミングで表現することを目指します。
プログラミングをしない方も、適宜飛ばせばなんとか読めると思います。

用語や定義については、主に Raymond Smullyan「記号論理学: 一般化と記号化」の第 7 章「命題論理入門」に記載されているものを使い、
解説と実装をしていきます。
プログラミング言語としては C# (.NET) を利用します。

Raymond Smullyan「記号論理学: 一般化と記号化」
「記号論理学: 一般化と記号化」
高橋 昌一郎, Raymond Smullyan, 川辺 治之

 

論理結合子

任意の命題 p は、真偽値 (truth value) を持ちます。真偽値とは、真 (T) または偽 (F) です。
そして、「かつ」「または」「ならば」などを一般化した概念を論理結合子 (logical connective) と呼び、次のものがあります。

  1. 否定 (~p)
    真偽値を反転させます。
  2. 論理積 (p ∧ q)
    p および q がともに真のときに限り真です。
  3. 論理和 (p ∨ q)
    p または q のうち、少なくとも一方が真のときに真です。
  4. 含意 (p ⇒ q)
    ~p ∨ q と同義です。
  5. 同値 (p ≡ q)
    (p ⇒ q) ∧ (q ⇒ p) と同義です。計算機においては、XOR 演算の否定と同義です。

論理式

さらに、論理式を次のように定義します。

  • 命題変数 (propositional variable): 命題を表す変数。値は、真または偽となりうる。
  • 命題定数 (propositional constant): 真を表す t および偽を表す f
  • 論理式 (formula): 命題変数、命題定数、論理結合子により構成される式。

ここまでの内容をもとに、論理式を実装します。 論理結合子は、否定のみが単項で、他はすべて二項です。

式を簡単に記述できるようにするため、静的メソッドを用意します。

 

恒真式および矛盾

論理式のうち、含まれる命題変数がどのような真偽値の組合せになっても真であるものを恒真式 (tautology) と呼びます。
逆に、含まれる命題変数がどのような真偽値の組合せになっても偽であるものを矛盾論理式または矛盾 (contradiction) と呼びます。
簡単な例を挙げると、p ∨ ~p は恒真式で、p ∧ ~p は矛盾です。

論理式が恒真式あるいは矛盾であるかどうかを判定するためのメソッドは、次のように実装できます。

以上のコードをまとめて、Blaze (GitHub) というライブラリとして公開しています。 NuGet でインストールできます。

 

ライブラリを使ってみる

さて、このライブラリを利用すると、命題論理に関する定理を証明できるようになります。
実際にいくつかやってみましょう。
Visual Studio で新規のコンソール アプリケーション プロジェクトを作成して、NuGet で Blaze をインストールします。
そして次のようなコードを記述します。

恒真式であるということは、その論理式がつねに成り立つということです。
このコードを実行することで、三段論法、背理法、対偶などを証明できます。

PropositionsConsole

 

次回は、このライブラリを利用してもう少し高度な問題を解いてみます。

次回: 命題論理の複雑な問題を解く (C#)

作成したライブラリ
Blaze (GitHub)
Blaze (NuGet Gallery)

作成したサンプル
PropositionsConsole (GitHub)

バージョン情報
.NET Framework 4.5

参照
Raymond Smullyan「記号論理学: 一般化と記号化」