　世界標準MIT教科書 Python言語によるプログラミングイントロダクション第2版: データサイエンスとアプリケーション

の写経をしたので記録を公開しておきます。この記事では、16章をまとめておきます。

　Jupyter Labで実行したものをペーストしているので、ちょっとした正誤表になるかもしれません。

　以下、使用モジュール。

import random
import statistics as stat
import pylab
pylab.rc('figure', dpi=400)

写経

code 16.1

def rollDie():
    return random.choice([1, 2, 3, 4, 5, 6])

def checkPascal(numTrials):
    numWins = 0
    for i in range(numTrials):
        for j in range(24):
            d1 = rollDie()
            d2 = rollDie()
            if d1 == 6 and d2 == 6:
                numWins += 1
                break
    print('Probability of winning = ', numWins/numTrials)

checkPascal(1000000)

Probability of winning =  0.49219

code 16.2

class CrapsGame(object):
    def __init__(self):
        self.passWins, self.passLosses = 0, 0
        self.dpWins, self.dpLosses, self.dpPushes = 0, 0, 0
        
    def playHand(self):
        throw = rollDie() + rollDie()
        if throw == 7 or throw == 11:
            self.passWins += 1
            self.dpLosses += 1
        elif throw == 2 or throw == 3 or throw == 12:
            self.passLosses += 1
            if throw == 12:
                self.dpPushes += 1
            else:
                self.dpWins += 1
        else:
            point = throw
            while True:
                throw = rollDie() + rollDie()
                if throw == point:
                    self.passWins += 1
                    self.dpLosses += 1
                    break
                elif throw == 7:
                    self.passLosses += 1
                    self.dpWins += 1
                    break
    
    def passResults(self):
        return (self.passWins, self.passLosses)
    
    def dpResults(self):
        return (self.dpWins, self.dpLosses, self.dpPushes)

code 16.3

def crapSim(handsPerGame, numGames):
    games = []
    for t in range(numGames):
        c = CrapsGame()
        for i in range(handsPerGame):
            c.playHand()
        games.append(c)
    
    pROIPerGames, dpROIPerGames = [], []
    for g in games:
        wins, losses = g.passResults()
        pROIPerGames.append((wins - losses)/float(handsPerGame))
        wins, losses, pushes = g.dpResults()
        dpROIPerGames.append((wins - losses)/float(handsPerGame))
    
    meanROI = str(round(100*stat.mean(pROIPerGames)/numGames, 4)) + '%'
    sigma = str(round(100*stat.stdev(pROIPerGames), 4)) + '%'
    print('Pass:', 'Mean ROI =', meanROI, 'Std. Dev. =', sigma)
    meanROI = str(round(100*stat.mean(dpROIPerGames)/numGames, 4)) + '%'
    sigma = str(round(100*stat.stdev(dpROIPerGames), 4)) + '%'
    print('Don\'t pass:', 'Mean ROI =', meanROI, 'Std. Dev. =', sigma)

crapSim(20, 20)

Pass: Mean ROI = 0.15% Std. Dev. = 18.0933%
Don't pass: Mean ROI = -0.3125% Std. Dev. = 19.3904%

crapSim(1000000, 10)

Pass: Mean ROI = -0.1492% Std. Dev. = 0.1136%
Don't pass: Mean ROI = -0.1286% Std. Dev. = 0.1083%

crapSim(20, 1000000)

Pass: Mean ROI = -0.0% Std. Dev. = 22.3581%
Don't pass: Mean ROI = -0.0% Std. Dev. = 22.044%

code 16.4

class CrapsGame(object):
    def __init__(self):
        self.passWins, self.passLosses = 0, 0
        self.dpWins, self.dpLosses, self.dpPushes = 0, 0, 0
        
    def playHand(self):
        pointsDict = {4: 1/3, 5: 2/5, 6: 5/11, 8: 5/11, 9: 2/5, 10: 1/3}
        throw = rollDie() + rollDie()
        if throw == 7 or throw == 11:
            self.passWins += 1
            self.dpLosses += 1
        elif throw == 2 or throw == 3 or throw == 12:
            self.passLosses += 1
            if throw == 12:
                self.dpPushes += 1
            else:
                self.dpWins += 1
        else:
            if random.random() <= pointsDict[throw]:
                self.passWins += 1
                self.dpLosses += 1
            else:
                self.passLosses += 1
                self.dpWins += 1
    
    def passResults(self):
        return (self.passWins, self.passLosses)
    
    def dpResults(self):
        return (self.dpWins, self.dpLosses, self.dpPushes)

def rollDie():
    return random.choice([1, 1, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6])

crapSim(1000000, 10)

Pass: Mean ROI = 0.2299% Std. Dev. = 0.0927%
Don't pass: Mean ROI = -0.508% Std. Dev. = 0.0978%

code 16.5

def throwNeedles(numNeedles):
    inCircle = 0
    for Needles in range(numNeedles): # テキストだと、1, 1+numNeedlesって書いてあったけど、for文の中でNeedles変数を使用するところなさそうだし、これでいいはず。
        x = random.random()
        y = random.random()
        if (x*x + y*y)**0.5 <= 1:
            inCircle += 1
    return 4*(inCircle/float(numNeedles))

def getEst(numNeedles, numTrials):
    estimate = []
    for t in range(numTrials):
        piGuess = throwNeedles(numNeedles)
        estimate.append(piGuess)
    sDev = stat.stdev(estimate)
    curEst = stat.mean(estimate)
    print('Est. =' + str(round(curEst, 5)) + ',', 'Std. dev. = ' + str(round(sDev, 5)) + ',', 'Needles = ' + str(numNeedles)) # テキストはprint関数の括弧が閉じてない気がする。
    return (curEst, sDev)

def estPi(precision, numTrials):
    numNeedles = 1000
    sDev = precision
    while sDev >= precision/1.96:
        curEst, sDev = getEst(numNeedles, numTrials)
        numNeedles *= 4 # じれったいので4倍で増やすことにした
    return curEst

estPi(0.01, 100)

Est. =3.13016, Std. dev. = 0.05046, Needles = 1000
Est. =3.14358, Std. dev. = 0.03472, Needles = 2000
Est. =3.14173, Std. dev. = 0.02282, Needles = 4000
Est. =3.14166, Std. dev. = 0.01758, Needles = 8000
Est. =3.14183, Std. dev. = 0.01209, Needles = 16000
Est. =3.14296, Std. dev. = 0.00865, Needles = 32000
Est. =3.14094, Std. dev. = 0.0062, Needles = 64000
Est. =3.14217, Std. dev. = 0.00492, Needles = 128000





3.142170625

curEst ± precisionの線をプロットしたヒストグラムを書いてみた。確かに、平均の近くに95%以上のデータが含まれていそうだということが確認できた。

def throwNeedles(numNeedles):
    inCircle = 0
    for Needles in range(numNeedles): # テキストだと、1, 1+numNeedlesって書いてあったけど、for文の中でNeedles変数を使用するところなさそうだし、これでいいはず。
        x = random.random()
        y = random.random()
        if (x*x + y*y)**0.5 <= 1:
            inCircle += 1
    return 4*(inCircle/float(numNeedles))

def getEst2(numNeedles, numTrials, precision):
    """precisionはestPi関数から参照することを前提としている"""
    estimate = []
    for t in range(numTrials):
        piGuess = throwNeedles(numNeedles)
        estimate.append(piGuess)
    sDev = stat.stdev(estimate)
    curEst = stat.mean(estimate)
    title = 'Est. =' + str(round(curEst, 5)) + ',' + 'Std. dev. = ' + str(round(sDev, 5)) + ',' + 'Needles = ' + str(numNeedles)
    pylab.figure()
    pylab.title(title)
    pylab.hist(estimate)
    pylab.axvline(curEst - precision)
    pylab.axvline(curEst)
    pylab.axvline(curEst + precision)
    return (curEst, sDev)

def estPi(precision, numTrials):
    numNeedles = 1000
    sDev = precision
    while sDev >= precision/1.96:
        curEst, sDev = getEst2(numNeedles, numTrials, precision=precision)
        numNeedles *= 4 # じれったいので16倍で増やすことにした
    return curEst