2016-04-12

March 2016

Volume 31 Number 3

C# - 離散事象型シミュレーション: 人口増加の例

歴史的に見て、シミュレーションはさまざまな科学技術の発展に貢献してきました。人体をシミュレーションする医学モデルは、人体構造の研究に寄与します。「World of Warcraft」などの PC シミュレーションゲームはファンタジーの世界を生み出し、フライトシミュレーターはパイロットの地上訓練に役立ちます。テロリストによる攻撃や深刻な感染病の流行など、可能性のある危機への対策の調査にも、さまざまなシミュレーションプログラムが利用されています。映画「ジュラシックパーク」での恐竜のシミュレーションは、シミュレーションの幅広い応用範囲やその可能性を示唆しています。

シミュレーションとは、モデルを設計して現実のシステムやプロセスを模倣するテクニックです。このモデルにシステムの機能や動作をすべてカプセル化します。シミュレーションでは、時間をかけてモデル化したシステムを実行します。シミュレーションの設計には、以下の複数の工程があります。

モデル化するシステムの定義。対象となる問題を調べ、環境の特性を見極めて、達成すべき目標を定めます。
モデルの系統化。すべての変数と論理関係を定義して、必要なフローダイアグラムを作成します。
データの定義。目的の結果を得るためにモデルに必要なデータを定義します。
モデルのコンピューター実装の作成。
実装モデルの検証。実装したモデルが、設計の要件を満たすかどうかを確認します。
現実と対比した検証。シミュレーターがシミュレーションする現実のシステムを実際に表現できているかどうかを確認します。
テスト。シミュレーターを使用して目的のデータを生成できるかどうかテストします。
結果の分析と解釈。シミュレーターの結果を分析および解析して、得られた結果に基づいて決断を下します。
ドキュメント化。作成したモデルとツールに使用したシミュレーターを文書化します。

シミュレーションには、一般的に連続プロセス型と離散事象型があります。たとえば、気象システムのシミュレーションでは、すべての要素が絶えず変化するため、追跡が連続的に行われます。その結果、時間と温度変化の関係は連続曲線で表現されることになります。一方、飛行機の離着陸は特定の時点で行われるため、これをシミュレーションする場合、その時点またはその事象のみを考慮して、その他のものは無視できます。この後者のシミュレーションを離散事象型シミュレーション (DES) といいます。そして、これが今回のテーマです。

離散事象型シミュレーション

DES は、ある時点から次の時点へと時系列に順序付けた個々の事象のシーケンスとして、システムまたはプロセスをモデル化します。そのため、DES シミュレーションでは、通常、時間が実時間よりもはるかに短くなります。

DES の開発時には、以下の 6 つの主要要素を考慮することになります。

現実のシステムの要素を表すオブジェクト。オブジェクトにはプロパティがあります。オブジェクトは、事象に関連します。オブジェクトはリソースを消費します。そして、オブジェクトは時間の経過に従ってキューを出入りします。先ほどの飛行機の離着陸シナリオでオブジェクトに相当するのは、飛行機です。医療システムでオブジェクトになり得るのは、患者や臓器です。倉庫システムでは、在庫製品がオブジェクトになるでしょう。オブジェクトは、他のオブジェクトまたはシステムと相互作用すると想定され、シミュレーション中であればいつでも作成される可能性があります。

各オブジェクト固有の特徴を表すプロパティ。プロパティは、シミュレーションで行われるさまざまなシナリオへの応答を決定するために格納される、各オブジェクト固有の特徴 (サイズ、着陸時刻、性別、価格など) です。プロパティの値は変更できます。

システム (特にオブジェクト) で発生する可能性のある事象。たとえば、飛行機の着陸、倉庫への製品の到着、特定の病気の症状などです。事象が発生する順序は決まっておらず、同じ事象が再発生することもあります。

オブジェクトにサービスを提供するリソース。たとえば、空港の滑走路、倉庫の保管区画、病院に勤務する医者などです。リソースに空きがなく、オブジェクトがそのリソースを必要とする場合、そのオブジェクトをキューに配置して、リソースが利用可能になるまで待機しなければなりません。

オブジェクトを整理するキュー。現在空きのないリソースが解放されるまで待機するために、オブジェクトを整理するパイプです。キューは最大容量が決まっています。キューには、先入れ先出し法 (FIFO)、後入れ先出し法 (LIFO)、または何らかの条件や優先度 (病気の進行、燃料消費など) に基づく方法など、さまざまな呼び出し方式があります。

シミュレーションに不可欠な時間 (実生活と同様)。時間を計測するには、シミュレーション開始時にクロックを作動します。その後、特定の期間 (発着時刻、輸送時間、特定の症状が治癒するまでの時間など) の追跡にこのクロックを使用します。このように時間を追跡するのが基本で、その結果、事象の次の発生時刻を把握できるようになります。

シミュレーションのプログラミングは複雑になる可能性があるため、開発が簡単になるように、シミュレーションパラダイムのすべての要件を盛り込んだ言語を作成する試みが何回も行われています。このような言語の 1 つが、1960 年代に Ole Johan と Kristen Nygaard によって考案された SIMULA です。この言語は、現在を代表するプログラミングパラダイムとしてオブジェクト指向プログラミング (OOP) を初めて導入しています。最近では、シミュレーションの作成時に開発者が必要なものを組み込むことができる、パッケージ、フレームワーク、ライブラリの作成に注目が集まっています。このようなライブラリは、C#、C++、Java、Python のような一般的な言語から呼び出せるようになっています。

Nance は自著「A History of Discrete Event Simulation Programming Languages」(離散事象型シミュレーションプログラミング言語の歴史、英語) の中で、DES プログラミング言語が満たすべき 6 つの最低要件を提案しています (http://eprints.cs.vt.edu/archive/00000363/、英語)。

乱数の生成
統一ランダム変数以外の変数を許可するプロセス変換器
オブジェクトの作成、操作、削除を容易にするリスト処理
モデルの動作の要約を提供する統計分析
膨大なデータの表示と意思決定を支援するレポートの作成
時系列メカニズム

こうした要件は C# ですべて満たすことができるため、今回は、C# で開発した人口増加の離散事象型シミュレーションを紹介します。

人口増加の離散事象型シミュレーション

個体群の増加は、時間や空間の変化や環境との相互作用に応じて、個体群 (動植物) の数が変化する傾向を研究する際に考慮される多くの側面の 1 つです。個体群とは同時期に生息する同種の生命体のグループで、同じ空間に存在することもあれば、特定の特性によって区分されることもあります。

個体群増加の研究が重要なのは、個体群増減のしくみを正しく理解することで、科学者は、生物多様性の保全、資源の利用、気候変動、汚染、医療、輸送の必要性などの変化について、より正確に予測できるようになります。また、個体群の増加または減少について考慮する際に重要な側面となる、生命体どうしの相互作用や環境との作用についての洞察を得ることができます。

今回は、人口増加の DES を取り上げます。目標は、時間の経過と人口増加の関係を観察し、最終結果 (人口の規模、高齢者、若年者など) を分析することによって統計を得ることです。人口は、m 人の男性と n 人の女性から始まり、各個人には年齢が設定されています。当然、m 人と n 人は 0 以上にし、見当違いなシミュレーションにならないようにしなければなりません。このモデルのオブジェクトは個人です。

個人は、パラメーター λ = 18 年のポアソン分布で特定の年齢に達すると、別の個人と結婚を意識するようになります (現時点では、ポアソン関数の正規分布や指数分布について完全に理解している必要はありません。詳しくは後ほど説明します)。

独身で、異性との年齢差が 5 歳以下の場合にのみ結婚でき、その確率は 50% 未満です。図 1 は、2 人の関係が終了する確率の例を示しています。

図 1. 関係が終了する確率

平均年齢	確率
14-20	0.7
21-28	0.5
29+	0.2

結婚した 2 人は、パラメーター λ = 8 年の指数分布で特定の年数が経過すると、子供を授かる可能性があります。

女性は、パラメーター µ = 28 年と σ² = 8 年の正規 (鐘状) 分布に従った年齢の場合に妊娠できます。すべての女性は、パラメーター µ = 2 年と σ² = 6 年の正規分布の人数の子供を授かります (パラメーター µ は平均年齢を、パラメーター σ² は年齢の変動性の尺度を表します)。

すべての人間は、パラメーター λ = 70 年のポアソン分布の平均寿命を持ちます (λ は、平均寿命の期待値を表します)。

ここまでの説明から、以下に示す複数の事象が特定されます。

関係の開始
関係の終了
妊娠
出産
死亡

各事象には、事象が発生する時期を決定する離散型ランダム変数が付属します。

確率分布

離散型ランダム変数とは、有限または可算無限の一連の値を持つ変数です。つまり、値は有限または無限の値の並び (1, 2, 3 ....) としてリストすることができます。離散型ランダム変数の場合、その確率分布は、各値に確率を割り当てるグラフ、表、または式になります。確率の総和は 1 にならなければなりません。個々の確率は 0 ～ 1 の値でなければなりません。たとえば、等面の (すべての面が同じ確率で出現する) サイコロを振った場合、起こり得る結果を表す離散型ランダム変数 X は、確率分布 X(1) = 1/6, X(2) = 1/6, ...., X(6) = 1/6 を保持します。すべての面が等しい確率で出現するため、それぞれの面に割り当てられる確率は 1/6 です。

パラメーター l と µ は、対応する分布の平均 (期待値) を示します。平均は、ランダム変数が平均して取る値を表します。つまり、総和 E = [(起こり得る各結果) × (その結果が得られる確率)] で、この場合の E が平均を示します。サイコロの場合の平均は、E = 1/6 + 2/6 + 3/6 + 4/6 + 5/6 + 6/6 = 3.5 になります。結果 3.5 は、実のところ、サイコロが取り得るすべての値の中間で、サイコロを何回も振った場合の期待値になります。

パラメーター σ2 は分布の分散を示します。分散はランダム変数として可能性のある値の離散状態を表し、常に非負の値になります。0 に近い小さな分散は、値どうしが接近していて、平均付近に集まっていることを示し、1 に近い大きな分散は、値どうしが大きく離れ、平均から遠く離散していることを示します。

ポアソン分布は離散型の分布で、単位時間あたりの事象数に関する確率を表します (図 2 参照)。ポアソン分布は、通常、事象が発生する確率が小さく、事象が発生する機会の数が多い場合に適用されます。書籍のミスプリントの数、ビジネスセンターに到着する顧客数、信号で止まる車の台数、特定の年齢層の 1 年あたりの死者数はすべて、ポアソン分布を適用できる例です。

図 2. パラメーター λ = 18 のポアソン分布

指数分布は、ポアソン過程における事象間の時間を示します (図 3 参照)。たとえば、特定の期間中にビジネスセンターに到着した顧客数を示すポアソン過程に取り組んでいる場合に、最初の顧客が到着するまでの時間を示すランダム変数に注目するとします。この目的には指数関数が有効です。指数関数は物理学の過程に適用することもでき、たとえば、λ で粒子の経過年数の比率を示して、粒子の寿命を表すことができます。

図 3. パラメーター λ = 8 の指数分布

正規分布は、左右に偏らずに、中央値付近に集まる傾向がある確率を示します (図 4 参照)。正規分布は対称で、平均値付近にピークを 1 つ持つ鐘状の濃度曲線になります。分布の 50% は平均の左に位置し、残りの 50% が右に位置します。標準偏差は鐘状曲線の広がりまたは胴回りを示し、標準偏差が小さくなるほどデータの集中度が高まります。平均と標準偏差はどちらも、正規分布のパラメーターとして指定する必要があります。血圧、身長、測定誤差など、多くの自然現象は正規分布に忠実に従います。

図 4. パラメーター µ = 28 と σ2 = 8 年の正規分布

ここからは、今説明した DES を、C# など、人気のある、洗練され優れた言語で実装する方法を示します。

実装

このシミュレーションを開発するため、OOP パラダイムのメリットをすべて活用します。これは、できる限り読みやすいコードにするという考えが根底にあります。科学アプリケーションは複雑になり、理解するのが難しくなりがちなので、他の開発者がコードを理解できるよう、できる限り簡潔になるように努めます。今回は、C# のコンソールアプリケーションとして、図 5 に示す構造で開発します。

図 5. シミュレーションアプリケーションの構造

パスを論理的にすることは重要です。つまり、名前空間の構造が常識的に理解できるようにします (Simulation.DES.PopulationGrowth)。

Events フォルダーに含まれる Event.cs ファイルは、以下のように、シミュレーションで発生する可能性のあるすべての事象を表す列挙型を定義します。

namespace DES.PopulationGrowth.Events
{
public enum Event
  {
Capable Engaging,
Birth EngageDisengage,
GetPregnant,
ChildrenCount,
TimeChildren,
    Die
  }
}

Objects フォルダーは、個人に関連するすべてのクラスを保持します。これらのクラスのコードは最も OOP のメリットを利用します。ここには、Individual クラス、Male クラス、および Female クラスという個人に関連する 3 つのクラスがあります。最初のクラスが抽象クラスで、他の 2 つのクラスはこの抽象クラスから継承します (男性も女性も個人です)。コーディングのほとんどは、Individual クラスで行います。図 6 に、Individual クラスのプロパティとメソッドを示します。

図 6. Individual クラスのプロパティとメソッド

各プロパティの説明を以下に示します。

Age: 個人の年齢。
Couple: 個人が独身の場合は Null、結婚している場合は相手の個人を示します。
Engaged: 結婚している場合は true、していなければ false です。
LifeTime: 個人が死亡する年齢です。
RelationAge: 結婚できる年齢。
TimeChildren: シミュレーションにおける個人が子供を授かるまでの時間 (年齢ではありません)。

各メソッドは以下のとおりです。

Disengage: 2 人の個人の関係を終了します。
EndRelation: 図 1 の確率表に従って、関係が終了するかどうかを判断します。
FindPartner: 特定の個人のパートナー候補を検索します。
SuitablePartner: 特定の個人が関係構築の開始に適しているかどうか (年齢差や異性かなど) を判断します。
SuitableRelation: ある個人が関係構築を開始できるかどうかを、独身かどうか、関係構築を開始できる年齢かどうかに基づいて判断します。
ToString: 個人情報を文字列として表現するためにオーバーライドします。

クラスは、すべてのプロパティとコンストラクターを定義することから始めます。コンストラクターは、その個人の年齢を設定するだけです。

public abstract class Individual
  {
public int Age { get; set; }
public int RelationAge{ get; set; }
public int LifeTime{ get; set; }
public double TimeChildren{ get; set; }
public Individual Couple { get; set; }
protected Individual(int age)
{
  Age = age;
}

SuitableRelation メソッドと SuitablePartner メソッド、および Engaged プロパティは、ごく簡単な論理テストを実行するだけです。

public bool SuitableRelation()
{
return Age >= RelationAge&& Couple == null;
}
public bool SuitablePartner(Individual individual)
{
return ((individual is Male && this is Female) ||
(individual is Female && this is Male)) &&Math.Abs(individual.Age - Age) <= 5;
}  
public bool Engaged
{
get { return Couple != null; }
}

Disengage メソッドは、カップルの Couple フィールドに null を設定した後、個人の Couple フィールドに null を設定して関係を終了させます。また関係が終了するため、子供を授かるまでの時間も 0 に設定します。

public void Disengage()
{
Couple.Couple = null;
Couple = null;
TimeChildren = 0;
}

EndRelation メソッドは、基本的には確率表を解釈して、関係が終了する可能性を判断します。そのため、統一ランダム変数を使用します。この変数は範囲 [0, 1] の乱数値を作成します。これは関係が終了する割合に相当します。distributions 辞書は Simulation クラス (後ほど説明) で作成します。この辞書は、(事象と確率分布の) ペアを保持して、各事象と分布を関連付けます。

public bool EndRelation(Dictionary<Event, IDistribution> distributions)
{
var sample =
  ((ContinuousUniform) distributions[Event.BirthEngageDisengage]).Sample();
if (Age >= 14 && Age <= 20 && sample <= 0.7)
return true;
if (Age >= 21 && Age <= 28 && sample <= 0.5)
return true;
if (Age >= 29 && sample <= 0.2)
return true;
return false;
}

図 7 の FindPartner メソッドは、個人のインスタンスのパートナー候補を探します。このメソッドは入力として、人口リスト、シミュレーションの現在時間、および distributions 辞書を受け取ります。

図 7. FindPartner メソッド

public void FindPartner(IEnumerable<Individual> population, int currentTime,
  Dictionary<Event, IDistribution> distributions)
{
foreach (var candidate in population)
if (SuitablePartner(candidate) &&
candidate.SuitableRelation() &&
((ContinuousUniform) distributions[Event.BirthEngageDisengage]).Sample() <= 0.5)
{
// Relate them
candidate.Couple = this;
Couple = candidate;
// Set time for having child
var childTime = ((Exponential)  distributions[Event.TimeChildren]).Sample()*100;
// They can have children on the simulated year: 'currentTime + childTime'.
candidate.TimeChildren = currentTime + childTime;
TimeChildren = currentTime + childTime;
break;
}
}

最後に、ToString メソッドは、個人の文字列表現を定義します。

public override string ToString()
{
Return string.Format("Age: {0} Lifetime {1}", Age, LifeTime);
}

Male クラスは非常にシンプルです。

public class Male: Individual
{
public Male(int age) : base(age)
{
}
public override string ToString()
{
return base.ToString() + " Male";
}
}

Female クラスでは、妊娠、出産なども扱うため、もう少し複雑です。Female クラスの定義は、すべてのプロパティとコンストラクターの宣言から始まります。

public class Female :Individual
{
public bool IsPregnant{ get; set; }
public double PregnancyAge{ get; set; }
public double ChildrenCount{ get; set; }
public Female(int age) : base(age)
{
}
}

Female クラスのプロパティを以下に示します。

IsPregnant: この女性が妊娠しているかどうかを判断します。
PregnancyAge: 女性が妊娠できる年齢を判断します。
ChildrenCount: 生まれてくる子供の数を示します。

このクラスに含まれるメソッドは以下のとおりです。

SuitablePregnancy: この女性が妊娠できるかどうかを判断します。
GiveBirth: 女性の出産を示し、人口に新しい個人が加わります。
ToString: オーバーライド: 女性の情報を文字列として表現するために使用します。

SuitablePregnancy は、女性のインスタンスが妊娠する条件をすべて満たしている場合に true を出力する、テストメソッドです。

public bool SuitablePregnancy(intcurrentTime)
  {
return Age >= PregnancyAge && currentTime <= TimeChildren && ChildrenCount > 0;
}

図 8 に示す GiveBirth メソッドは、新しい個人を人口に加えて初期化するコードです。

図 8. GiveBirth メソッド

public Individual GiveBirth(Dictionary<Event, IDistribution> distributions, int currentTime)
  {
var sample =
  ((ContinuousUniform) distributions[Event.BirthEngageDisengage]).Sample();
var child = sample > 0.5 ? (Individual) newMale(0) : newFemale(0);
// One less child to give birth to
ChildrenCount--;
child.LifeTime = ((Poisson)distributions[Event.Die]).Sample();
child.RelationAge = ((Poisson)distributions[Event.CapableEngaging]).Sample();
if (child is Female)
{
(child as Female).PregnancyAge =
  ((Normal)distributions[Event.GetPregnant]).Sample();
(child as Female).ChildrenCount =
  ((Normal)distributions[Event.ChildrenCount]).Sample();
}
if (Engaged && ChildrenCount> 0)
{
TimeChildren =
  currentTime + ((Exponential)distributions[Event.TimeChildren]).Sample();
Couple.TimeChildren = TimeChildren;
}
else
TimeChildren = 0;
IsPregnant = false;
return child;
  }

統一サンプルを生成し、最初に新しい個人がそれぞれ 50% (0.5) の確率で女性または男性のどちらになるかを決定します。ChildrenCount 値を 1 デクリメントし、この女性が出産する可能性のある子供の数を 1 人減らします。残りのコードでは、個人の初期化と一部変数のリセットを行っています。

ToString メソッドは、女性を文字列として表現する方法を変えます。

public override string ToString()
  {
return base.ToString() + " Female";
}

個人に関連するすべての関数を個々のクラスに配置したため、Simulation クラスは大幅に簡潔になります。コードの冒頭で、プロパティとコンストラクターを宣言します。

public class Simulation
  {
public List<Individual> Population { get; set; }
public int Time { get; set; }
private int _currentTime;
private readonly Dictionary<Event, IDistribution> _distributions ;

プロパティと変数は名前が示すとおりで、その一部は既に説明しました。Time は、シミュレーションの開始からの継続時間 (年数) を表し、_currentTime はシミュレーションでの現在年を表します。今回の場合、コンストラクター (図 9 参照) に各個人のランダム変数の初期化が含まれているため、やや複雑です。

図 9. Simulation クラスのコンストラクター

public Simulation(IEnumerable<Individual> population, int time)
{
Population = new List<Individual>(population);
Time = time;
_distributions = new Dictionary<Event, IDistribution>
{
{ Event.CapableEngaging, new Poisson(18) },
{ Event.BirthEngageDisengage, new ContinuousUniform() },
{ Event.GetPregnant, new Normal(28, 8) },
{ Event.ChildrenCount, new Normal(2, 6) },
{ Event.TimeChildren, new Exponential(8) },
{ Event.Die, new Poisson(70) },
                                 };
foreach (var individual in Population)
{
// LifeTime
individual.LifeTime = ((Poisson) _distributions[Event.Die]).Sample();
// Ready to start having relations
individual.RelationAge = ((Poisson)_distributions[Event.CapableEngaging]).Sample();
// Pregnancy Age (only women)
if (individual is Female)
                {
(individual as Female).PregnancyAge = ((Normal) _distributions[Event.GetPregnant]).Sample();
(individual as Female).ChildrenCount = ((Normal)_distributions[Event.ChildrenCount]).Sample();
}
}
}

最後に、Execute メソッド (図 10 参照) で、すべてのシミュレーションロジックを実行します。

図 10 Execute メソッド

public void Execute()
{
while (_currentTime< Time)
{
// Check what happens to every individual this year
for (vari = 0; i<Population.Count; i++)
{
var individual = Population[i];
// Event -> Birth
if (individual is Female&& (individual as Female).IsPregnant)
Population.Add((individual as Female).GiveBirth(_distributions,
  _currentTime));
// Event -> Check whether someone starts a relationship this year
if (individual.SuitableRelation())
individual.FindPartner(Population, _currentTime, _distributions);
// Events where having an engaged individual represents a prerequisite
if (individual.Engaged)
{
// Event -> Check whether arelationship ends this year
if (individual.EndRelation(_distributions))
individual.Disengage();
// Event -> Check whether a couple can have a child now
if (individual is Female &&
(individual as Female).SuitablePregnancy(_currentTime))
(individual as Female).IsPregnant = true;
}
// Event -> Check whether someone dies this year
if (individual.Age.Equals(individual.LifeTime))
{
// Case: Individual in relationship (break relation)
if (individual.Engaged)
individual.Disengage();
Population.RemoveAt(i);
}
individual.Age++;
_currentTime++;
}
}

外側のループの繰り返しは、シミュレーションの経過年数を表します。内側のループは、特定の年に発生する可能性がある事象に対処します。

時間経過に従って人口が増加するようすを見るために、新しいコンソールアプリケーションをセットアップします (図 11 参照)。

図 11. 時間経過に応じた人口の表示

static void main()
{
var population = new List<Individual>
{
new Male(2),
new Female(2),
new Male(3),
new Female(4),
new Male(5),
new Female(3)
};
var sim = new Simulation(population, 1000);
sim.Execute();
// Print population after simulation
foreach (var individual in sim.Population)
Console.WriteLine("Individual {0}", individual);
Console.ReadLine();
}

図 12 に、1,000 年後の人口を示します。

図 12. 1,000 年後の人口

今回は、離散事象型シミュレーションを開発して、時間経過に応じて人口が増加するようすを調べました。このオブジェクト指向のアプローチは、読者が必要に応じて試行や改良を行えるように、読みやすく、わかりやすいコードにすることを心掛けました。

Arnaldo Pérez Castaño は、キューバを拠点とするコンピューター科学者です。彼は、『JavaScript Fácil』(Marcombo S.A、2013 年)、『HTML y CSS Fácil』(Marcombo S.A、2014 年)、『Python Fácil』(Marcombo S.A、2015 年) というプログラミング書籍シリーズを執筆しています。Visual Basic、C#、.NET Framework、および人工知能を専門とし、nubelo.com (英語) でフリーランサーとして独自のサービスを提供しています。映画や音楽にも情熱を注いでいます。連絡先は arnaldo.skywalker@gmail.com (英語のみ) です。

この記事のレビューに協力してくれたマイクロソフト技術スタッフの James McCaffrey に心より感謝いたします。

次の方法で共有

C# - 離散事象型シミュレーション: 人口増加の例

離散事象型シミュレーション

人口増加の離散事象型シミュレーション

確率分布

実装

その他のリソース