※本記事は、Stanford大学のAA228V講座「Validation of Safety Critical Systems」におけるImportance Samplingに関する講義内容を基に作成されています。講義の詳細情報やコース内容は公式ウェブサイト（https://aa228v.stanford.edu/ ）でご覧いただけます。また、講義で使用されているテキストブックは https://algorithmsbook.com/validation/ で入手可能です。
本記事では講義内容を要約しております。なお、本記事の内容は原著作者の見解を正確に反映するよう努めていますが、要約や解釈による誤りがある可能性もありますので、正確な情報については原講義をご参照ください。
登壇者は、Stanford大学のポスドク研究員であるSydney Katz氏（詳細は https://sydneymkatz.com/ ）と、航空宇宙工学准教授および（兼任）コンピュータサイエンス准教授のMykel Kochenderfer氏（詳細は https://mykel.kochenderfer.com/ ）です。このコースについて詳しく知りたい方や受講を希望される方は、Stanford Onlineのウェブサイト（https://online.stanford.edu/courses/a... ）をご参照ください。
1. 前回の復習: 失敗分布からのサンプリング
1.1 Markov Chain Monte Carlo (MCMC)の課題
Sydney Katz：前回の講義では、失敗分布からサンプルを抽出しようとしていました。私たちは失敗分布の正規化された密度は計算できませんが、分子（unnormalized density）は計算できることを説明しました。そこで、この非正規化密度からサンプルを抽出するための異なる方法について議論していました。具体的には棄却サンプリングとMarkov Chain Monte Carlo（MCMC）を取り上げ、MCMCについての議論の途中で終わりました。
MCMCの主な課題の一つは、理論的には無限のサンプル数の極限では失敗分布からのサンプルを得られますが、実際には常に有限のサンプル集合を扱わなければならないという点です。実用的な観点から、この問題にどう対処するかについて考える必要があります。
私たちはいくつかの対処方法について話しました。一つはバーンイン期間を設けることで、最初の数サンプルを捨てるというものです。また、サンプル間の相関を減らすために「間引き（thinning）」を行い、例えば5番目や10番目のサンプルだけを保持するという方法もあります。
最後に議論していたのは「スムージング（smoothing）」というアイデアです。これが今日の最初のパートで深く掘り下げるトピックです。私たちは前回、バイモーダル（二つのモードを持つ）失敗分布の例を見ていました。そこでナイーブなMCMCを適用したとき、開始したモードからのサンプルしか得られないという問題が生じました。これは次のサンプルを提案するために使用していた分布（カーネルと呼ばれる）が、もう一方の失敗モードに対してどうしても非常に低い確率しか割り当てられないためでした。
理論的には無限のサンプルを取れば、最終的にはモード間のジャンプが起こりますが、有限のサンプル数では、単一のモードからのサンプルしか得られない可能性が高いです。この問題を解決するために、失敗分布をスムージングして、モード間の移動を助けるアプローチを考えました。
1.2 スムージング手法の詳細説明
Sydney Katz：スムージング手法の詳細に入りましょう。まず、delta(τ)と呼ぶ距離関数を定義します。この関数は軌道が失敗の場合は0を返し、それ以外は正の値を返します。これは軌跡の「失敗までの距離」を表しています。一つの簡単な方法としては、ロバストネス値の最大値と0を取る方法があります。ロバストネス値が失敗でない場合は正になるので、これによって失敗でない場合は正の値、失敗の場合は0という条件を満たせます。
この距離関数を使って、失敗分布の非正規化密度を書き換えることができます。以前は「τがΨに含まれない」という条件を使っていましたが、これを「delta(τ)が0に等しい」という同等の条件に置き換えることができます。
ここで重要なのは、この指示関数（indicator function）が非常に小さな分散を持つ正規分布と類似していると考えられることです。指示関数は、delta(τ)が0のときに1、それ以外で0となりますが、これを0を中心とした小さな分散を持つ正規分布で近似できます。この近似を使って、非正規化密度の指示関数を滑らかな近似に置き換えることができます。
この操作の利点は、失敗ではない軌道に対しても非ゼロの確率を割り当てることになります。特に失敗に近い成功には高い確率、失敗から遠い成功には低い確率が割り当てられるようになります。興味深いのは、イプシロン（ε）というパラメータを調整して、分布のスムージングの程度をコントロールできる点です。εが0のとき、スムージングはなく、元の指示関数に戻ります。εが無限大に近づくと、正規分布は無限に広がり実質的に一様分布となるため、非正規化密度は公称軌道分布に近づきます。
実際に使う際には、スムージングによって失敗でないサンプルも含まれるようになりますが、これは問題ではありません。最終的なサンプル集合では、失敗でないサンプルを棄却すれば良いのです。これは前回説明した棄却サンプリングと同様のアプローチです。スムーズ化された分布を提案分布として使い、失敗でないものを棄却するのです。この方法は非直感的に思えるかもしれませんが、数学的に正当化できます。
1.3 多様体間ジャンプの問題と解決方法
Sydney Katz：ここで実際の例を使ってスムージングの効果を見てみましょう。前回の講義で扱ったバイモーダル失敗分布の例に戻ります。上のグラフの赤い線は、サンプリングしようとしている真の失敗分布を示しています。紫の線はイプシロン値が0.05の少しだけスムージングされた分布です。
このスムージングによって、分布の端が滑らかになり、以前はゼロだった領域に非ゼロの確率が割り当てられるようになりました。しかし、このわずかなスムージングでも、サンプル列を見ると、まだもう一方の失敗モードへジャンプできていないことがわかります。
次に棄却サンプリングを適用します。これはすべてのサンプルから失敗でないものを除外する処理です。そうすると失敗分布からのサンプルが残りますが、やはり片方のモードからのサンプルしかありません。問題は、もう片方のモードへのジャンプがまだ大きすぎることにあります。
そこでスムージングのパラメータをもう少し増やしてみましょう。これにより、非常に滑らかな分布が得られます。ここで見ると失敗モード間を行き来するジャンプが発生し、両方の側からのサンプルが得られるようになりました。唯一の欠点は、最終的に棄却することになる非失敗サンプルも取得することになるため、一部のサンプルが無駄になることです。しかし、これは可能な失敗モードの完全な空間を効率的に探索するために必要な犠牲と言えます。
中間的なケースも見てみましょう。ここではモード間を行き来する頻度が少し足りないことがわかります。限られたサンプルサイズでは、もう一方の失敗モードからいくつかサンプルを得られますが、完全に均等にはなりません。結局、このパラメータはシステムごとに調整する必要があるものです。
質問がありました。「カーネル分布を広げるだけでも同様の効果が得られるのでは？」確かにそれも可能です。しかし一般的には、特に数百次元の高次元空間では、どのようなカーネルがモード間のジャンプを可能にするかを設計するのは明らかではありません。また、カーネルを大きくしすぎると、他の面でパフォーマンスが低下する可能性もあります。
もう一つの質問は「良いイプシロン値を見つけたかどうかをどう判断するか」というものでした。研究論文では、ナイーブな棄却サンプリングを非常に長時間（例えば30日間）実行して基準を得て、その後、改良手法が数時間でより良い結果を出すことを示します。実世界では、システムに関するドメイン知識に基づいて、何かが欠けているように見えるかどうかを判断する必要があります。
1.4 スムージングパラメータ(ε)の調整方法
Sydney Katz：スムージングパラメータ(ε)の調整は、失敗分布のサンプリングにおいて重要な役割を果たします。先ほど示した例で見たように、εの値によって結果は大きく変わります。
εが0の場合、スムージングは全く行われません。これは元の指示関数に相当し、失敗でない軌道には確率が割り当てられません。この場合、MCMCは最初に入ったモードから抜け出すことができず、バイモーダル分布の両方のモードからサンプルを得ることができません。
一方、εが無限大に近づくと、正規分布は無限に広がり、実質的に一様分布になります。この極端な場合、私たちの非正規化密度は公称軌道分布に近づきます。つまり、スムージングのない元の分布とスムージングが無限大の分布の間で、適切なバランスを見つける必要があります。
実践的には、εの値は以下のトレードオフを考慮して調整する必要があります：
- εが小さすぎると：モード間のジャンプが起こりにくく、単一モードからのサンプルしか得られない
- εが大きすぎると：非失敗サンプルが多くなりすぎ、効率が低下する
最適なεの値はシステムによって異なり、特定のシステムに対して経験的に調整する必要があります。例えば、先ほど見た例では、完全に両方のモードからサンプルを得るためには大きめのεが必要でしたが、これは非失敗サンプルも多く取得することを意味します。
良いεを見つけたかどうかを判断する方法については、真の失敗分布が分からない場合は難しい問題です。研究では、ベースラインとして非常に長時間の実行結果と比較しますが、実際の応用では、システムに関するドメイン知識に基づいて判断することになります。
重要なのは、εを調整することで、失敗モード間の移動のしやすさと、無駄なサンプル（非失敗サンプル）の数のバランスを取ることです。この調整は、システムの特性やサンプリングの目的に応じて行う必要があります。
多くの場合、いくつかの異なるε値でテストを行い、サンプルの多様性と効率のバランスが最も良い値を選択するというアプローチが効果的です。
1.5 高次元問題へのスケーリング
Sydney Katz：ここまで紹介した手法は1次元のガウス分布だけでなく、より高次元の問題にもスケールできます。例えば、振り子システムは約50次元の問題ですが、これに対しても失敗分布からサンプルを抽出するようにスケールアップできます。
教科書の例題を見ていただければ、コードがどのように動作するかの詳細がわかりますが、実際に動作することを示すために簡単に紹介します。これは振り子の失敗分布からサンプリングしようとしている例です。ここでもスムージングを適用しています。スムーズな分布では失敗でないサンプルもいくつか得られていますが、棄却処理を行うと失敗分布からのサンプルが残ります。
ここで重要なのは、これらのサンプルはただの1次元の点ではなく、完全な軌道であるということです。つまり、初期状態とその軌道のすべての摂動をサンプリングしています。新しいサンプルを提案するたびに、完全に新しい軌道を提案して、それを受け入れるか棄却するかを決定します。
また、振り子が右に倒れるか左に倒れるかの確率は等しいはずですが、サンプルの分布が少しアンバランスに見えます。これは、先ほどの1次元の例と同様、サンプルがモード間を十分に行き来していないためです。スムージングパラメータをもう少し増やすことで、これを改善できます。実行には少し時間がかかりますが、スムージングを強くすることで、非失敗軌道が増え、振り子の失敗分布のより良い表現が得られます。このコードはノートブックに含まれているので、興味があればご確認ください。
高次元問題へのスケーリングを助ける他の方法もあります。これについては、おそらくMarkov Chain Monte Carloについての別の授業で詳しく学ぶことになるでしょう。特に役立つのは、尤度の勾配を使用してカーネルを定義することです。我々は単純なガウスカーネルを定義していましたが、これは基本的に空間をランダムウォークするものです。
しかし一般的に、高尤度領域を低尤度領域よりも頻繁にサンプリングしたいと考えます。そこで、勾配を利用して尤度が増加する方向にランダムウォークを導く方法があります。時々尤度が減少する方向にも進みますが、勾配を使ってそのウォークを導くことができます。詳しくは教科書の3.6.3.3節で説明しています。
また、Hamiltonian Monte Carlo（HMC）という非常に有名な手法や、現在多くの人が使用している最先端の手法であるNUTS（No-U-Turn Sampler）もあります。これらはすべてMCMCをより効率的にするために使用できます。
スケールアップを助けるもう一つのパラダイムは確率的プログラミングです。これは前回のパラメータ推定で少し触れましたが、基本的に確率モデルをコンピュータプログラムとして指定できるようにするものです。そして、確率プログラミング言語がそれを理解し、自動的に推論を実行します。これは高度なトピックであり、クイズには出ませんが、この分野を深く掘り下げたい場合は非常に興味深いものです。
確率プログラミングがいかに優れているかを少し紹介したいと思います。教科書には、Turing.jlというJulia言語のライブラリを使ったアルゴリズムがあります。基本的には、モデル関数を作成するだけでよく、システムのロールアウト方法をモデル化するだけです。初期状態が初期状態分布からサンプリングされ、各摂動が摂動分布からサンプリングされると指定します。そして、サンプリングしたい分布の密度が失敗分布のためのスムーズ化された密度であることをTuringに伝えるだけです。
すると、Turingが自動的にサンプリングを行ってくれます。このような構文をすべて理解する必要はありませんが、これらの要素を簡単に指定するだけで、Turingが推論を自動的に行ってくれることがわかります。Hamiltonian Monte CarloやNUTSなどの実装もTuringに含まれており、このコードブロックを実行するだけで失敗分布からサンプルを抽出できます。
Pythonを使いたい場合は、PyMC、Pyro、Stanなどの同様のライブラリがあります。Stanが最も古くからあるものですが、C++で書かれており、Pythonとの統合は少し難しい場合があります。
2. 失敗確率の推定
2.1 失敗確率の数学的定義
Sydney Katz：ここまでの内容をまとめると、我々は失敗分析を偽陽性（falsification）から始めました。プロジェクト1では、とにかく失敗を見つけようとしました。次に、失敗の完全な分布からサンプリングする方法について説明しました。これはやや難しい作業でした。そして今からは、さらに一歩進んで、失敗確率の推定を試みます。
失敗確率とは何か、数学的に考えてみましょう。失敗確率は次のように定義されます：公称軌道分布からの軌道に対して、それが失敗かどうかの期待値です。第3章から覚えているかもしれませんが、1か0の値をとる指標の期待値を取ると確率が得られます。
さらに、この期待値を積分として書き下すことができます。失敗確率は、すべての可能な軌道に対して、それが失敗かどうかを公称分布下での尤度で重み付けした積分となります。
興味深いことに、この式は前回の講義で説明した失敗分布の分母と正確に一致します。つまり、失敗確率は失敗分布の正規化定数なのです。ここに明確な関連性があります。失敗確率は、非正規化された失敗分布密度の下の面積に等しいことになります。
しかし、前回話したように、この正規化定数は一般的に厳密に計算することが非常に難しいものです。これが、失敗分布からサンプリングするための手法全体を開発する必要があった理由です。この講義では、その分母、つまり失敗確率を推定する方法について説明します。
今日と木曜日の講義で、いくつかの方法について説明します。まず直接推定法から始めます。この方法は実際にこのクラスですでに何度か見てきました。次に、直接推定法の改良版である重要度サンプリング（importance sampling）に進みます。これはこのクラスの基本的な概念の一つで、非常に有用で興味深いものです。皆さんがこの講義から得るものがあるとすれば、これがその一つになることを願っています。その後、適応的重要度サンプリングなど、重要度サンプリングの拡張について説明します。
前回、ファジング（fuzzing）がプロジェクト1の出発点になりうると言いましたが、重要度サンプリングはプロジェクト2の出発点となるでしょう。プロジェクト2ではプロジェクト1で使用したシステムの失敗確率を推定することになります。そして直接推定法のベースラインを上回る必要がありますが、その一つの方法が重要度サンプリングになります。
2.2 失敗分布との関係性（正規化定数）
Sydney Katz：失敗確率と失敗分布の関係性についてさらに詳しく見ていきましょう。先ほど説明したように、失敗確率は失敗分布の正規化定数になっています。これは非常に重要な関係性です。
失敗分布を思い出してください。我々は失敗している軌道の分布を表現しようとしていました。この分布の非正規化密度は、公称分布の密度に失敗かどうかを示す指示関数（indicator function）を掛けたものでした。しかし、この密度を正規化するには、すべての軌道に対してこの非正規化密度を積分する必要があります。その結果得られるのが、我々が計算しようとしている失敗確率に他なりません。
数学的に表現すると：
P(失敗) = ∫ p(τ) × [τが失敗] dτ
ここで p(τ) は公称軌道分布の密度で、[τが失敗] は軌道が失敗の場合に1、そうでない場合に0を返す指示関数です。この積分が正確に失敗確率であり、同時に失敗分布の正規化定数でもあります。
しかし、この正規化定数を直接計算することは一般的に困難です。特に、複雑なシステムでは解析的に計算することはほぼ不可能で、数値的な近似が必要になります。これが前回の講義で、失敗分布からサンプリングする方法について多くの時間を費やした理由です。
失敗分布の正規化定数が失敗確率であるという事実は、我々の問題に対する異なるアプローチを示唆しています。失敗分布からサンプリングする際に開発した技術を利用して、その正規化定数を推定することができるのです。そして、その正規化定数こそが我々が知りたい失敗確率なのです。
この講義では、この正規化定数、つまり失敗確率を推定するための方法を詳しく見ていきます。特に重要度サンプリングというテクニックは、この問題に対する非常に効果的なアプローチを提供します。それは公称分布とは異なる分布からサンプリングすることで、特に稀な失敗事象に対する推定効率を大幅に向上させることができます。
全体として、失敗確率が失敗分布の正規化定数であるという理解は、我々の問題に対する理論的な洞察を提供するだけでなく、実用的な推定アルゴリズムの開発にも直接つながっています。
3. 直接推定法
3.1 ベルヌーイ分布としてのモデル化
Sydney Katz：では、直接推定法の説明に入りましょう。公称システムからサンプルを生成する過程を2段階のプロセスに分けて考えてみます。ここでは特に複雑なことはありません。
第1段階：システムのロールアウトを実行します。例えば、倒立振り子のロールアウトを行います。 第2段階：その結果が失敗であれば1、そうでなければ0とラベル付けします。
これを何度も繰り返すと、1と0のデータセットが得られます。このプロセスは、パラメータp_failによって定義されるベルヌーイ分布からサンプルを生成することと同等だと考えることができます。
これは非常に高度な言い方に聞こえるかもしれませんが、実際には単純なことを言っています。パラメータp_failがあり、システムをサンプリングするたびに、確率p_failで失敗、確率(1-p_fail)で非失敗になります。これはパラメータp_failを持つベルヌーイ分布そのものです。
このように定式化する理由は、問題をパラメータ推定問題として捉えることができるからです。これは第2章でシステムモデリングを行った際に使用した手法と同じです。つまり、第2章で使用したのと同じテクニックを使って、失敗確率を推定することができます。
具体的に、第2章で説明したパラメータ推定の手法は、最尤推定（Maximum Likelihood Estimation）とベイズ推定（Bayesian Estimation）でした。この二つのアプローチを使って、ベルヌーイ分布のパラメータp_fail、つまり失敗確率を推定することができます。
このモデル化の利点は、すでに学んだパラメータ推定の枠組みをそのまま適用できることです。つまり、システムを実行して失敗か非失敗かを記録するという単純なプロセスを、統計的に厳密なパラメータ推定問題として定式化できるのです。
3.2 最尤推定（MLE）アプローチ
Sydney Katz：ベルヌーイ分布のパラメータ推定に最尤推定（Maximum Likelihood Estimation）を適用してみましょう。公称分布からサンプルを取得し、それらが失敗かどうかのデータセットが与えられた場合、ベルヌーイ分布の最尤推定量は直感的に予想される通りのものになります。
もし街中で誰かに声をかけて、「ここにデータがあるから、失敗確率はどのくらいだと思うか？」と尋ねたら、おそらくその人はこう答えるでしょう：失敗した軌道の数を、サンプリングした軌道の総数で割る。これが正に最尤推定量です。
数式で表すと：
p̂_fail = n_fail / m
ここで、n_failは失敗した軌道の数、mはサンプリングした軌道の総数です。
この推定量は非常に直感的で理解しやすいですが、その性質について少し考えてみましょう。推定量を評価する際に重要な特性が3つあります：不偏性（unbiasedness）、一貫性（consistency）、分散（variance）です。
不偏性とは、推定量の期待値が真の値に等しいという性質です。数式で表すと： E[p̂_fail] = p_fail
つまり、もし我々が多くの異なる推定を行い、それらの平均を取ると、真の失敗確率に近づくということです。
一貫性は、サンプル数が無限に増えるにつれて、推定値が真の値に収束する性質です： lim{m→∞} p̂fail = p_fail
最後に、分散は推定値の真の値周りの散らばりを表す量です。これらの特性を視覚的に理解するために、いくつかの例を見てみましょう。
この最尤推定アプローチは非常にシンプルであり、計算も容易です。しかし、後ほど見るように、特に失敗が稀なシステムの場合には効率の問題が生じることがあります。それでも、直接的でわかりやすいため、多くの場合でベースラインとして使用されます。
3.3 推定器の特性（不偏性、一貫性、分散）
Sydney Katz：推定器の特性について、もう少し詳しく見ていきましょう。不偏性、一貫性、分散という3つの特性があります。
不偏性は、推定量の期待値が真の値に等しいという性質です。つまり、もし様々なp_fail推定値の期待値をとると、それは真の失敗確率になるはずです。
一貫性は、サンプル数が無限に近づくにつれて、推定値が真の値に収束するという性質です。つまり、サンプル数mが無限に近づくにつれて、我々の推定値p̂_failは真の失敗確率p_failに等しくなるということです。
分散は、推定値の真の値周りの散らばりを表す量です。
これらの概念を画像で明確にしましょう。ここに3つの異なる推定器の例があります。すべてのプロットにおいて、真の失敗確率はこの点線で示されています。推定値の平均は濃い線で、複数の試行での推定値の散らばりは陰影付きの領域で示されています。横軸は右に行くほどサンプル数が増えていることを表しています。
最初の推定器は不偏ではありません。なぜなら、どの時点でも推定値の期待値（青い線）が真の失敗確率ではないからです。例えばここを見ると、期待値は真の失敗確率の上にあります。つまり、この推定器には偏りがあります。しかし、サンプル数が増えるにつれて真の値に収束しているので、一貫性はあります。
2つ目の推定器は、期待値が真の失敗確率に等しいので不偏です。しかし、サンプル数が増えても分散が減少せず、真の値に収束していないので一貫性がありません。
3つ目の推定器は不偏であり、かつ一貫性もあります。期待値は真の失敗確率を中心に対称に分布しており、サンプル数が増えるにつれて分散が減少しています。
分散については、サンプル数が少ない場合は高い分散があります。つまり、少数のサンプルしか実行しない場合、非常に異なる推定値が得られる可能性があります。一方、多くのサンプルを実行すると、真の失敗確率に近い推定値が得られる可能性が高くなります。
理想的には、推定器は不偏で一貫性があり、低分散であることが望ましいですが、常にそうとは限りません。
モンテカルロ直接推定器のこれらの特性を見てみましょう。この推定器は不偏です。つまり、期待値は真の失敗確率を与えます。そしてこの推定器の分散は、ベルヌーイ分布の分散をサンプル数で割ったものになります：
Var[p̂_fail] = p_fail(1-p_fail)/m
この式から、サンプル数mが無限に近づくと分散はゼロに近づくことがわかります。これは推定器が一貫性を持つことを意味します。分散がゼロになると、無限サンプルの極限では真の値に正確に収束します。
これを経験的に示すには、各試行で5,000サンプルずつ10試行を実行しました。そして各試行を途中で切り取って、例えば最初の1,000サンプルだけを使った場合、失敗確率の平均推定値と、10回の試行間での推定値の散らばりを見ました。予想通り、サンプル数が増えるにつれて推定値の散らばりは小さくなります。
この分散から、次の2つの観察ができます：
- サンプル数が増えるにつれて、推定値の分散は減少します。これは直感的であり、より多くのサンプルを得ることで、より良い推定値が得られるということです。
- 失敗確率が減少するにつれて、推定値の分散は増加します。これは稀な失敗事象を持つシステムにとって悪いニュースです。
このような稀な失敗事象の問題に対処するために、重要度サンプリング技術について説明します。
3.4 稀な失敗イベントに対する課題
Sydney Katz：直接推定法の最大の課題は、稀な失敗イベントを持つシステムに対する効率の悪さです。これは安全クリティカルなシステムにおいて特に重要な問題です。なぜなら、そのようなシステムは設計上、失敗が非常に稀であることが求められるからです。
先ほど示した分散の式から、失敗確率p_failが非常に小さい場合、推定の分散は主にサンプル数mに依存することがわかります。つまり、非常に稀な失敗を正確に推定するためには、膨大な数のサンプルが必要になります。
例えば、失敗確率が10^-9のオーダー（10億分の1）のシステムを考えてみましょう。このような場合、1つの失敗を観測するだけでも10億回のシミュレーションが必要かもしれません。しかし、1つの失敗だけを観測しても、失敗確率の推定に大きな自信を持つことはできません。おそらく正確な推定のためには、10億回よりもさらに多くのシミュレーションが必要でしょう。これは非常に非効率的です。
この問題を具体的な例で見てみましょう。航空機衝突回避システムの失敗確率を推定したいとします。12月の航空安全データベースに飛行記録があり、そこに中空衝突の記録がなかったとします。この場合、最尤推定による失敗確率の推定値は何になるでしょうか？
そうです、ゼロになります。しかし、実際には失敗確率がゼロであるとは考えにくいです。実際、残念ながら先週も事故があったことから、明らかにゼロではありません。
このように、稀な失敗事象を持つシステムに対して直接推定法を適用すると、失敗の観測がない場合、失敗確率をゼロと推定してしまうという問題があります。これは明らかに正しくない推定です。
これらの課題を克服するために、次に説明するベイズ推定法では、失敗確率に関する分布を生成し、ゼロの推定値を避けることができます。さらに重要度サンプリングという技術は、稀な失敗の推定効率を大幅に向上させることができます。
稀な失敗事象に対する推定の難しさは、安全クリティカルなシステムの検証において中心的な課題の一つです。直接推定法の限界を理解することは、より高度な手法の必要性を認識する上で重要です。
4. ベイズ推定法
4.1 ゼロ失敗の問題への対応
Sydney Katz：最尤推定法の問題点として、ゼロ失敗の場合の推定値についてさらに考えてみましょう。航空機衝突回避システムの例では、12月のデータに中空衝突がなかった場合、最尤推定では失敗確率がゼロと推定されてしまいます。しかし、実際には失敗確率がゼロであるとは考えにくいです。
この問題に対処するために、ベイズ推定法を使うことができます。ベイズ推定では、失敗確率に関する分布を生成します。つまり、失敗確率がゼロである可能性もありますが、他の非ゼロの確率値についても、それらの確率を割り当てることができます。
これはどのように行うのでしょうか。フリスビー投げの例を思い出してください。そこでは多くのフリスビーを投げてデータセットを収集しました。θはフリスビーが同じ方向に着地する確率でした。nは同じ方向に着地した回数、mは投げた合計回数でした。
このフリスビーの例を失敗確率の推定に応用できます。フリスビーを投げる代わりに、システムをロールアウトして多くのデータを収集します。θは、フリスビーが同じ方向に着地する確率ではなく、失敗確率p_failになります。nは失敗軌道の数、mはシミュレートした軌道の総数です。
ベイズアプローチでは、最尤推定のように単一の値を推定するのではなく、確率分布として表現します。これにより、ゼロ失敗の場合でも、失敗確率に関する有用な情報を得ることができます。
ベイズ推定の主な利点は、ゼロ失敗の問題に対処できることです。観測された失敗がない場合でも、失敗確率がゼロである可能性は高いかもしれませんが、非ゼロである可能性も排除しません。これは特に安全クリティカルなシステムにおいて重要です。なぜなら、失敗確率がゼロであるというのは通常、非現実的な仮定だからです。
さらに、ベイズアプローチは事前知識を組み込むことができます。例えば、類似システムの過去の経験から、おおよその失敗確率の範囲がわかっている場合、その情報を事前分布として組み込むことができます。
これらの理由から、特に稀な失敗事象を持つシステムに対しては、最尤推定よりもベイズ推定が適していることが多いです。
4.2 ベータ分布による確率モデリング
Sydney Katz：ベイズ推定法では、失敗確率に関する事後分布を求めたいと考えます。つまり、データセットDが与えられたときの確率p_failの分布p(θ|D)を推定します。フリスビー投げの例では、θはフリスビーが同じ方向に着地する確率でしたが、ここでは失敗確率p_failを表します。
フリスビーの例では、θに関する事後分布がベータ分布に従うことを発見しました。失敗確率の推定でも同様のアプローチが使えます。フリスビーを投げる代わりにシステムのロールアウトを行い、同じ方向への着地回数の代わりに失敗軌道の数をカウントします。総フリップ回数は、シミュレートした軌道の総数に対応します。
ベータ分布は、二項分布（ベルヌーイ試行の繰り返し）のパラメータに対する共役事前分布として知られています。これは計算上の利点があり、事後分布も同じベータ分布の形になります。
具体的には、n回の失敗と(m-n)回の非失敗が観測されたとき、一様な事前分布を仮定すると、失敗確率p_failの事後分布はBeta(n+1, m-n+1)となります。ここでm-nは非失敗の回数です。
このベータ分布の意味を理解しましょう。パラメータは、成功（この場合は失敗）の回数+1と、失敗（この場合は非失敗）の回数+1です。例えば、50回のシミュレーションで8回の失敗があった場合、事後分布はBeta(9, 43)となります。
パラメータの+1は、一様な事前分布Beta(1, 1)を使用した結果です。異なる事前知識を持っている場合、異なるパラメータの事前分布を選択できます。
ベータ分布を使った確率モデリングの美しさは、新しいデータが得られるたびに、事後分布を簡単に更新できることです。現在の事後分布を次のステップの事前分布として使用すれば、逐次的にベイズ更新を行うことができます。
このアプローチにより、特にゼロ失敗を観測した場合でも、失敗確率に関する豊かな情報を持つ分布が得られます。分布全体を得ることで、単一の点推定ではなく、確率の範囲やその確信度について議論することができます。
4.3 信頼度の定量化方法
Sydney Katz：ベイズ推定の大きな利点の一つは、失敗確率に関する信頼度を定量化できることです。これにより、単なる点推定ではなく、より豊かな情報を得ることができます。
具体的な例で考えてみましょう。倒立振子をシミュレーションし、50回の軌道中に8回の失敗があったとします。最尤推定では、失敗確率は0.16（8/50）となります。しかし、ベイズアプローチを使うと、失敗確率の分布が得られます。この場合、事後分布はBeta(9, 43)となります。
では、知覚ノイズが非常に小さくなり、失敗がより稀になったケースを考えてみましょう。この場合、失敗を観測しないかもしれません。最尤推定では失敗確率はゼロと推定されますが、ベイズアプローチでは失敗確率の分布を推定します。この分布からは、確率がゼロである可能性が高いですが、もう少し高い確率である可能性も排除されていないことがわかります。
さらに多くの軌道をシミュレートして失敗が観測されなければ、この分布はどんどん先鋭化していきます。つまり、ゼロに近い確率値にますます集中していきます。
ベイズ推定の素晴らしい点は、このような分布に基づいて追加の主張ができることです。例えば、「失敗確率が0.01未満である確率は何か？」という問いに答えることができます。システムを展開する際に、失敗確率が1%未満であることを確認したい場合、この質問は非常に重要です。
このような問いに答えるには、事後分布の0.01より左側の確率質量を計算します。これは、累積分布関数（CDF）を使って計算できます。例えば、「失敗確率が0.01未満である確率は0.637」というような結果が得られます。
あるいは、「失敗確率がある値X未満であると95%確信できる値Xは何か」という問いにも答えることができます。この場合、分布の左側に95%の確率質量がある点を探します。例えば、「失敗確率が0.029未満であると95%確信できる」という結果が得られるかもしれません。
これらの定量的な信頼度の表現は、安全クリティカルなシステムの意思決定において非常に価値があります。例えば、まだ失敗を観測していないが、失敗確率がさらに低いことを証明したい場合、サンプル数を増やすことでこれらの信頼度の指標を改善できます。
サンプル数を50から150に増やしても失敗が観測されなかった場合、「失敗確率が0.01未満である確率」は0.637から0.78に向上します。同様に、95%信頼上限も改善されます。
このように、ベイズ推定は直接的な点推定よりも豊かな情報を提供し、特に稀な失敗事象を持つシステムの評価において強力なツールとなります。
4.4 ベイズ推定の実例と利点
Sydney Katz：実際のベイズ推定の適用例を見てみましょう。ここでは倒立振子システムのシミュレーションを行っています。50回の試行中に8回の失敗が観測された場合、最尤推定では失敗確率は約0.16（8/50）となります。
ここで知覚ノイズを非常に小さくし、失敗がより稀になるようにシステムを調整したとします。この場合、例えば50回の試行で失敗が全く観測されないかもしれません。最尤推定では失敗確率はゼロとなりますが、これは現実的ではありません。
ベイズアプローチを使うと、このようなケースでも失敗確率に関する分布が得られます。この分布を見ると、確率がゼロに近い値である可能性が高いものの、少し高い値である可能性も排除されていないことがわかります。視覚的には、分布がゼロ付近に集中していますが、少し右側にも広がりを持っています。
さらに多くの軌道をシミュレートしても失敗が観測されない場合、この分布はどんどん鋭くなっていきます。つまり、ゼロに近い確率値にますます強く集中していきます。
ベイズ推定の強力な点は、失敗確率に関する様々な主張を定量的に評価できることです。例えば、「失敗確率が0.01未満である確率」を計算できます。この例では、その確率は0.637です。これは、分布の0.01より左側の確率質量を計算することで得られます。
また、「失敗確率がある値X未満であると95%確信できる値X」も計算できます。この場合、分布の左側に95%の確率質量がある点を見つけます。この例では、その値は0.029です。つまり、「失敗確率が0.029未満であると95%確信できる」と言えます。
これらの指標をさらに改善したい場合、より多くのサンプルを収集することができます。試行回数を50から150に増やしても失敗が観測されなかった場合、「失敗確率が0.01未満である確率」は0.637から0.78に向上します。同様に、95%信頼上限も改善されます。
ベイズ推定の主な利点をまとめると：
- ゼロ失敗の問題に対処できる：失敗が観測されない場合でも、失敗確率に関する有用な情報を提供します。
- 信頼度を定量化できる：確率分布として表現することで、様々な主張に対する確信度を計算できます。
- 逐次的に更新できる：新しいデータが得られるたびに、事後分布を簡単に更新できます。
- 事前知識を組み込める：システムに関する事前の知識や経験を、事前分布の形で組み込むことができます。
これらの特性により、ベイズ推定は特に稀な失敗事象を持つ安全クリティカルなシステムの評価において、非常に価値のあるアプローチとなります。最尤推定のような単純な点推定よりも、より豊かな情報を提供し、より堅牢な意思決定をサポートします。
5. 重要度サンプリング
5.1 直接推定法との比較
Sydney Katz：ここまで見てきた直接推定法とベイズ推定法には、稀な失敗事象を持つシステムでの効率の悪さという共通の問題があります。そこで、重要度サンプリング（Importance Sampling）という手法を紹介します。これは直接推定法の改良版と考えることができます。
まず、直接推定法と重要度サンプリングを比較してみましょう。直接推定法では、公称軌道分布からサンプルを多数抽出し、失敗サンプルの数を総サンプル数で割ることで失敗確率を推定します。この方法の問題点は何でしょうか？これはかなり修辞的な質問で、このコースで何度も出てきましたが、稀な失敗事象を持つシステムでは非効率的だということです。
例えば、失敗確率が10^-9のオーダー（10億分の1）のシステムを考えてみましょう。このような場合、1つの失敗を観測するだけでも10億回のシミュレーションが必要かもしれません。しかし、1つの失敗だけを観測しても、失敗確率の推定に大きな自信を持つことはできません。おそらく正確な推定のためには、10億回よりもさらに多くのシミュレーションが必要でしょう。これは非常に非効率的です。
この問題に対処するために重要度サンプリングが登場します。重要度サンプリングの主な違いは、公称軌道分布からサンプルを抽出する代わりに、「提案分布」と呼ばれる別の分布からサンプルを抽出することです。この提案分布は、失敗軌道をより頻繁に生成するように選ばれます。
これは棄却サンプリングの提案分布と似たような概念です。重要度サンプリングの提案分布を適切に選ぶことで、より多くの失敗サンプルを観測できるようになります。これは稀な失敗の問題に対処する上で大きな助けとなります。
しかし、ここで注意が必要です。提案分布からサンプルを抽出して、単純に失敗の数を総サンプル数で割ることは、正しい失敗確率の推定にはなりません。なぜなら、実際の世界では公称分布からのサンプルが観測されるからです。
しかし、驚くべきことに、提案分布からのサンプルを適切に「再重み付け」することで、正しい失敗確率の推定が可能になります。これが重要度サンプリングの核心部分です。
直接推定法と重要度サンプリングの主な違いは：
- サンプリング分布：直接推定法は公称分布からサンプリングし、重要度サンプリングは提案分布からサンプリングします。
- 計算方法：直接推定法は単純に失敗サンプルの比率を計算し、重要度サンプリングはサンプルを再重み付けします。
- 効率性：直接推定法は稀な失敗に対して非効率的ですが、重要度サンプリングは適切な提案分布を使用することで大幅に効率を向上させることができます。
この重要度サンプリングの概念は、このクラスでの最も基本的かつ重要な概念の一つです。次のセクションでは、この手法の数学的基礎と実装方法について詳しく見ていきます。
5.2 提案分布からのサンプリング
Sydney Katz：重要度サンプリングの最初のステップは、「提案分布」からサンプルを抽出することです。この提案分布は、失敗軌道をより頻繁に生成するように選ばれます。ここでは、この提案分布をQと呼ぶことにします。
提案分布Qの選択が重要度サンプリングの性能を大きく左右します。理想的な提案分布は、失敗が発生しやすい領域により高い確率を割り当てるものです。これにより、限られたサンプル数でも多くの失敗サンプルを観測できるようになります。
提案分布の選択について視覚的に説明しましょう。この図では直接推定法と重要度サンプリングを比較しています。考えているのは、ある閾値より下（赤い領域）が失敗となる単純なガウスシステムです。
直接推定の場合、公称分布からサンプルを抽出し、失敗サンプルの数を総サンプル数で割ります。この図のx軸はサンプル数、y軸は失敗確率の推定誤差の絶対値を示しています。理想的には、この値は0に近いほど良く、分散も小さいほど良いです。
直接推定では、かなり大きな分散と誤差が見られます。500サンプルを使っても、推定値は真の失敗確率からかなり離れていることがあります。
重要度サンプリングでは、提案分布を変更することで、同じサンプル数でもより良い推定が可能になります。右側の図で提案分布を示しています。最初は公称分布と同じで、つまり直接推定と同じ結果になります。
しかし、提案分布の平均を左側（失敗領域側）に移動させると、失敗サンプルがより多く得られるようになります。この図の青い線は重要度サンプリングの推定誤差と分散を示しています。直接推定と比較して、分散が大幅に減少し、誤差も小さくなっています。
有効サンプルサイズという概念も有用です。直接推定の場合、有効サンプルサイズは観測された失敗の数です。提案分布を変更すると、より多くの失敗を観測できるようになりますが、これらの失敗は公称分布の下での正しい重みを持っていないため、有効サンプルサイズは観測された失敗の数よりも少なくなります。しかし、適切な提案分布を選べば、直接推定よりも大幅に有効サンプルサイズを増やすことができます。
提案分布をさらに調整して、失敗閾値に中心を置くと、さらに多くの失敗を観測でき、有効サンプルサイズが増加します。また、提案分布の分散を減らして、より多くのサンプルが失敗領域に入るようにすることもできます。
しかし、提案分布の選択には注意が必要です。失敗を多くサンプリングするだけでなく、「可能性の高い」失敗を「可能性の低い」失敗よりも頻繁にサンプリングすることが重要です。つまり、提案分布は失敗サンプルを増やすだけでなく、公称分布の下で確率の高い失敗をより多くサンプリングするべきです。
提案分布の選択は重要度サンプリングの性能を決定する重要な要素であり、次のセクションでは、サンプルをどのように再重み付けするかについて説明します。
5.3 重み付け手法の導出
Sydney Katz：重要度サンプリングの核心部分は、提案分布からのサンプルを適切に再重み付けして、正しい失敗確率の推定を得る方法です。この再重み付けの数学的導出を見ていきましょう。
まず、失敗確率の定義から始めます。失敗確率は、公称分布からの軌道が失敗かどうかの期待値として定義されます：
p_fail = E_{τ〜P}[τが失敗]
これを積分の形で書くと：
p_fail = ∫ p(τ) · [τが失敗] dτ
ここでp(τ)は公称密度で、[τが失敗]は軌道が失敗なら1、そうでなければ0を返す指示関数です。
ここで一見奇妙に思えるかもしれない操作をします。この式を1倍します。ただし、1の表現としてq(τ)/q(τ)を選びます。つまり：
p_fail = ∫ p(τ) · [τが失敗] · (q(τ)/q(τ)) dτ
これにより、提案分布q(τ)が式に導入されました。次に、式の要素を整理して：
p_fail = ∫ (p(τ)/q(τ)) · [τが失敗] · q(τ) dτ
この式を見ると、提案分布q(τ)からのサンプリングに基づく期待値の形になっていることがわかります：
p_fail = E_{τ〜Q}[(p(τ)/q(τ)) · [τが失敗]]
ここで重要な洞察が得られました。公称分布Pからサンプリングする代わりに、提案分布Qからサンプリングして、各サンプルを(p(τ)/q(τ))で重み付けすることで、同じ期待値を計算できるのです。
期待値を計算するには、多数のサンプルの平均を取ればよいので、失敗確率の推定値p̂_failは次のようになります：
p̂fail = (1/m) · ∑{i=1}^m (p(τi)/q(τi)) · [τ_i が失敗]
ここでτiは提案分布Qからのサンプルです。式を整理するために、(p(τi)/q(τ_i))をW_iと呼ぶことにします。すると：
p̂fail = (1/m) · ∑{i=1}^m W_i · [τ_i が失敗]
ここでW_iは「重要度重み」と呼ばれ、公称分布の下でのサンプルの相対的な尤度を表します。
驚くべきことに、この推定量は真の失敗確率の不偏推定量であり、一貫性も持っています。つまり、提案分布からサンプリングしているにもかかわらず、サンプル数が増えるにつれて真の失敗確率に収束するのです。
直接推定と重要度サンプリングを比較すると：
- 直接推定：Pからサンプリングし、p̂fail = (1/m) · ∑{i=1}^m [τ_i が失敗]
- 重要度サンプリング：Qからサンプリングし、p̂fail = (1/m) · ∑{i=1}^m W_i · [τ_i が失敗]
つまり、重要度サンプリングでは、サンプルが公称分布から来なかったため、すべてのサンプルを等しく扱うことはできません。真の分布での見え方に応じて、サンプルを再重み付けする必要があるのです。
この重み付け手法は重要度サンプリングの本質であり、稀な失敗事象の推定効率を大幅に向上させることができます。
5.4 不偏推定器としての証明
Sydney Katz：重要度サンプリングが提供する推定量が、真の失敗確率の不偏推定量であるという事実は非常に重要です。ここでは、その証明と推定器の実装について詳しく見ていきます。
重要度サンプリングの推定量は次のように表されます：
p̂fail = (1/m) · ∑{i=1}^m W_i · [τ_i が失敗]
ここでW_i = p(τi)/q(τi)は重要度重みです。この推定量が不偏であることを示すには、その期待値が真の失敗確率p_failに等しいことを証明する必要があります。
期待値の定義から：
E[p̂fail] = E[(1/m) · ∑{i=1}^m W_i · [τ_i が失敗]]
サンプルは互いに独立なので：
E[p̂fail] = (1/m) · ∑{i=1}^m E[W_i · [τ_i が失敗]]
提案分布Qからのサンプリングに基づく期待値なので：
E[W_i · [τ_i が失敗]] = ∫ (p(τ)/q(τ)) · [τが失敗] · q(τ) dτ
これを簡略化すると：
E[W_i · [τ_i が失敗]] = ∫ p(τ) · [τが失敗] dτ = p_fail
よって：
E[p̂fail] = (1/m) · ∑{i=1}^m p_fail = p_fail
これにより、重要度サンプリング推定量が不偏であることが証明されました。
また、サンプル数が無限大に近づくにつれて、推定量は真の失敗確率に収束するという一貫性も持ちます。つまり、提案分布からサンプリングしているにもかかわらず、理論的には真の分布における正確な推定が可能なのです。
この推定器をコードで実装すると次のようになります：
function importance_sampling(m, proposal_dist)
samples = zeros(m)
weights = zeros(m)
failures = zeros(m)
for i in 1:m
# 提案分布からサンプリング
samples[i] = rand(proposal_dist)
# 重要度重みを計算
weights[i] = pdf(nominal_dist, samples[i]) / pdf(proposal_dist, samples[i])
# 失敗かどうかをチェック
failures[i] = is_failure(samples[i])
end
# 失敗確率の推定値を計算
p_fail_estimate = sum(weights .* failures) / m
return p_fail_estimate
endこのコードでは、提案分布からサンプルを抽出し、各サンプルの重要度重みを計算し、失敗かどうかをチェックします。最後に、重み付けされた失敗の合計を総サンプル数で割って、失敗確率の推定値を返します。
重要度サンプリングの美しさは、適切な提案分布を選ぶことで、同じサンプル数でも直接推定よりもはるかに良い推定値が得られることです。特に稀な失敗事象を持つシステムでは、この手法は推定効率を大幅に向上させることができます。
また、重要度サンプリングの不偏性と一貫性は、単なる近似手法ではなく、理論的に正当化された手法であることを示しています。これにより、安全クリティカルなシステムの検証において、信頼性の高い失敗確率の推定が可能になります。
6. 最適な提案分布
6.1 分散最小化の観点
Sydney Katz：重要度サンプリングの効率は提案分布の選択に大きく依存します。そこで自然に生じる疑問は、「最適な提案分布は存在するのか？」ということです。この問いに答えるためには、まず「最適」の意味を定義する必要があります。
重要度サンプリングの文脈では、最適な提案分布とは、推定器の分散を最小化するものと定義できます。分散が小さいほど、少ないサンプル数でも信頼性の高い推定が可能になるからです。
重要度サンプリング推定器の分散は次のように表されます（この導出は教えなくてもよいですが、参考として示します）：
Var[p̂_fail] = (1/m) · (∫ (p(τ)²/q(τ)) · [τが失敗]² dτ - p_fail²)
この分散をq(τ)に関して最小化することが、最適な提案分布を見つけることになります。
特定の提案分布を試してみましょう。最適な提案分布の候補として、次のようなものを考えてみます：
q(τ) = (p(τ) · [τが失敗]) / p_fail
この候補の提案分布を分散の式に代入してみると：
Var[p̂_fail] = (1/m) · (∫ (p(τ)² · [τが失敗]² / (p(τ) · [τが失敗] / p_fail)) dτ - p_fail²)
整理すると：
Var[p̂_fail] = (1/m) · (∫ p(τ) · [τが失敗] · p_fail dτ - p_fail²)
さらに：
Var[p̂_fail] = (1/m) · (p_fail · ∫ p(τ) · [τが失敗] dτ - p_fail²)
積分部分はp_failなので：
Var[p̂_fail] = (1/m) · (p_fail · p_fail - p_fail²) = 0
驚くべきことに、この提案分布を使用すると、分散がゼロになります！これは、理論上、単一のサンプルだけで真の失敗確率を完全に正確に推定できることを意味します。
しかし、ここで重要な疑問が生じます。誰かが「失敗をまったくサンプリングしない提案分布を使えば、推定値の分散も常にゼロになるのではないか？」と尋ねるかもしれません。確かに、そのような分布を使うと常に失敗確率をゼロと推定することになり、分散はゼロになります。しかし、それは明らかに有用な推定ではありません。
分散を最小化するだけでなく、有用な推定を得るためには、失敗をサンプリングする必要があります。そして上記の候補分布は、分散を最小化しながらも、失敗をサンプリングする分布なのです。
次のセクションでは、この最適な提案分布が実は私たちがよく知るものであることを見ていきます。
6.2 失敗分布との関係
Sydney Katz：先ほど導出した最適な提案分布を改めて見てみましょう：
q(τ) = (p(τ) · [τが失敗]) / p_fail
この式をよく見ると、何か見覚えがあるかもしれません。そうです、これは失敗分布の正確な表現です！実は、分散を最小化する最適な提案分布は、失敗分布そのものなのです。
振り返ってみると、失敗分布は次のように定義されます：
- 非正規化密度：p(τ) · [τが失敗]
- 正規化定数：p_fail
したがって、正規化された失敗分布の密度は (p(τ) · [τが失敗]) / p_fail となり、これは先ほど導出した最適な提案分布と一致します。
この結果は直感的にも理解できます。失敗分布は、失敗する軌道のみに確率を割り当て、かつ公称分布における相対的な確率に比例して確率を割り当てます。つまり、公称分布で発生しやすい失敗に対してより高い確率を割り当てるのです。
理論的には、この最適な提案分布（つまり失敗分布）を使用すると、重要度サンプリングの推定量の分散はゼロになります。つまり、単一のサンプルで真の失敗確率を完全に正確に推定できるのです。
失敗分布からサンプルを抽出すると、すべてのサンプルが失敗であり、かつその重要度重みがすべて同じになります。具体的には、すべてのサンプルτについて：
W = p(τ)/q(τ) = p(τ)/((p(τ) · [τが失敗]) / p_fail) = p_fail / [τが失敗]
サンプルはすべて失敗であるため、[τが失敗] = 1であり、したがってW = p_failとなります。
つまり、失敗分布からサンプルを抽出すると、各サンプルの重要度重みは正確に失敗確率になります。これにより、1つのサンプルから直接失敗確率を読み取ることができ、推定量の分散はゼロになります。
この結果は理論的に美しく、重要度サンプリングの基礎を深く理解する上で重要です。しかし、次のセクションで見るように、実践においてはこの知見を直接適用することは難しいという問題があります。
6.3 円環的な問題（失敗確率が必要）
Sydney Katz：最適な提案分布が失敗分布であることがわかりましたが、ここで大きな問題に直面します。最適な提案分布を表現すると：
q(τ) = (p(τ) · [τが失敗]) / p_fail
この式には正規化定数としてp_fail（失敗確率）が含まれています。しかし、p_failは私たちが推定しようとしている値そのものです！つまり、失敗確率を推定するために最適な提案分布を使いたいのですが、その提案分布を定義するためには失敗確率が必要という循環的な問題が生じるのです。
これは非常に皮肉な状況です。最適な提案分布を知るためには、すでに答えを知っていなければならないのです。前回の講義で学んだように、失敗分布からのサンプリングは可能ですが、その正規化密度を計算することができません。そして正規化密度の計算には、正規化定数である失敗確率の知識が必要です。
質問がありました。「失敗分布からサンプルを抽出している場合、重要度重みは正確に失敗確率に等しくなるのではないか？」
その通りです。失敗分布からサンプリングすると、すべてのサンプルの重要度重みは失敗確率p_failに等しくなります。これは直観的に見ると、サンプル一つから失敗確率を直接読み取れることを意味します。これにより推定量の分散がゼロになるのです。しかし、重要度重みを計算するためには正規化密度が必要で、それには失敗確率の知識が必要となります。
もう一つの質問がありました。「棄却サンプリングなどの方法を使って、サンプルが棄却サンプリングから出てくる確率を何らかの方法で定量化できないのか？」
実際には、失敗分布からサンプリングできれば、その重要度重みは正確に失敗確率に等しくなります。しかし、正規化密度を計算するためには失敗確率が必要で、これが循環的な問題を引き起こします。この問題を理解するには少し時間がかかるかもしれません。
このジレンマにより、実践では最適な提案分布を直接使用することはできません。代わりに、次のセクションで説明するように、可能な限り失敗分布に近い提案分布を選択することを目指します。
6.4 提案分布の選択基準
Sydney Katz：最適な提案分布（失敗分布）を直接使用できないという問題に直面しているため、実践的なアプローチとしては、可能な限り失敗分布に近い提案分布を選択することを目指します。提案分布を選択する際の重要な基準について説明します。
まず、提案分布を選ぶ際に考慮すべき主な要素は以下の通りです：
- 失敗分布との近さ：分散を最小化するためには、提案分布が失敗分布に近いほど良いです。つまり、公称分布において可能性の高い失敗に対してより多くの確率質量を割り当てる分布が望ましいです。
- 計算可能性：提案分布からサンプリングできるだけでなく、その正規化密度も計算できる必要があります。これは重要度重みW = p(τ)/q(τ)を計算するために不可欠です。
- サポートの包含：提案分布のサポート（非ゼロ確率を持つ領域）は、元の分布のサポートを含む必要があります。特に、失敗が発生する可能性のある領域すべてに非ゼロの確率を割り当てる必要があります。
これらの基準を満たす提案分布を設計するためのいくつかのアプローチがあります：
まず、失敗分布に似た、既知の計算可能な分布を選ぶことができます。例えば、多変量ガウス分布などのパラメトリックな分布を使用し、そのパラメータを調整して失敗領域に高い確率を割り当てることができます。
先ほどの質問に関連して、提案分布が失敗分布のサポートを完全にカバーする必要があります。棄却サンプリングでは、提案分布が目標分布を「上から押さえる」必要がありましたが、重要度サンプリングではそのような制約はありません。ただし、目標分布が非ゼロ確率を持つ領域では、提案分布も非ゼロ確率を持つ必要があります。
例として、前回の講義で見たバイモーダル（二峰性）の失敗分布を考えてみましょう。このような場合、提案分布も両方のモードをカバーする必要があります。前回説明したスムージング技術は、このような複雑な失敗分布からのサンプリングに役立ちますが、重要度サンプリングのためには、計算可能な正規化密度を持つ分布に近似する必要があります。
次のセクションでは、失敗分布からサンプルを取得し、それに基づいて計算可能な提案分布を構築する実践的なアプローチについて説明します。これにより、理論的な最適性と実践的な実装可能性のバランスを取ることができます。
7. 実践的アプローチ
7.1 失敗分布からのサンプル取得
Sydney Katz：最適な提案分布が失敗分布であることはわかりましたが、その正規化密度を直接計算することはできません。そこで実践的なアプローチとして、まず失敗分布からサンプルを取得し、それを基に計算可能な提案分布を構築する方法を考えます。
失敗分布からのサンプル取得には、前回の講義で説明した手法を使用できます。具体的には以下の手法があります：
- マルコフ連鎖モンテカルロ（MCMC）法：非正規化密度からサンプリングするための強力な手法です。前回見たように、スムージングを適用してモード間の移動を促進することで、多様なサンプルを取得できます。
- 棄却サンプリング：失敗領域を含む提案分布を設計し、失敗していないサンプルを棄却する方法です。効率は提案分布の選択に大きく依存します。
- ハイブリッド手法：スムージングを適用したMCMCと棄却サンプリングを組み合わせるなど、複数の手法を組み合わせることもできます。
これらの手法を使って、まず失敗分布からの代表的なサンプル集合を取得します。このステップの目標は、可能な限り多様で代表的な失敗サンプルを収集することです。特に複数の失敗モードがある場合、それらすべてを適切に表現するサンプルが必要です。
例えば、倒立振子システムでは、左に倒れる失敗と右に倒れる失敗の両方からサンプルを取得することが重要です。同様に、航空機衝突回避システムでは、異なるタイプの衝突シナリオ（正面衝突、斜め衝突など）からサンプルを取得する必要があります。
実践において、失敗分布からのサンプリングは計算コストが高い場合があります。特に複雑なシステムや高次元の状態空間を持つシステムでは、代表的なサンプルを得るために大量の計算が必要になることがあります。しかし、このステップに投資する価値はあります。なぜなら、良い提案分布を構築することで、その後の重要度サンプリングでの失敗確率推定の効率が大幅に向上するからです。
また、サンプリング効率を向上させるために、システムに関するドメイン知識を活用することも重要です。例えば、特定のパラメータが失敗に大きく影響することがわかっている場合、そのパラメータに焦点を当ててサンプリングすることで、効率的に失敗サンプルを収集できます。
次のステップでは、これらの失敗サンプルを使って、計算可能な密度を持つ分布にフィッティングする方法を見ていきます。
7.2 正規化密度計算可能な分布へのフィッティング
Sydney Katz：前のセクションで失敗分布からサンプルを取得したら、次のステップはこれらのサンプルを使って、正規化密度を計算できる分布にフィッティングします。このアプローチでは、失敗分布からサンプルを取得する方法は知っていても、その正規化密度を計算できないという問題を解決します。
具体的には、失敗分布からのサンプルを使って、正規化密度を計算できる既知の分布族（例えばガウス分布など）にフィッティングします。この分布族は、次の条件を満たす必要があります：
- 計算可能な正規化密度：重要度重みW = p(τ)/q(τ)を計算するために、提案分布q(τ)の正規化密度を計算できる必要があります。
- サンプリング可能：提案分布からサンプルを容易に抽出できる必要があります。
- 表現力：失敗分布の主要な特徴（例えば複数のモードなど）を適切に表現できる必要があります。
フィッティングには様々な方法がありますが、一般的なアプローチは最尤推定（MLE）や最小二乗法などの標準的な統計的手法を使用することです。例えば、失敗サンプルをガウス分布にフィッティングする場合、サンプルの平均と分散を計算し、それをガウス分布のパラメータとして使用します。
より複雑な失敗分布の場合は、単一のガウス分布では十分な表現ができないことがあります。そのような場合は、混合ガウス分布などのより柔軟なモデルを使用できます。混合ガウス分布は、複数のガウス分布の重み付き和であり、複数のモードを持つ分布を表現することができます。これは、バイモーダルやマルチモーダルな失敗分布に特に有用です。
フィッティングの品質は、得られるサンプルの代表性と数に大きく依存します。十分な数のサンプルが得られない場合、フィッティングされた分布は失敗分布の一部の側面しか捉えられないかもしれません。したがって、可能な限り多くの代表的なサンプルを収集することが重要です。
また、分布族の選択も重要です。例えば、失敗分布が強い非対称性を持つ場合、対称なガウス分布では適切に表現できない可能性があります。そのような場合は、ログ正規分布や他の非対称分布を検討することがあります。
フィッティングの評価には、実際にフィッティングした分布からサンプリングし、元の失敗サンプルとの類似性を視覚的に比較することや、適合度の統計的指標を使用することができます。
重要なのは、フィッティングされた分布が失敗分布の重要な特性を捉えていることであり、完璧な一致を求める必要はありません。提案分布が失敗分布に近いほど、重要度サンプリングの効率は向上しますが、計算可能性と表現力のバランスを取ることが重要です。
7.3 ガウス分布の活用例
Sydney Katz：実際にガウス分布を使った提案分布の構築例を見てみましょう。ガウス分布は計算が容易で柔軟性があるため、提案分布として広く使われています。ここでは、Plutoノートブックを使った実装例を紹介します。
このノートブックでは、まず失敗分布からサンプルを取得します。前回の講義で説明したMCMCサンプリングを使って、単純なガウス分布のシステムから失敗サンプルを抽出しています。このシステムでは、ある閾値より下の領域が失敗と定義されています。
次に、取得したサンプルに対してガウス分布をフィッティングします。このノートブックでは、単純にfit関数を呼び出してサンプルに正規分布をフィッティングしています。
# 失敗分布からサンプルを取得（MCMCを使用）
failure_samples = sample_from_failure_distribution(n_samples)
# ガウス分布をフィッティング
proposal_dist = fit(Normal, failure_samples)グラフで結果を見ると、赤い線が真の失敗分布で、その上にプロットされているのがフィッティングしたガウス分布です。サンプルもプロットされており、フィッティングの品質を視覚的に確認できます。
このフィッティングされたガウス分布を提案分布として使用した重要度サンプリングは、直接推定法と比較してどれだけ良い性能を示すでしょうか。実験結果を見ると、同じサンプル数でも重要度サンプリングの方が明らかに優れた性能を示しています。
特に、真の失敗分布が単一のモードを持ち、比較的単純な形状である場合、ガウス分布は非常に良い近似を提供します。しかし、真の失敗分布が複数のモードを持つ場合や、複雑な形状を持つ場合は、単一のガウス分布では十分に表現できないことがあります。
そのような場合は、複数のガウス分布の混合モデルを使うことができます。例えば、前回見たバイモーダルな失敗分布の場合、2つのガウス分布の混合モデルを使うことで、両方のモードを適切に表現できます。
# 混合ガウスモデルをフィッティング
gmm = fit(MixtureModel, failure_samples, 2) # 2つのコンポーネントガウス分布を使用する利点としては、以下が挙げられます：
- 計算効率が高い：正規化密度の計算やサンプリングが効率的に行えます。
- 汎用性がある：多くのシステムの失敗分布は、ガウス分布や混合ガウス分布で十分に近似できることが多いです。
- 実装が簡単：多くの言語やライブラリでガウス分布のサポートが充実しています。
ただし、ガウス分布は対称性があり、裾が指数関数的に減少するという特性があります。一方、実際の失敗分布は非対称であったり、裾が重い（つまり極端な値の確率が高い）場合があります。そのような場合は、t分布や他の裾の重い分布の使用を検討する必要があるかもしれません。
最終的に、どの分布族を選ぶかは、具体的なシステムの特性や得られる失敗サンプルの性質によって異なります。目標は、計算可能性を維持しながら、失敗分布をできるだけ正確に表現することです。
8. 高度な手法への展望
8.1 適応的重要度サンプリング
Sydney Katz：ここまで説明した重要度サンプリングの基本的なアプローチには、さらなる改良が可能です。次の講義で詳しく説明する予定ですが、その一つが適応的重要度サンプリング（Adaptive Importance Sampling）です。
適応的重要度サンプリングの基本的なアイデアは、サンプリングの過程で得られる情報を活用して、提案分布を動的に改良していくというものです。従来の重要度サンプリングでは、提案分布は最初に固定され、その後変更されません。一方、適応的アプローチでは、サンプリングが進むにつれて提案分布のパラメータを更新していきます。
具体的には、以下のようなプロセスで実行されます：
- 初期の提案分布q_0(τ)を選択します。これは前に説明したように、ガウス分布などの計算可能な分布であり、失敗領域に重点を置いたものであることが望ましいです。
- この提案分布からいくつかのサンプルを抽出し、重要度重みを計算します。
- これらのサンプルと重みを使用して、提案分布のパラメータを更新します。例えば、ガウス分布を使用している場合、重み付き平均と分散を計算して、新しいガウス分布のパラメータとして使用します。
- 更新された提案分布からさらにサンプルを抽出し、このプロセスを繰り返します。
このアプローチの利点は、サンプリングが進むにつれて提案分布が徐々に真の失敗分布に近づく可能性があることです。特に、初期の提案分布が失敗分布から大きく離れている場合に効果的です。
また、適応的重要度サンプリングは、複数の失敗モードを持つシステムや、失敗分布の形状が事前にわかりにくいシステムにも有効です。サンプリングが進むにつれて、これまで探索されていなかった失敗モードを発見し、提案分布を調整してそれらのモードもカバーするようになります。
ただし、適応的アプローチには注意点もあります。提案分布の更新方法によっては、特定のサンプリング領域に過度に集中してしまい、他の重要な領域を見落とす可能性があります。これを防ぐためには、探索と活用のバランスを取ることが重要です。
また、適応的重要度サンプリングでは、サンプル間に依存関係が生じるため、理論的な解析が複雑になることがあります。特に、推定量の不偏性や一貫性の証明が、標準的な重要度サンプリングよりも難しくなります。
これらの課題にもかかわらず、適応的重要度サンプリングは実践において非常に強力なツールであり、特に複雑な失敗分布を持つシステムに対して大きな利点をもたらします。これについては次回の講義でさらに詳しく説明しますが、プロジェクト2に取り組む際に有用な技術となるでしょう。
8.2 確率的プログラミングの活用
Sydney Katz：失敗確率の推定をさらに効率化するもう一つのアプローチとして、確率的プログラミング（Probabilistic Programming）の活用があります。これは前回のパラメータ推定の際にも少し触れましたが、確率モデルをコンピュータプログラムとして指定できるようにするパラダイムです。
確率的プログラミングの主な特徴は、確率モデルを明示的にコードとして表現し、そのモデルに対する推論を言語やライブラリが自動的に処理してくれることです。これは失敗確率の推定において非常に強力なツールとなります。
確率的プログラミングを失敗確率の推定に活用する際の大きな利点は、複雑な推論アルゴリズムを自分で実装する必要がないことです。例えば、マルコフ連鎖モンテカルロ（MCMC）やハミルトニアンモンテカルロ（HMC）などの高度なサンプリング手法を、確率的プログラミング言語やライブラリが自動的に提供してくれます。
このクラスでは、特にJulia言語のTuring.jlライブラリを使った例を紹介しています。これは高度なトピックであり、クイズには出ませんが、失敗確率の推定をさらに深く探求したい方には非常に興味深いアプローチです。
確率的プログラミングを用いた失敗確率の推定の基本的なステップは以下の通りです：
- システムのモデルを確率的プログラムとして表現します。これには、初期状態分布、摂動分布、システムのダイナミクスなどが含まれます。
- 失敗の定義を確率的プログラムの中に組み込みます。例えば、特定の状態変数が閾値を超えた場合に失敗と定義できます。
- 確率的プログラミング言語やライブラリの推論機能を使用して、失敗確率を推定します。多くの場合、これはMCMCやその変種を使用して、事後確率分布からサンプリングすることを意味します。
教科書のアルゴリズムを見ると、Turing.jlでモデルを定義する方法が示されています。基本的にはmodel関数を定義し、システムのロールアウトをモデル化します。初期状態が初期状態分布からサンプリングされ、各摂動が摂動分布からサンプリングされるという形です。そして、失敗分布のためのスムーズ化された密度を、サンプリングしたい分布の密度としてTuringに伝えます。
@model function failure_model(simulator)
# 初期状態のサンプリング
s ~ initial_state_distribution()
# 摂動のサンプリング
w = Vector{Float64}(undef, simulator.n_steps)
for i in 1:simulator.n_steps
w[i] ~ disturbance_distribution()
end
# 軌道の生成とスムーズ化された失敗密度の計算
tau = (s, w)
trajectory = rollout(simulator, tau)
robustness = evaluate_robustness(trajectory)
distance = max(0, robustness)
# スムーズ化された失敗密度を目標分布として使用
Turing.@addlogprob! smooth_failure_density(distance, epsilon)
endこのモデルを定義したら、Turingの推論機能を使用して失敗確率を推定できます：
# MCMCを使用してサンプリング
chain = sample(failure_model(simulator), NUTS(), 1000)
# 失敗確率の推定
failure_samples = filter(tau -> is_failure(rollout(simulator, tau)), get_samples(chain))
p_fail_estimate = length(failure_samples) / length(get_samples(chain))このように、確率的プログラミングは複雑な推論タスクを簡潔に表現し、自動的に実行するための強力なツールです。特に、失敗確率の推定のような複雑な問題に対して、実装の複雑さを大幅に減らしつつ、最先端の推論アルゴリズムの恩恵を受けることができます。
8.3 Julia言語とTuring.jlの使用例
Sydney Katz：確率的プログラミングの具体例として、Julia言語とTuring.jlライブラリの使用例を詳しく見ていきましょう。このアプローチは、失敗確率の推定のような複雑な推論タスクを簡潔に表現し、自動化するための強力な方法です。
教科書に掲載されているアルゴリズムを見ると、失敗分布からサンプリングするためのTuring.jlモデルの実装が示されています。このコードの構文をすべて理解する必要はありませんが、重要なのは、システムのダイナミクスを確率的プログラムとして簡単に指定できるという点です。
@model function failure_model(simulator)
# 初期状態のサンプリング
s ~ initial_state_distribution()
# 摂動のサンプリング
w = Vector{Float64}(undef, simulator.n_steps)
for i in 1:simulator.n_steps
w[i] ~ disturbance_distribution()
end
# 軌道の計算
tau = (s, w)
# スムーズ化された失敗密度を目標分布として使用
Turing.@addlogprob! smooth_failure_density(tau, simulator, epsilon)
endこのモデルでは、まず初期状態を初期状態分布からサンプリングし、続いて各ステップの摂動を摂動分布からサンプリングします。そして、これらのサンプリング値を使ってシステムの軌道を計算し、スムーズ化された失敗密度をモデルの対数確率に追加します。
このモデルを定義した後、Turingが提供する様々なサンプリングアルゴリズムを使って、このモデルからサンプリングできます：
# No-U-Turn Sampler (NUTS)を使用してサンプリング
chain = sample(failure_model(simulator), NUTS(), 1000)
# または他のサンプラーを使用することも可能
# chain = sample(failure_model(simulator), HMC(), 1000)
# chain = sample(failure_model(simulator), MH(), 1000)ここでは、No-U-Turn Sampler（NUTS）、ハミルトニアンモンテカルロ（HMC）、メトロポリス・ヘイスティングス（MH）など、様々なサンプラーを選択できます。これらは全て異なる特性を持ち、問題に応じて最適なものを選べます。
サンプリングが完了したら、得られたサンプルから失敗確率を推定できます：
# 失敗サンプルを抽出
failure_samples = filter(tau -> is_failure(rollout(simulator, tau)), get_samples(chain))
# 失敗確率の推定
p_fail_estimate = length(failure_samples) / length(get_samples(chain))この例からわかるように、確率的プログラミングを使うことで、複雑な失敗確率の推定問題を数行のコードで表現し、解くことができます。自分でサンプリングアルゴリズムを実装する必要はなく、Turingが提供する高度なアルゴリズムをそのまま利用できます。
Pythonユーザーのために、同様の機能を提供するPythonライブラリもあります。PyMC、Pyro、Stanなどが一般的です。Stanが最も古くからあるものですが、C++で書かれており、Pythonとの統合は少し難しい場合があります。
確率的プログラミングの大きな利点は、モデルと推論を明確に分離できることです。モデルを一度定義すれば、様々な推論アルゴリズムを簡単に試すことができ、それぞれの性能を比較できます。また、モデルの変更も容易であり、様々な仮定や条件の下での失敗確率を簡単に推定できます。
このアプローチは特に、複雑なシステムや高次元の状態空間を持つシステムの失敗確率推定において強力です。これらの場合、従来の手法では実装が複雑になりがちですが、確率的プログラミングを使うことで、実装の複雑さを大幅に減らしつつ、効率的な推論が可能になります。
