乱数は LCG（固定 seed=42）で決定論的。再サンプリングで seed を更新します。
青実線=ターゲット p(x)／緑破線=包絡 M·q(x)／緑丸=受理点／赤×=棄却点
青棒=受理点のヒストグラム（正規化）／白線=ターゲット p(x)
棄却サンプリングは、提案 q(x) と定数 M ≥ sup p(x)/q(x) を選び、X ~ q(x), U ~ U(0,1) として受理判定 U < p(X)/(M·q(X)) を行います。
受理確率（任意の試行が受理される確率）：
$$P(\text{accept}) = \int \frac{p(x)}{M\,q(x)}\,q(x)\,dx = \frac{1}{M}\int p(x)\,dx = \frac{1}{M}$$
受理された標本の条件付き分布：
$$f_{X\,|\,\text{accept}}(x) = p(x)$$
期待される受理サンプル数：
$$\mathbb{E}[N_{\text{accept}}] = \frac{N}{M}$$
M が大きいほど包絡は厚くなり安全ですが、受理率は 1/M に比例して下がり計算量が増えます。
棄却サンプリングシミュレーターとは
🙋
確率分布 p(x) からサンプルを取るって、コンピュータは普通に乱数を生成すればいいだけじゃないんですか？なんで「棄却」とか必要なんですか？
🎓
いい質問だ。一様乱数や正規乱数みたいに「式が逆関数で解ける分布」なら直接サンプリングできるんだけど、ベイズ統計の事後分布なんかは正規化定数が分からなかったり、形が複雑だったりして直接生成できないことが多いんだ。そんなとき棄却サンプリングは「サンプリングしやすい提案 q(x) をかぶせて、当たりだけ拾う」というアイデアでこれを回避する。上のシミュレーターのデフォルトでは Beta(2,5) を一様分布で覆っているよ。
🙋
緑の破線が「M·q(x)」で、青い実線がターゲットなんですね。なんで M を掛けるんですか？
🎓
それは「M·q(x) ≥ p(x) を全ての x で満たす」必要があるからだ。つまり提案を上に押し上げて、青い山を完全に覆う「天井」を作るんだ。Beta(2,5) の最大値は x=0.2 で約2.458 だから、q=1（一様）なら M≥2.458 必要。デフォルトは安全側で M=3.0 にしてある。M を 2 に下げてみると、青線が緑破線を突き抜けるはずだよ。そうなると下の正しい分布から外れてしまう。
🙋
受理率が 33% って表示されてますけど、これって「1/M」と一致してますね。偶然じゃないんですか？
🎓
偶然じゃない、これが棄却サンプリングの核心定理だ。受理確率を積分すると ∫p(x)/(M·q(x))·q(x)dx = (1/M)∫p(x)dx = 1/M になる。だから 2000 試行で受理は約 666、率は約 33%。M を 5 にすれば 20%、10 にすれば 10% と一直線に効率が落ちる。シミュレーターの「理論最大効率 1/M」カードと「受理率」カードを見比べてごらん、ほぼ一致するはずだ。
🙋
「ターゲット分布」を 1（2峰）に切り替えてみました。受理率がガクッと下がりますね…！
🎓
そう、2峰や複雑な形状になると、提案ガウスでカバーするには M を大きく取らざるを得ない。提案幅 σ を小さくしすぎると裾でカバー不足、大きくしすぎると無駄が増える。実務ではこれが大問題で、高次元では受理率が指数的にゼロに近づく「次元の呪い」が発生する。そこで MCMC（マルコフ連鎖モンテカルロ）やハミルトニアン MC のような、もっと賢い手法が必要になるんだ。棄却サンプリングはその出発点であり、限界も同時に教えてくれる教科書的な手法なんだよ。
よくある質問
M を最小にしたいのですが、どう求めればよいですか？
理論的には M = sup_x p(x)/q(x) が最適で、これを採れば受理率が最大化されます。Beta(2,5) と一様提案では p の最頻値 x=0.2 で比が 2.458 となるため M=2.458 が最適です。実装では p と q を解析的に微分して上限を求めるか、グリッド評価で近似します。形状が分からない場合は適応的に M を更新する適応棄却サンプリング (ARS) を使うのが安全です。
ターゲットの正規化定数が分からなくても使えますか？
使えます。p(x) = p̃(x)/Z で Z が未知でも、M を sup p̃(x)/q(x) として比較すれば受理確率は p̃(x)/(M·q(x)) で計算でき、得られる標本は依然 p(x) に厳密に従います。ベイズ事後分布のように尤度×事前で Z が不明な場合に特に有用で、これがメトロポリス・ヘイスティングス法の前身となる重要な性質です。
高次元でなぜ効率が悪くなるのですか？
高次元では p(x) の質量が極めて狭い領域に集中し、提案 q(x) を p に完全にフィットさせるのが極端に難しくなるためです。比 p/q の上限が次元 d に対して指数的に増加し、受理率 1/M も指数的に減少します。例えば 10 次元のガウスターゲットを別のガウスで覆おうとすると、平均が少しずれただけで M が天文学的な値になります。高次元では MCMC、変分推論、HMC など別系統の手法が必須です。
なぜ受理サンプルだけが正しく p(x) に従うのですか？
受理イベントの条件付き分布を計算すると証明できます。X ~ q かつ U < p(X)/(M·q(X)) のもとで X の条件付き密度は f(x|accept) = q(x)·[p(x)/(M·q(x))] / (1/M) = p(x) となり、ちょうどターゲットになります。M·q(x) ≥ p(x) という条件が「天井で覆ったランダム点が一様にばらまかれる」幾何的解釈を保証し、これがその一様点を縦切りで見たものが p(x) になる根拠です。
実世界での応用
ベイズ統計と確率モデリング： 事後分布 p(θ|D) ∝ p(D|θ)p(θ) は通常正規化定数が解析的に求まらないため、棄却サンプリングを使ってパラメータ θ のサンプルを生成します。低〜中次元のモデルや、より高度な MCMC のウォームアップとして広く利用され、ベイズ推論の基礎構成要素となっています。
計算物理・分子シミュレーション： マクロな物理量を確率分布からのサンプリングで評価する重み付きモンテカルロ法では、ボルツマン分布 p(x) ∝ exp(-E(x)/kT) からの標本生成に棄却サンプリングが使われます。エネルギーランドスケープが単純な場合に有効で、メトロポリス法の理論的背景にもなっています。
CAE と信頼性解析： 製品の信頼性評価では、入力パラメータの分布（材料強度・寸法公差）から標本を生成し、有限要素法で応答を求める確率論的設計が行われます。入力分布が標準形でない場合に棄却サンプリングで標本生成し、得られた応答分布から破壊確率を推定します。
機械学習と生成モデル： 確率モデルの正規化定数評価、エネルギーベースモデルからの標本生成、レアイベント推定などで利用されます。GAN や拡散モデルの理論的解析にも棄却サンプリングの概念が登場し、現代の生成モデル研究の理解にも不可欠です。
よくある誤解と注意点
最も多い誤解は、「受理率が低くてもサンプルさえ得られれば品質は同じ」と考えてしまう ことです。確かに受理された標本は理論上厳密に p(x) に従いますが、受理率 1/M が極端に低いと、希望のサンプル数を得るのに膨大な計算コストがかかります。Beta(2,5) の例で M=10 を選ぶと受理率は 10%、つまり同じサンプル数を得るのに M=2.458 と比べ約 4 倍の試行が必要です。M はギリギリの値を選ぶのが基本で、シミュレーターで M スライダーを動かして「受理率」カードを見ながら最適点を体感してみてください。
次に多いのが、提案分布 q(x) の選び方を軽視する ミスです。q が p と形状的に大きく異なると、たとえ M·q(x) ≥ p(x) を満たしても、提案の大部分が p の山から外れた「無駄な領域」に落ち、棄却ばかりになります。シミュレーターで「ターゲット分布」を 1（2峰）にして「提案ガウス幅」を 0.5 にしてみてください。中央付近に提案が集中し、両峰の標本が極端に少なくなる様子が観察できます。良い q は「形状が p に近く、裾も p より十分広い」ものです。
最後に、条件 M·q(x) ≥ p(x) が破られた場合の影響 を理解しておくこと。これは全 x で成り立たねばならない必須条件で、1 点でも破られると受理確率が 1 を超え、得られる標本は p(x) から系統的に偏ります。シミュレーターで M=2 にしてターゲットを Beta(2,5) のままにすると、ピーク付近で青線が緑破線を突き抜け、ヒストグラムがピーク付近で過小評価されるのが見えます。実装ではグリッド評価や数値最適化で sup p/q を保守的に評価し、必ず安全係数を掛けて M を設定します。