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()

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

以上

関連する投稿:

• 2016-07-18:Think Bayes　ベイズの定理のメモ 2
• 2016-07-09:Think Bayes　ベイズの定理のメモ
• 2017-01-21:Yahoo!Japan ファイナンス から株価の取得
• 2017-02-01:GDPを支える要素のメモとグラフ
• 2020-04-11:新型コロナウイルス（covid-19）陽性患者数のグラフ 1