【自主セミナー】ベイズ統計オフラインセミナー第4回

開催日付: 2018/1/8

今回やった内容の範囲: ThinkBayes 第3章「推論」

今回は初めて「オフライン」でやりました！場所はyahooロッジ。

今回の資料はこちらです。(まぁ、相変わらずのクオリティですが…)

ベイズ統計セミナーVol4 from まえすとろ

www.slideshare.net

前回やってから1ヶ月半ほど経っているのでベイズ的に問題を解くプロセスをおさらいしてからサイコロ問題を解いていき、本題の「事後確率から推定値を計算する」というところに入りました。

この辺り、パラメータ推定とかでも結構ネックになってくる部分になってくるじゃないかと思います。「尤度」に推定したいパラメータが入っていたりとかするんです。これがよくベイズ統計でググって出てくる「パラメータが確率分布」というところに相当します。(実際はパラメータベースで立てた仮説の確率分布が、そうなっているんですが！)

あと、ここで大事なのは「事前確率分布の設定は適切に！」というところでしょうか。

信用区間とかもちょっとやったけど、やっぱり理解しきれていないなー。。。と思いました。もっとこの辺りを勉強して、信用区間を出してどう意思決定につなげるかを深く研鑽する必要があると痛感。

ちなみに、今回のセミナーで取り扱った機関車問題、テキストや資料となりそうなサイトがあったのでいかにURLを貼り付けておきます。

store.doverpublications.com

Fifty Challenging Problems in Probability

うーん、この辺りになるとちょっとずつ骨がある内容になって噛み砕くのに時間がかかりますね。。

ちなみに、べき乗則が $(\frac{1}{x})^\alpha \ (x=1,2,…)$ という形をしていたので、番外編として「リーマンゼータ関数」を紹介して、 $1 \leq \alpha$ における収束性について紹介しました。(一応結論だけ言うと、 $\alpha = 1$ のときは発散して、 $1 \le \alpha$ では収束します。)

そして、いかに私が今回の資料用のグラフを描画するに使ったコードを掲載します。(だいぶ初心者なコードですが) 以下のコードをライブラリとして実装して見ました。

github.com

#coding: utf-8

from code.thinkbayes import Suite
from code import thinkplot

class TrainCase1(Suite):

    def Likelihood(self, data, hypo):
        if hypo < data:
            return 0
        else:
            return 1.0 / hypo

class TrainCase2(Suite):

    def __init__(self, hypos, alpha=1.0):
        Suite.__init__(self)
        for hypo in hypos:
            self.Set(hypo, hypo**(-alpha))
        self.Normalize()

    def Likelihood(self, data, hypo):
        if hypo < data:
            return 0
        else:
            return 1.0 / hypo


# 適当に上限を決めて、一様分布でベイズ更新を行う。
def case1():
    thinkplot.PrePlot(1)

    for n0 in [500, 1000, 2000]:
        hypos = xrange(1, n0)
        suite = TrainCase1(hypos)
        suite.Update(60)
        suite.name = "N0 = " + str(n0)
        print "N0 = " + str(n0) + " mean."
        print suite.Mean()
        thinkplot.Pmf(suite)

    thinkplot.Save(root='train1',
                   xlabel='Number of trains',
                   ylabel='Probability',
                   formats=['pdf'])


# べき乗則に従った形でベイズ更新を行う。
def case2():
    thinkplot.PrePlot(1)

    for n0 in [500, 1000, 2000]:
        hypos = xrange(1, n0)
        suite = TrainCase2(hypos, 0.999)
        suite.Update(60)
        suite.name = "H = " + str(n0)
        print "H = " + str(n0) + " mean."
        print suite.Mean()
        thinkplot.Pmf(suite)

    thinkplot.Save(root='train2',
                   xlabel='Number of trains',
                   ylabel='Probability',
                   formats=['pdf'])

if __name__ == '__main__':
    #case1()
    case2()

とりあえず推論と信用区間について自分的にも少し理解を深めた会でした。

次回はThinkBayes 第4章。2月下旬ぐらい開催を予定。