基本的な定義を邪魔にならないようにしましょう。マルコフ連鎖モンテカルロ（MCMC）法では、分布を計算できなくても、分布からサンプルを計算できます。

これは何を意味するのでしょうか？バックアップして、モンテカルロサンプリングについて話しましょう。

モンテカルロサンプリング

モンテカルロ法とは何ですか？

「モンテカルロ法、またはモンテカルロ実験は、数値結果を得るために繰り返しランダムサンプリングに依存する幅広いクラスの計算アルゴリズムです。」 （（ ウィキペディア ）

それを分解しましょう。

例1：面積の推定

以下に示すような不規則な形をしていると想像してください。

次のうち、IoTに関する懸念事項ではないものはどれですか？

そして、あなたはこの形で囲まれた領域を決定する任務を負っています。使用できる方法の1つは、形状に小さな正方形を作成し、正方形を数えることです。これにより、面積をかなり正確に概算できます。しかし、それは難しく、時間がかかります。

モンテカルロサンプリングが救助に！

まず、形状の周囲に既知の領域の大きな正方形、たとえば50cm2を描画します。次に、この正方形を「吊るして」、その形状にランダムにダーツを投げ始めます。

次に、長方形のダーツの総数と、関心のある形状のダーツの数を数えます。使用された「ダーツ」の総数が100で、そのうちの22が形状内に収まったと仮定します。これで、面積は次の簡単な式で計算できます。

形状の面積=正方形の面積*（形状のダーツの数）/（正方形のダーツの数）

したがって、私たちの場合、これは次のようになります。

形状の面積= 50 * 22/100 = 11 cm2

「ダーツ」の数に10を掛けると、この近似は実際の答えに非常に近くなります。

形状の面積= 50 * 280/1000 = 14 cm2

これは、モンテカルロサンプリングを使用して、上記のような複雑なタスクを分解する方法です。

大数の法則

投げるダーツが多いほど、面積の近似値が近くなりました。これは、 大数の法則 ：

「大数の法則は、同じ実験を何度も実行した結果を表す定理です。法律によれば、多数の試行から得られた結果の平均は期待値に近いはずであり、より多くの試行が実行されるにつれて近づく傾向があります。」

これにより、次の例である有名なモンティホール問題に移ります。

例2：モンティホール問題

ザ・ モンティホール問題 非常に有名な頭の体操です：

「ドアは3つあり、1つは後ろに車があり、他のドアは後ろに山羊がいます。ドアを選択すると、ホストが別のドアを開き、その後ろに山羊がいることを示します。それから彼はあなたがあなたの決定を変えたいかどうかあなたに尋ねます。あなたは？どうして？何故なの？'

最初に頭に浮かぶのは、切り替えても切り替えなくても勝つ可能性は同じだということですが、そうではありません。同じことを示す簡単なフローチャートを作成しましょう。

車がドアの後ろにあると仮定します3：

したがって、切り替えた場合は2/3回勝ち、切り替えなかった場合は1/3回しか勝ちません。

サンプリングを使用してこれを解決しましょう。

wins = [] for i in range(int(10e6)): car_door = assign_car_door() choice = random.randint(0, 2) opened_door = assign_door_to_open(car_door) did_he_win = win_or_lose(choice, car_door, opened_door, switch = False) wins.append(did_he_win) print(sum(wins)/len(wins))

assign_car_door()関数は、ドア0、1、または2を選択する単なる乱数ジェネレーターであり、その後ろに車があります。 assign_door_to_openを使用する後ろに山羊がいて、選択したドアではないドアを選択すると、ホストがドアを開きます。 win_or_lose戻り値 true または false 、あなたが車に勝ったかどうかを示すには、ドアを切り替えたかどうかを示すブール「スイッチ」が必要です。

このシミュレーションを1,000万回実行してみましょう。

クレジットカードの暗証番号を解読する方法

• 切り替えない場合の勝率：0.334134
• 切り替えた場合の勝率：0.667255

これは、フローチャートが示した答えに非常に近いものです。

実際、このシミュレーションを実行すればするほど、答えは真の値に近づき、大数の法則が検証されます。

この表からも同じことがわかります。

シミュレーションの実行	切り替えた場合の勝率	切り替えない場合の勝率
10	0.6 0.6	0.2
10 ^ 2	0.63	0.33
10 ^ 3	0.661	0.333
10 ^ 4	0.6683	0.3236
10 ^ 5	0.66762	0.33171
10 ^ 6	0.667255	0.334134
10 ^ 7	0.6668756	0.3332821

ベイジアンの考え方

「統計のより古典的なバージョンとして知られている頻度主義者は、確率がイベントの長期的な頻度であると想定しています（したがって、授与されたタイトル）。」

「ベイズ統計は、確率のベイズ解釈に基づく統計の分野における理論であり、確率はイベントに対するある程度の信念を表します。信念の程度は、以前の実験の結果など、イベントに関する事前の知識、またはイベントに関する個人的な信念に基づいている可能性があります。」 -から ハッカーのための確率的プログラミングとベイジアン手法

これは何を意味するのでしょうか？

頻度主義的な考え方では、長期的には確率を調べます。頻度主義者が自動車事故が発生する可能性が0.001％であると言う場合、それは、無限の自動車旅行を考慮すると、それらの0.001％が交通事故で終わることを意味します。

私たちが以前の信念から始めるので、ベイジアンの考え方は異なります。私たちが0の信念について話す場合、それはあなたの信念がイベントが決して起こらないということを意味します。逆に、1の信念は、それが起こると確信していることを意味します。

次に、データの観察を開始したら、データを考慮に入れるようにこの信念を更新します。これをどのように行うのですか？を使用して ベイズの定理 。

ベイズの定理

....（1）

分解してみましょう。

• P(A | B)イベントBが与えられた場合のイベントAの確率を示します。これは事後確率であり、Bは観測されたデータです。したがって、観測されたデータを考慮して、イベントが発生する確率は基本的に何であるかを示しています。
• P(A)は事前であり、イベントAが発生するという私たちの信念です。
• P(B | A)は尤度であり、Aが真であると仮定してデータを観測する確率はどのくらいですか。

例として、がん検診を見てみましょう。

がんの確率

患者がマンモグラムを完成させに行き、マンモグラムが陽性に戻ったとしましょう。患者が実際に癌にかかっている確率はどれくらいですか？

確率を定義しましょう：

• 女性の1％が乳がんを患っています。
• マンモグラムは、実際に存在する時間の80％で癌を検出します。
• マンモグラムの9.6％は、実際には癌がないのに癌があると誤って報告しています。

したがって、マンモグラムが陽性に戻った場合、80％の確率で癌にかかっていると言えば、それは間違いです。がんになることはまれな出来事である、つまり、女性の1％だけが乳がんにかかっていることを考慮に入れません。これを事前確率としてとる必要があります。ここでベイズの定理が役立ちます。

P（C + | T +）=（P（T + | C +）* P（C +））/ P（T +）

• P(C+ | T+)：これは、検査が陽性であったことを考えると、癌が存在する確率です。これが私たちが関心を持っていることです。
• P(T+ | C+)：これは、癌がある場合にテストが陽性である確率です。これは、上記のように、80％= 0.8に等しくなります。
• P(C+)：これは事前確率であり、個人が癌を患う確率であり、1％= 0.01に相当します。
• P(T+)：これは、何があってもテストが陽性である確率であるため、2つの要素があります。 P（T +）= P（T + | C-） P（C-）+ P（T + | C +） P（C +）
• P(T+ | C-)：これは、検査が陽性に戻ったが癌がない確率です。これは0.096で与えられます。
• P(C-)：これは癌を患う確率が1％であるため、癌を患っていない確率です。これは99％= 0.99に等しくなります。
• P(T+ | C+)：これは、あなたが癌を持っているとすると、検査が陽性に戻った確率であり、これは80％= 0.8に等しいです。
• P(C+)：これは癌になる確率であり、1％= 0.01に相当します。

これらすべてを元の数式に代入します。

したがって、マンモグラムが陽性に戻ったことを考えると、患者が癌を患う可能性は7.76％です。最初は奇妙に思えるかもしれませんが、それは理にかなっています。この検定では、9.6％の確率で誤検出が発生するため（かなり高い）、特定の母集団には多くの誤検出が発生します。まれな病気の場合、陽性の検査結果のほとんどは間違っています。

ここで、モンティホール問題を再検討し、ベイズの定理を使用して問題を解決しましょう。

モンティホール問題

事前定義は次のように定義できます。

• 最初にドアAを選択するとします。
• P（H）=⅓、⅓。 3つのドアすべての1/3。つまり、ドアが開いて山羊が現れる前に、車がいずれかのドアの後ろにある可能性が同じでした。

可能性は次のように定義できます。

• P(D|H)、イベントDは、モンティがドアBを選択し、その後ろに車がないことです。
• P(D|H) =車がドアAの後ろにある場合は½-ドアBを選択する確率は50％、ドアCを選択する確率は50％であるため。
• P(D|H)車がドアBの後ろにある場合、= 0です。これは、車がドアBの後ろにある場合、ドアBを選択する可能性が0％であるためです。
• P(D|H) = 1車がCの後ろにあり、Aを選択した場合、彼がドアBを選択する確率は100％です。

P(H|D)に関心があります。これは、他のドアの1つの後ろにヤギが表示されている場合に、車がドアの後ろにある確率です。

s corp ccorpとllcの違い

後部P(H|D)から、AからドアCに切り替えると、勝つ確率が高くなることがわかります。

メトロポリス-ヘイスティングス

それでは、これまでに取り上げたすべてを組み合わせて、その方法を理解してみましょう。 メトロポリス-ヘイスティングス 動作します。

Metropolis-Hastingsは、ベイズの定理を使用して、直接サンプリングが困難な複雑な分布の事後分布を取得します。

どうやって？基本的に、空間からさまざまなサンプルをランダムに選択し、式（1）の確率の比率P（B）を調べているため、新しいサンプルが最後のサンプルよりも後方からのものである可能性が高いかどうかを確認します。キャンセル：

P（acceptance）=（P（newSample）*新しいサンプルの可能性）/（P（oldSample）*古いサンプルの可能性）

新しい各サンプルの「尤度」は、関数によって決定されます f 。それが理由です f サンプリングしたい後部に比例する必要があります。

θ 'を受け入れるか拒否するかを決定するには、新しく提案されたθ'ごとに次の比率を計算する必要があります。

ここで、θは古いサンプルです。P(D| θ)はサンプルθの尤度です。

これをよりよく理解するために例を使用しましょう。データはあるが、それに適合する正規分布を見つけたいとしましょう。

X〜N（平均、標準）

次に、平均と標準の両方の事前確率を定義する必要があります。簡単にするために、両方とも平均1と標準2の正規分布に従うと仮定します。

平均〜N（1、2）
std〜N（1、2）

ここで、いくつかのデータを生成し、真の平均と標準がそれぞれ0.5と1.2であると仮定しましょう。

import numpy as np import matplotlib.pyplot as plt import scipy.stats as st meanX = 1.5 stdX = 1.2 X = np.random.normal(meanX, stdX, size = 1000) _ = plt.hist(X, bins = 50)

PyMC3

まず、というライブラリを使用しましょう PyMC3 データをモデル化し、平均と標準の事後分布を見つけます。

import pymc3 as pm with pm.Model() as model: mean = pm.Normal('mean', mu = 1, sigma = 2) std = pm.Normal('std', mu = 1, sigma = 2) x = pm.Normal('X', mu = mean, sigma = std, observed = X) step = pm.Metropolis() trace = pm.sample(5000, step = step)

コードを見てみましょう。

最初に、平均と標準の事前分布、つまり平均1と標準2の正規分布を定義します。

x = pm.Normal（ 'X'、mu =平均、シグマ= std、観測値= X）

この行では、関心のある変数を定義します。以前に定義した平均と標準を取り、以前に計算した観測値も取ります。

結果を見てみましょう：

ロボットを構築してプログラムする

_ = plt.hist(trace['std'], bins = 50, label = 'std') _ = plt.hist(trace['mean'], bins = 50, label = 'mean') _ = plt.legend()

stdは1.2を中心とし、平均は1.55です-かなり近いです！

メトロポリス・ヘイスティングスを理解するために、最初からやってみましょう。

ゼロから

提案の生成

まず、メトロポリス・ヘイスティングスが何をしているのかを理解しましょう。この記事の前半で、私たちは言いました、 「メトロポリス・ヘイスティングスは、空間からさまざまなサンプルをランダムに選択します。」 しかし、どのサンプルを選択するかをどうやって知るのでしょうか？

これは、提案配布を使用して行います。これは、現在受け入れられているサンプルを中心とした正規分布であり、STDは0.5です。 0.5がハイパーパラメータである場合、微調整できます。提案のSTDが多いほど、現在受け入れられているサンプルからさらに検索されます。これをコーディングしましょう。

分布の平均と標準を見つける必要があるため、この関数は現在受け入れられている平均と現在受け入れられている標準を取得し、両方の提案を返します。

def get_proposal(mean_current, std_current, proposal_width = 0.5): return np.random.normal(mean_current, proposal_width), np.random.normal(std_current, proposal_width)

提案の承認/拒否

それでは、提案を受け入れるか拒否するロジックをコーディングしましょう。それは2つの部分に分かれています： 前 そしてその 可能性 。

まず、提案が前から来る確率を計算しましょう。

def prior(mean, std, prior_mean, prior_std): return st.norm(prior_mean[0], prior_mean[1]).pdf(mean)* st.norm(prior_std[0], prior_std[1]).pdf(std)

それはただかかります 平均 そして 時間 これがどの程度可能性があるかを計算します 平均 そして 時間 から来た 事前配布 私たちが定義したこと。

尤度を計算する際に、私たちが見たデータが提案された分布からのものである可能性を計算します。

def likelihood(mean, std, data): return np.prod(st.norm(mean, std).pdf(X))

したがって、各データポイントについて、提案された分布に由来するそのデータポイントの確率を乗算します。

ここで、現在の平均と標準、および提案された関数に対してこれらの関数を呼び出す必要があります。 平均 そして 時間 。

prior_current = prior(mean_current, std_current, prior_mean, prior_std) likelihood_current = likelihood(mean_current, std_current, data) prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std) likelihood_proposal = likelihood(mean_proposal, std_proposal, data)

コード全体：

def accept_proposal(mean_proposal, std_proposal, mean_current, std_current, prior_mean, prior_std, data): def prior(mean, std, prior_mean, prior_std): return st.norm(prior_mean[0], prior_mean[1]).pdf(mean)* st.norm(prior_std[0], prior_std[1]).pdf(std) def likelihood(mean, std, data): return np.prod(st.norm(mean, std).pdf(X)) prior_current = prior(mean_current, std_current, prior_mean, prior_std) likelihood_current = likelihood(mean_current, std_current, data) prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std) likelihood_proposal = likelihood(mean_proposal, std_proposal, data) return (prior_proposal * likelihood_proposal) / (prior_current * likelihood_current)

それでは、すべてを結び付けて必要な事後サンプルを生成する最終関数を作成しましょう。

後部の生成

ここで、上記で記述した関数を呼び出して、後部を生成します。

def get_trace(data, samples = 5000): mean_prior = 1 std_prior = 2 mean_current = mean_prior std_current = std_prior trace = { 'mean': [mean_current], 'std': [std_current] } for i in tqdm(range(samples)): mean_proposal, std_proposal = get_proposal(mean_current, std_current) acceptance_prob = accept_proposal(mean_proposal, std_proposal, mean_current, std_current, [mean_prior, std_prior], [mean_prior, std_prior], data) if np.random.rand()

改善の余地

ログはあなたの友達です！

思い出してくださいa * b = antilog(log(a) + log(b))およびa / b = antilog(log(a) - log(b)).

これは、非常に小さい確率を処理し、乗算するとゼロになるため、対数確率を追加して問題を解決するため、私たちにとって有益です。

したがって、以前のコードは次のようになります。

def accept_proposal(mean_proposal, std_proposal, mean_current, std_current, prior_mean, prior_std, data): def prior(mean, std, prior_mean, prior_std): return st.norm(prior_mean[0], prior_mean[1]).logpdf(mean)+ st.norm(prior_std[0], prior_std[1]).logpdf(std) def likelihood(mean, std, data): return np.sum(st.norm(mean, std).logpdf(X)) prior_current = prior(mean_current, std_current, prior_mean, prior_std) likelihood_current = likelihood(mean_current, std_current, data) prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std) likelihood_proposal = likelihood(mean_proposal, std_proposal, data) return (prior_proposal + likelihood_proposal) - (prior_current + likelihood_current)

受け入れの確率のログを返すので：

if np.random.rand()

になる：

if math.log(np.random.rand())

新しいコードを実行して、結果をプロットしてみましょう。

cvv2017で働くクレジットカード番号

_ = plt.hist(trace['std'], bins = 50, label = 'std') _ = plt.hist(trace['mean'], bins = 50, label = 'mean') _ = plt.legend()

ご覧のとおり、 時間 は1.2を中心とし、 平均 1.5で。やりました！

前進する

これまでに、メトロポリス・ヘイスティングスがどのように機能するかを理解し、どこで使用できるのか疑問に思われるかもしれません。

上手、 ハッカーのためのベイズ法 は確率的プログラミングを説明する優れた本です。ベイズの定理とその応用についてもっと知りたい場合は、 ベイズだと思う アレン・B・ダウニーの素晴らしい本です。

読んでいただきありがとうございます。この記事がベイズ統計の驚くべき世界を発見するきっかけになることを願っています。

基本を理解する

MCMCは何の略ですか？

MCMCは、確率分布からサンプリングするためのメソッドのクラスであるマルコフ連鎖モンテカルロの略です。

意思決定に対するベイズのアプローチは何ですか？

意思決定へのベイジアンアプローチでは、最初に事前確率から始めます。これがあなたの信念です。次に、データが入ってくると、そのデータを組み込んでこれらの事前確率を更新し、事後確率を取得します。

ベイズモデルとは何ですか？

ベイズモデルは、確率を使用してモデル内のすべての不確実性を表す統計モデルです。