Think Bayes ベイズの定理のメモ 3

続き

■サイコロ問題
4面のサイコロ、6面のサイコロ、8面のサイコロ、12面のサイコロ、20面のサイコロの入った箱を持っているとする。
箱からサイコロをランダムに選んで振って6が出たとする。各サイコロを降った確率はどうなるか。

from thinkbayes import Suite
class Dice(Suite):
 def Likelihood(self, data, hypo):
     if hypo < data:
        return 0
     else : 
        return 1.0/hypo

def main():
 suite = Dice([4,6,8,12,20])
 suite.Update(6)
 suite.Print()

if __name__ == '__main__':
    main()	

答え

4 0.0
6 0.392156862745
8 0.294117647059
12 0.196078431373
20 0.117647058824

■機関車問題
「ある鉄道会社では、機関車に1...Nという番号をつけている。ある日、60番という機関車を目撃したとすると、
鉄道会社は何台の機関車を所有しているのか推定せよ。」

これはサイコロの問題と同じになります。

したは60番を見た後に30番、90番、130番・・・・・と続いたとき

import thinkplot
from thinkbayes import Suite
class Train(Suite):
 def Likelihood(self, data, hypo):
     if hypo < data:
        return 0
     else : 
        return 1.0/hypo

def main():
  hypos = xrange(1, 1001)
  suite = Train(hypos)
  for data in [60,30,90,130,100,20,30,50,10,40,50,100,20,10,30,40,80,90,50,100]:
   suite.Update(data)
  print suite.Mean()
  thinkplot.PrePlot(1)
  thinkplot.Pmf(suite)

### ファイルに書き出したいときは以下を実行
#  thinkplot.Save(root='train1',
#                   xlabel='Number of trains',
#                   ylabel='Probability',
#                   formats=['pdf', 'eps'])


if __name__ == '__main__':
    main()	

■事前確率をべき乗分布で更新
1台の機関車を所有する会社と1000台を所有する会社は同じ確率じゃないだろ
1台の機関車をもつ零細企業は沢山あるのにたいして1000台を所有する大企業は
べき乗分布になるだろう
ということで、事前確率をべき乗分布にした場合
ついでに仮説を500にしたときと1000と2000にしたとき。

import thinkbayes
import thinkplot
from thinkbayes import Pmf, Percentile
from dice import Suite

class Train(Suite):

 def __init__(self, hypos, alpha=1.0):
        Pmf.__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 Mean(suite):
    total = 0
    for hypo, prob in suite.Items():
        total += hypo * prob
    return total

def MakePosterior(high, dataset):
    hypos = xrange(1, high+1)
    suite = Train(hypos)
    suite.name = str(high)

    for data in dataset:
        suite.Update(data)

    thinkplot.Pmf(suite)
    return suite


def main():
    dataset = [60,30,90,130,100,20,30,50,10,40,50,100,20,10,30,40,80,90,50,100]
    for high in [500, 1000, 2000]:
        suite = MakePosterior(high, dataset)
        print high, suite.Mean()

if __name__ == '__main__':
    main()

これからさらに信用区間をとったりするんですが、それは割愛で。

以上

Think Bayes ベイズの定理のメモ 2

前回の続きです。

前回やった問題のコードを書き直して一般的なものにする。

■クッキーの問題
「クッキーの入った2つのボウルがある。
ボウル1 ? バニラクッキー 30枚 :チョコレートクッキー 10枚
ボウル2 ? バニラクッキー 20枚 :チョコレートクッキー 20枚
どちらかのボウルをランダムに選びクッキーを1つランダムに選んだ場合、
それがボウル1だという確率を求めよ。」

from thinkbayes import Pmf
class Cookie(Pmf):
    def __init__(self, hypos):
        Pmf.__init__(self)
        for hypo in hypos:
            self.Set(hypo, 1)
        self.Normalize()

    def Update(self, data):
        for hypo in self.Values():
            like = self.Likelihood(data, hypo)
            self.Mult(hypo, like)
        self.Normalize()

    mixes = {
        'Bowl 1':dict(vanilla=0.75, chocolate=0.25),
        'Bowl 2':dict(vanilla=0.5, chocolate=0.5),
        }

    def Likelihood(self, data, hypo):
        mix = self.mixes[hypo]
        like = mix[data]
        return like


def main():
    hypos = ['Bowl 1', 'Bowl 2']
    pmf = Cookie(hypos)
    pmf.Update('vanilla')
    for hypo, prob in pmf.Items():
        print hypo, prob
if __name__ == '__main__':
    main()

■20年前の早稲田大学の入試問題
「5回に1回の割合で帽子を忘れるくせのあるK君が、正月に A、B、C 3軒を順に年始回りをして家に帰ったとき、
帽子を忘れてきたことに気がついた。2軒目の家 B に忘れてきた確率を求めよ。」

from thinkbayes import Pmf
class Waseda(Pmf):
    def __init__(self, hypos):
        Pmf.__init__(self)
        for hypo in hypos:
            self.Set(hypo, 1)
        self.Normalize()

    def Update(self, data):
        for hypo in self.Values():
            like = self.Likelihood(data, hypo)
            self.Mult(hypo, like)
        self.Normalize()


    def Likelihood(self, data, hypo):
        if hypo == 'A':
            return 1/5.0
        elif hypo == 'B':
            return 4/25.0
        elif hypo == 'C'
            return 16/125.0


def main():
    hypos = ['A', 'B', 'C']
    pmf = Waseda(hypos)
    pmf.Update('B')
    for hypo, prob in pmf.Items():
        print hypo, prob
if __name__ == '__main__':
    main()

■喫煙者の推定の問題
「男性10人、女性7人が一室でパーティーを開いた。男子の喫煙者は5人、女性は3人である。
部屋に入ったら煙草の吸殻が1本、灰皿の上にあった。このとき、吸った人が女性である確率を求めなさい(煙草の吸い回しはしていないと仮定する)」

from thinkbayes import Pmf
class Tabaco(Pmf):
    def __init__(self, hypos):
        Pmf.__init__(self)
        for hypo in hypos:
            self.Set(hypo, 1)
        self.Normalize()

    def Update(self, data):
        for hypo in self.Values():
            like = self.Likelihood(data, hypo)
            self.Mult(hypo, like)
        self.Normalize()

    mixes = {
        'Men':dict(Smoke=1/2.0, NoSmoke=1/2.0),
        'Women':dict(Smoke=3/7.0, NoSmoke=4/7.0),
        }

    def Likelihood(self, data, hypo):
        mix = self.mixes[hypo]
        if hypo == 'Men':
         like = 10/17.0 * mix[data] 
         return like
        elif hypo == 'Women':
           like = 7/17.0 * mix[data] 
           return like


def main():
    hypos = ['Men', 'Women']
    pmf = Tabaco(hypos)
    pmf.Update('Smoke')
    for hypo, prob in pmf.Items():
        print hypo, prob
if __name__ == '__main__':
    main()	

=====フレームワークをカプセル化したとき=====

上のコードの共通化部分をカプセル化したとき。

■クッキーの問題

from thinkbayes import Suite
class Cookie(Suite):
    mixes = {
        'Bowl 1':dict(vanilla=0.75, chocolate=0.25),
        'Bowl 2':dict(vanilla=0.5, chocolate=0.5),
        }

    def Likelihood(self, data, hypo):
        mix = self.mixes[hypo]
        like = mix[data]
        return like

def main():
    suite = Cookie(['Bowl 1', 'Bowl 2'])
    suite.Update('vanilla')
    suite.Print()

if __name__ == '__main__':
    main()

■20年前の早稲田大学の入試問題

from thinkbayes import Suite
class Waseda(Suite):

    def Likelihood(self, data, hypo):
        if hypo == 'A':
            return 1/5.0
        elif hypo == 'B':
            return 4/25.0
        elif hypo == 'C':
            return 16/125.0


def main():
    suite = Waseda(['A', 'B', 'C'])
    suite.Update('B')
    suite.Print()

if __name__ == '__main__':
    main()

■喫煙者の推定の問題

from thinkbayes import Suite
class Smoke(Suite):

    mixes = {
        'Men':dict(Smoke=1/2.0, NoSmoke=1/2.0),
        'Women':dict(Smoke=3/7.0, NoSmoke=4/7.0),
        }

    def Likelihood(self, data, hypo):
        mix = self.mixes[hypo]
        if hypo == 'Men':
         like = 10/17.0 * mix[data] 
         return like
        elif hypo == 'Women':
           like = 7/17.0 * mix[data] 
           return like


def main():
    suite = Smoke(['Men', 'Women'])
    suite.Update('Smoke')
    suite.Print()

if __name__ == '__main__':
    main()	

喫煙者の推定の問題はもっとうまい方法があるかも。
今回は完全にメモ。

*7/24追記
喫煙者の問題で、事前確率を尤度の計算にもってきていたのを
事前確率のところに変更しました。

from thinkbayes import Suite
class Smoke(Suite):

 def __init__(self, hypos):
        Pmf.__init__(self)
        self.Set('Men', 10/17.0)
        self.Set('Women', 7/17.0)
        self.Normalize()

 def Update(self, data):
        for hypo in self.Values():
            like = self.Likelihood(data, hypo)
            self.Mult(hypo, like)
        self.Normalize()
 mixes = {
        'Men':dict(Smoke=1/2.0, NoSmoke=1/2.0),
        'Women':dict(Smoke=3/7.0, NoSmoke=4/7.0),
        }

 def Likelihood(self, data, hypo):
        mix = self.mixes[hypo]
        if hypo == 'Men':
         like = mix[data] 
         return like
        elif hypo == 'Women':
           like = mix[data] 
           return like

def main():
    suite = Smoke(['Men', 'Women'])
    suite.Update('Smoke')
    suite.Print()

if __name__ == '__main__':
    main()	

以上

Think Bayes ベイズの定理のメモ

プログラミングスキルをつかってベイズ統計を理解するという面白そうな本があったので購入してみました。そのメモです。

この本ではthinkbayes.pyで定義されるクラスや関数を使っています。
以下のリンクからthinkbayes.pyというモジュールをダウンロードして使います。
http://thinkbayes.com/thinkbayes.py

・使い方
正六面体のサイコロの分布結果

from thinkbayes import Pmf

pmf = Pmf ()
for x in [1,2,3,4,5,6]:
    pmf.Set(x, 1/6.0)
    
//Setメソッドが各値に確率1/6を設定している。
//各サイコロの目の確率の出力方法は以下

print pmf.Prob(1)
0.166666666667

print pmf.Prob(2)
0.166666666667

print pmf.Prob(3)
0.166666666667

・ベイズの定理

p(H|D) = p(H)p(D|H)/p(D)

p(H) : データを見る前の仮設の確率。事前確率(prior probability)
p(H|D) : データを見た後の仮設の確率。事後確率(posterior probability or posterior)
p(D|H) : 仮説の下でのデータの確率。尤度(likelihood)
p(D) : どのような仮説でもデータの得られる確率。正規化定数(normalizing constant)

・クッキー問題
「クッキーの入った2つのボウルがある。
ボウル1 – バニラクッキー 30枚 :チョコレートクッキー 10枚
ボウル2 – バニラクッキー 20枚 :チョコレートクッキー 20枚
どちらかのボウルをランダムに選びクッキーを1つランダムに選んだ場合、
それがボウル1だという確率を求めよ。」

ステップ1 事前確率の設定-仮説に確率を割り当てる
 仮説B1 ボウル1からバニラクッキーを選ぶ確率 1/2
 仮説B2 ボウル2からバニラクッキーを選ぶ確率 1/2

from thinkbayes import Pmf
pmf = Pmf()
pmf.Set('Bowl 1', 0.5)
pmf.Set('Bowl 2', 0.5)

ステップ2 分布の更新(各事前確率に対応する尤度を掛ける)
 ボウル1からクッキーを取り出す確率 3/4
 ボウル2からクッキーを取り出す確率 1/2

pmf.Mult('Bowl 1', 0.75)
pmf.Mult('Bowl 2', 0.5)

Multメソッドは指定された仮説の仮設を取り出し、そこに指定された尤度を掛ける。

ステップ3 最正規化しボウル1の事後確率を求める。

pmf.Normalize()
print pmf.Prob('Bowl 1')

答え 0.6

このモジュールを使ってほかの問題も解いてみます。
以下は、ネット上で見つけた問題です。

・20年前の早稲田大学の入試問題
「5回に1回の割合で帽子を忘れるくせのあるK君が、正月に A、B、C 3軒を順に年始回りをして家に帰ったとき、
帽子を忘れてきたことに気がついた。2軒目の家 B に忘れてきた確率を求めよ。」
http://d.hatena.ne.jp/pashango_p/20090809/1249805193

ステップ1 事前確率の設定-
A Aに帽子を忘れてきた確率1/3
B Bに帽子を忘れてきた確率1/3
C Cに帽子を忘れてきた確率1/3

ステップ2 
A Aに帽子を忘れて確率(1/5)
B Aで帽子を忘れない確率(4/5)× Bで帽子を忘れる確率(1/5)
C Aで帽子を忘れない確率(4/5)× Bで帽子を忘れない確率(4/5) × Cで帽子を忘れる確率(1/5)

from thinkbayes import Pmf
pmf = Pmf()
pmf.Set('A', 1/3.0)
pmf.Set('B', 1/3.0)
pmf.Set('C', 1/3.0)

pmf.Mult('A', 1/5.0)
pmf.Mult('B', 4/25.0)
pmf.Mult('C', 16/125.0)

pmf.Normalize()
print pmf.Prob('B')

答え 0.327868852459

・喫煙者の推定の問題
「男性10人、女性7人が一室でパーティーを開いた。男子の喫煙者は5人、女性は3人である。
部屋に入ったら煙草の吸殻が1本、灰皿の上にあった。このとき、吸った人が女性である確率を求めなさい(煙草の吸い回しはしていないと仮定する)」
 http://www.yasuienv.net/BayesExamples.htm

1.事前確率
Men  :10/17
Women :7/17

2.分布の更新
Men  :1/2
Women :3/7

pmf = Pmf()
pmf.Set('Men', 10/17.0)
pmf.Set('Women', 7/17.0)

pmf.Mult('Men', 1/2.0)
pmf.Mult('Women', 3/7.0)

pmf.Normalize()
print pmf.Prob('Women')

答え 0.375

・雨の日のお誘いの問題
「酒好きのAさんはB氏をよくお酒に誘う。統計をとると、雨の降っていない日に誘うと、B氏は5回中4回誘いに応じ、
雨の降っている日に誘うと、5回中3回誘いに応じることが分かった。B氏がAさんの誘いに応じたとき、雨が降っていない確率を求めよ。
雨が降った日と降らない日の割合は7:1とする」
http://www.yasuienv.net/BayesExamples.htm

1.事前確率
sunny 7/8
rainy 1/8

2.尤度
sunny 4/5
rainy 3/5

pmf = Pmf()
pmf.Set('sunny', 7/8.0)
pmf.Set('rainy', 1/8.0)

pmf.Mult('sunny', 4/5.0)
pmf.Mult('rainy', 3/5.0)

pmf.Normalize()
print pmf.Prob('sunny')

答え 0.903225806452

このモジュールを使えば簡単に確率の問題が解けることが分かりました。

以上

return top