やたらと長いタイトルになったが表題の通り。メモがてらなので内容は期待しないように。「導出ラブ」の人は下記のページが勉強になるぞ。

拙者の場合…過程や…！方法なぞ…！どうでもよいのだァーッ！

というわけで、頑張って実装しませう。ちなみに、趣味的な問題ではなくて、真面目な問題ついて、正規分布の平均及び分散をデータから推定したいあなたはPyMCでもStanでもいいから、そういうの使おうね！時間の無駄だから。無駄無駄無駄無駄ァーーーッ！

とはいえ、全く何も分からないとコードをかけないので、ベイズ統計のおさらい。

ベイズ統計

ベイズ統計とは、ベイズという苗字の男が編み出した、ふわりとした統計の技術だ。主に下記の定理にその思想が集約されている。これにより、「観察と信念の更新」というアドホックな哲学により、神羅万象（神の意志）を知ることができるというありがたい定理なので、ちゃんと拝むように（発案者のトーマスベイズは牧師だったそうです…昔は宗教家も賢「かった」んだねぇ）。

ちなみに、理系のおねいさんやお兄さんたちは、小文字のpを使うことが多いぞ。ふんわりと、連続分布の時は小文字、みたいな業界ルールがある気がするけど、式の意味はおんなじ。

さて、このベイズの定理だけれども、原因Aと結果Bから、AによってBが生じる因果を推定するという、ブッチ神父も真っ青な深遠な定理である。AをパラメータΘ、BがデータD、パラメータに対するデータの尤度をP(D|Θ)とすることで、この定理を使うとアーラ不思議、（あるデータDのもとでの）パラメータΘの分布P(Θ|D)が分かっちゃうのだ。なんかすごそうだぜ。

パラメータ推定

というわけで、「割り算」と「掛け算」しかないので簡単そうに見えるが、実は問題はそう簡単ではないのである。普通、パラメータ推定というと、こういうのを想像する：

Θ' = optimize(model, Θ, data)

なんかよく分からないけど、データを入れるとパラメータΘが更新されて、より尤もらしい、新しいパラメータΘ'が手に入る、みたいなイメージである。ところが、ベジタリアン…じゃねーや、ベイジアンと呼ばれる人々は、このようには考えない。そのΘっていうのにも、不確実性があるよね？と言い出すのである。たとえば、正規分布の平均を推定する場合、頭がまとも（あるいはちょっと程度の低い）私たちは、μとσがわからんのだから、新しいμ'とσ'をくれよ、そしたらコード書けるから、コード！と思うのだが、「その場合、パラメータ推定に関する不確実性を無視していることになるよね？」などと屁理屈をこねてくるので、もう、どうしたらいいか分からない。

事前分布

というわけで、パラメータΘについて、何らかの分布を考えようと言い出すのがベイジアンと呼ばれる人々である。これは、ある程度の「お約束」があって「何とか分布の事前分布は、何とか」と暗記することになる。もちろん、「どうしゅつ」することもできるが、それがどんだけのハードワークか知りたい人は、冒頭でリンクしたページを見てみような。やりたくなくなるから。

今回の場合、求めたい分布は正規分布なので、パラメータといえば平均μと分散σ^2と相場が決まっているが、なにしろベイジアンなので、もはやこのμとか、σがただの数値だとは考えない。そもそも、分からないものなのに3とか4とか言うのはおかしいでしょ？というわけで、求めようとする正規分布の平均や分散も、別のなんかの分布に従う確率変数だ、と考えねばならぬ。

それが、事前分布（ここでエコーがかかる）。

ちなみに、正規分布の平均μの場合、実は平均μの分布も正規分布であることが知られている。母集団の平均を推定するとかってとき、中心極限定理から（手元のデータから推定された）平均μ'は正規分布に従う確率変数とみなして～などという理屈ものと、信頼区間を計算したことを思い出そう。よって、求めたい正規分布N(μ,σ)の「パラメータμ」は正規分布するはずなので、μの事前分布は正規分布だ！でも、求めたい正規分布のパラメータμが従う正規分布の平均と分散は、どうやって決めればよいのだろう？

実は、これがベイズ統計は主観が入り込むからなんとか、などとな癖をつけられる原因である。ぶっちゃけ、正規分布でさえあれば何でもよく、それは「あなた」のμに対する信念（あるいは、何らかの知識）を反映していると考える。もし、私のように頭が空っぽなら平均0、分散1でいいや、と思えばいいのである。ただし、データが十分でない場合、推定結果には事前分布で設定した分のバイアスがかかるので注意。

ちなみに、平均μは正規分布に従うとして、分散σはどんな事前分布にすればいいのだろうか？と疑問に思うところだが、私の綿密な調査の結果、それは逆ガンマ分布とよばれるありがたい分布であることが判明している。そして、逆ガンマ分布はscipy.statsにあるので、一石二鳥である。てか、正規分布の分散なので、本来逆chi2分布（カイ二乗検定とかに出てくるアイツ）に従うというのがほんとうなのだが、逆ガンマ分布はパラメータ次第で逆chi2分布になるので、逆ガンマ分布を使うのが良いらしい（そしてscipy.stats.invgammaという実装もすでにある）。

さて、こんがらがってきたので問題を整理する。

前提：ここに正規分布Nに従って分布するデータがあるが、平均μと分散σの両方が不明である。
要請：よって、データからNのパラメータμとσを求めたい。
要求：そのためには、まずμとσが従うであろう事前分布を決める必要がある。

この二段階ロケット的な推定が、ベイズ統計のデフォルトなやり方である。そして、μとσが従う分布の「パラメータ～を～、決めるのは、あな～た～」である。もし、パラメータμとσについて、既にあらかた見当がついているのなら、パラメータμが従う正規分布（求めたい正規分布Nとは別のやつね、念のため）の平均と分散は、貴方の当て推量を反映したものにすればよい。たとえば、Nの平均μが、なんとなく5くらいと思っているのなら、μは平均5くらいの正規分布に従うべきである。そして、その当て推量が、どの程度確からしいかを分散にすればよさそうだ。とまぁ、こういう小難しいことをする。

導出

さて、私が記事を書くときは、「可能な限り、数式は避けて通る」方式を貫こうと思っているのだが、一応どんな計算なのかくらいは知っておいた方がいいので、一番上で引用したページからちょっと大事そうな部分だけ抜粋。

まず、これな。これが超大事。細かいことは元のページを見てもらうとして、このπというのが、パラメータの事前分布を示す。今、我々は謎の正規分布Nの平均μと分散はσを推定するため、その事前分布を考えたのでした。これがその事前分布である！なぜ、右辺が掛け算になっているのかというと、私たちは二つのパラメータμとσを同時に推定しようとしているからなのだ（と思われる）。そして、平均が従う分布は分散に左右されそうなのでπ(μ|σ^2)という分布になり、分散は分散なのでπ(σ^2)という形にして分離すると…なるほどねー！そして、平均μと分散σの同時分布が知りたいので（同時分布は確率の掛け算）それを掛け算にすればオッケーということだそうです。では、このπの形を具体的に見てみよう。

平均μの事前分布：

分散σの事前分布：

へー。あー、上のやつはなんとなく正規分布の式に見えないこともないわ。2πがあるし、expって書いてるもん。ただ下のやつな。これが意味不明だけど、引用元のページによるとInv-Scaled-Chi-Squaredというらしい。なんか、新型のガンダムみたいな名前だね。ただ、こういうのは、ほら、scipyとかにあるから。ライブラリが…だから、分かってなくともよいと思われ！

さて、事前分布の式が分かったのはいいとして、これをどうすればいいのだろう？そう、これを「更新」しなければならない。幸いなことに、ベイジアンたちも、いや、「更新にも不確実性が～」と言い出すことはない（し、そこに不確実性は考えなくてよい、たぶん）。よって、例えばデータX=np.array([1,6,-5,4,-9,0])のようなものが手に入ったら、そこからこの分布たちのパラメータを更新することができる。求めたい正規分布Nが従うμ…が従う正規分布の平均とか、分散を更新するというわけである。いやー、めんどうだね！

というわけで、純真無垢な僕らは「どうしゅつ」についてはこのくらいにしよう。欲しいのは、更新の計算式である。

事前分布の更新

さて、具体的にどういう計算をするのかだけれども、それはありがたい外人が書いてくれているので、下記のページをご覧ください。

冒頭で紹介したページと違って、導出はないんだけど、こっちは事前分布が逆ガンマ分布になるとちゃんと書かれているので、実装するならこっちの記事を見た方が良いと思われ。さっきの分散の事前分布の式を見たとき、私のうすぼんやりの第六感が、これは「なんとか分布」によく似ていると囁いたので、さらに調査してみたところ、やっぱりなんとか分布（逆ガンマ分布）だとわかった。さすが私。調査力抜群。ちなみに、その噂の逆ガンマ分布というのは下記の通り。

さあ、先ほどのInverse-Scaled-Chi-Squaredと人力パターンマッチしてみよう。なんとなく、同じ式に見えるはずだ！よくわからなくとも、α=v0/2くらいは分かるよね。見たまんま。で、ベータはσ^2*v0/2っぽいね。うん。というわけで、そういうことである。そして、こいつらの更新式は次の通り：

このxバーみたいなやつはサンプル平均、S^2はサンプルから計算した分散にすればよいとのこと。いやー、これでやっと計算できますわ。

実装

ながかったぜ。さて、これで実装の準備が整ったのでコード書くなり。

from scipy import stats
import matplotlib.pyplot as plt
import numpy as np

def plot_norm(ax, mu, sigma, xmin=None, xmax=None):
    xmin = xmin if xmin else -3 * sigma + mu
    xmax = xmax if xmax else  3 * sigma + mu
    xs = np.linspace(xmin, xmax, 100)
    ys = stats.norm(mu, sigma).pdf(xs)
    ax.plot(xs, ys)

def plot_invgamma(ax, alpha, beta, **kwargs):    
    xmin, xmax = 0, 20
    xs = np.linspace(xmin, xmax, 250)
    ys = stats.invgamma(alpha, scale=beta).pdf(xs)
    #ys = inv_gamma(xs, alpha, beta)
    ax.plot(xs, ys, **kwargs)

class NormPrior:

    def __init__(self, m0, k0, v0, vs0, sigma0):
        self.sigma = sigma0
        self.m = m0
        self.k = k0
        self.v = v0
        self.vs = vs0

    # 偉い外人が教えてくれた小難しい更新式
    def update(self, x_hat, x_S, n):
        m0, k0, v0, vs0 = self.m, self.k, self.v, self.vs
        self.k = k0 + n
        self.m = (k0 * m0 + n * x_hat) / self.k
        self.v = v0 + n
        self.vs = vs0 + (n-1)*x_S + k0*n*((x_hat-m0)**2) / self.k

    # これ、標準偏差なのでルートを取るのが正しい（はず）
    def mean(self):
        return stats.norm(self.m, np.sqrt(self.sigma/self.k))

    # さっき見比べたとき、α=v0/2で、β=v0*σ^2/2だったのでした
    def var(self):
        return stats.invgamma(self.v/2, self.vs/2)

    def plot(self):
        f, ax = plt.subplots(1,2)
        plot_norm(ax[0], self.m, np.sqrt(self.sigma / self.k))
        plot_invgamma(ax[1], self.v/2, self.vs/2)

このNormPriorというのが、これまでのありがたい数学的真理を、numpyのコードに落とし込んだものである！いやー、わかりやすい、見やすい。実際コードにすると、量自体はこんなものなのだが、調べたりなんだりで小一時間かかるのが、この手の分野である。どっかの小うるさい人々が、なにかにつけて「データを出して議論しよう」などというのを耳にする昨今だけど、我々（三流）ぷろぐらまーも「コードを出して」議論するのがいいね。間違ってたらすぐにわかるから。

ちなみに、statsのinvgammaだけど、最初のパラメータがα、二番目がβで良い模様です。なんかねー、ドキュメントがあいまい。αについてはちゃんと書かれてるんだけど、βがなかったことにされてるもの。そういうとこがねー、scipy.statsのドキュメントはいただけないね（ご意見番調で喋ってるつもり）。

では、動作確認してみるなり。

# 真の分布は平均7で分散4（標準偏差2）
n = stats.norm(7.0, 2.0)
# 事前分布のパラメータは「てきとう」
# 本来、手元にあるデータから標本平均何かを求めて設定してもよく、
# そういうのを経験ベイズとか言うらしい（あいまい）
prior = NormPrior(0.0, 1.0, 2.0, 2.0, 1.0)
prior.plot()
plt.show()

samp = n.rvs(10)
prior.update(np.mean(samp), np.var(samp), 10)
prior.plot()
plt.show()

samp = n.rvs(10)
prior.update(np.mean(samp), np.var(samp), 10)
prior.plot()
plt.show()

samp = n.rvs(100)
prior.update(np.mean(samp), np.var(samp), 100)
prior.plot()
plt.show()

一応、こういうのは間違ってても気づかないパターンが結構あるので、ちゃんとあっているかどうか図示すると分かりやすいぞ。今回、全くもって適当に設定した事前分布を、平均7.0、分散4.0（標準偏差2.0）の正規分布から発生させたデータで更新してみた。もし、私の計算が正しければ、時速140キロを超えたところで、何かが起きるはずだ…じゃねーや、事前分布のピークがそれぞれ7.0とか4.0に近づくはずである。

あー、適当な感じになってるね。これが、私が独断と偏見により設定した事前分布。平均μの分布は、N(0,1)みたいな感じになっていますね。分散は、0.6とかそのくらいになってます（実際のデータの分散は4.0だから大外れですね）。では、データを使って事前分布を更新するなり。

おー、良い感じに近づいているぞ。

左の図は、スケールが変わっているのでわかりずらいが、X軸を見るとμが従う分布の幅が狭まっているのが分かる。ようするに、沢山（といっても30サンプルくらい）データを食わして、分布を更新したので、平均μの存在範囲がどんどん狭まってきているということのようです。先ほどの更新式を見てもらうと分かる通り、食わしたデータの数nが大きければ、それだけμの事前分布の分散が小さくなるようになっています。分散の事前分布（逆ガンマ分布）についても、存在範囲が狭まり真値(4.0)にちゃんと近づいているみたいです。

まー、多分これであってるんじゃないかな？

まとめ

正規分布の平均分散の両方を、ベイズ更新する場合は、このような手順を踏むようです。おさらいすると：

まず、パラメータの洗い出しをする。今回の場合、未知の正規分布の平均μと分散σ^2。
それぞれのパラメータについて、事前分布を求める（ググる）。
事前分布をデータから更新するための更新式を求める（ググる）。
実装して、グラフを描いて、あってそうか確かめる。

以上です！