Python言語によるプログラミングイントロダクション 第2版 写経（15章 part 1）

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

の写経をしたので記録を公開しておく。 この記事では、15章のcode 15.9までをまとめておく。

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

写経

import random

import scipy as sp
import scipy.integrate
import scipy.stats

import numpy as np

import statistics as stat

import math

メモ

　matplotlibにはrcParamsという、グラフのフォーマット（フォントサイズ、線の太さ、色等）を設定しているグローバル変数がある。以下のコードのpylab.rcParamsでは、rcParamsのfigure.dpiを400に変更している。デフォルトは100だそうだ。400にすると出力されるグラフが大分滑らかになる。

　変更方法はpylab.rc('figure', dpi=400)と、pylab.rcParam['figure.dpi']=400の二つある。pylab.rcを使用する方は、第一引数にgroupを指定する。この例では、figureというgroupを指定して、figureグループのdpiの値を設定している。また、同グループの別の値を同時に設定することができるようだ。以下の例ではlineグループのlinewidthとcolorの二つの値を設定している。

例：

mpl.rc('lines', linewidth=4, color='g')

参考：Customizing Matplotlib with style sheets and rcParams

import pylab

pylab.rc('figure', dpi=400)

code 15.1

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

def rollN(n):
    result = ''
    for i in range(n):
        result = result + str(rollDie())
    print(result)

rollN(10)

3162435546

code 15.2

def flip(numFlips):
    heads = 0
    for i in range(numFlips):
        if random.choice(('H', 'T')) == 'H':
            heads += 1
    return heads/numFlips

def flipSim(numFlipsPerTrial, numTrials):
    fracHeads = []
    for i in range(numTrials):
        fracHeads.append(flip(numFlipsPerTrial))
    mean = stat.mean(fracHeads)
    return mean

print('Mean = ', flipSim(100, 1000))

Mean =  0.49982

code 15.3

def regressToMean(numFlips, numTrials):
    fracHeads = []
    for t in range(numTrials):
        fracHeads.append(flip(numFlips))
    
    extremes, nextTrials = [], []
    for i in range(len(fracHeads) - 1):
        if fracHeads[i] < 0.33 or fracHeads[i] > 0.66:
            extremes.append(fracHeads[i])
            nextTrials.append(fracHeads[i+1])
    pylab.plot(range(len(extremes)), extremes, 'ko', label = 'Extreme')
    pylab.plot(range(len(nextTrials)), nextTrials, 'k^', label = 'Next Trial')
    pylab.axhline(0.5)
    pylab.ylim(0, 1)
    pylab.xlim(-1, len(extremes) + 1)
    pylab.xlabel('Extreme Example and Next Trial ')
    pylab.ylabel('Fraction Heads ')
    pylab.title('Regression to the Mean ')
    pylab.legend(loc = 'best')
    
regressToMean(15, 40)

code 15.4

def flipPlot(minExp, maxExp):
    ratios, diffs, xAxis = [], [], []
    for exp in range(minExp, maxExp + 1):
        xAxis.append(2**exp)
    for numFlips in xAxis:
        numHeads = 0
        for n in range(numFlips):
            if random.choice(('H', 'T')) == 'H':
                numHeads += 1
        numTails = numFlips - numHeads
        try:
            ratios.append(numHeads/numTails)
            diffs.append(abs(numHeads - numTails))
        except ZeroDivisionError:
            continue
    pylab.title('Difference Between Heads and Tails')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('Abs(#Heads - #Tails)')
    pylab.plot(xAxis, diffs, 'k')
    pylab.figure()
    pylab.title('Heads/Tails Ratios')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('#Heads/#Tails')
    pylab.plot(xAxis, ratios, 'k')
    
random.seed(0)
flipPlot(4, 5)

pylab.figureの参考： matplotlib.pyplot.figure

matplotlib の figure（図） と axes（座標軸）－python

出力する図のサイズや解像度などを設定できるらしい。引数を渡さなければ、rcParamで設定されている値を参照するようだ。

code 15.4

def flipPlot(minExp, maxExp):
    ratios, diffs, xAxis = [], [], []
    for exp in range(minExp, maxExp + 1):
        xAxis.append(2**exp)
    for numFlips in xAxis:
        numHeads = 0
        for n in range(numFlips):
            if random.choice(('H', 'T')) == 'H':
                numHeads += 1
        numTails = numFlips - numHeads
        try:
            ratios.append(numHeads/numTails)
            diffs.append(abs(numHeads - numTails))
        except ZeroDivisionError:
            continue
    pylab.title('Difference Between Heads and Tails')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('Abs(#Heads - #Tails)')
    pylab.plot(xAxis, diffs, 'bo')
    pylab.semilogx()
    pylab.semilogy()
    pylab.figure()
    pylab.title('Heads/Tails Ratios')
    pylab.xlabel('Number of Flips')
    pylab.ylabel('#Heads/#Tails')
    pylab.plot(xAxis, ratios, 'bo')
    pylab.semilogx()
    
random.seed(0)
flipPlot(4, 20)

code 15.5は省略。標準偏差はstat.stdevにより求めることにした。

code 15.8

def CV(X):
    try:
        return stat.stdev(X)/stat.mean(X)
    except ZeroDivisionError:
        return float('nan')

code 15.6, 15.7, 15.9

def makePlot(xVals, yVals, title, xLabel, yLabel, style, logX = False, logY = False):
    pylab.figure()
    pylab.title(title)
    pylab.xlabel(xLabel)
    pylab.ylabel(yLabel)
    pylab.plot(xVals, yVals, style)
    if logX:
        pylab.semilogx()
    if logY:
        pylab.semilogy()

def runTrial(numFlips):
    numHeads = 0
    for n in range(numFlips):
        if random.choice(('H', 'T')) == 'H':
            numHeads += 1
    numTails = numFlips - numHeads
    return (numHeads, numTails)

def flipPlot2(minExp, maxExp, numTrials):
    ratiosMeans, diffsMeans, ratiosSDs, diffsSDs = [], [], [], []
    ratiosCVs, diffsCVs, xAxis = [], [], []
    for exp in range(minExp, maxExp + 1):
        xAxis.append(2**exp)
    for numFlips in xAxis:
        ratios, diffs = [], []
        for t in range(numTrials):
            numHeads, numTails = runTrial(numFlips)
            ratios.append(numHeads/float(numTails))
            diffs.append(abs(numHeads - numTails))
        ratiosMeans.append(stat.mean(ratios))
        diffsMeans.append(stat.mean(diffs))
        ratiosSDs.append(stat.stdev(ratios))
        diffsSDs.append(stat.stdev(diffs))
        ratiosCVs.append(CV(ratios))
        diffsCVs.append(CV(diffs))
    numTrialsString = '(' + str(numTrials) + ' Trials)'
    title = 'Mean Heads/Tails Ratios' + numTrialsString
    makePlot(xAxis, ratiosMeans, title, 'Number of flips', 'Mean Heads/Tails', 'ko', logX = True)
    title = 'SD Heads/Tails Ratios' + numTrialsString 
    makePlot(xAxis, ratiosSDs, title, 'Number of Flips', 'Standard Deviation', 'ko', logX = True, logY = True)
    title = 'Mean abs(#Heads - #Tails)' + numTrialsString
    makePlot(xAxis, diffsMeans, title, 'Number of Flips', 'Mean abs(#Heads - #Tails)', 'ko', logX = True, logY=True)
    title = 'SD abs(#Heads - #Tails)' + numTrialsString
    makePlot(xAxis, diffsSDs, title, 'Number of Flips', 'Standard Deviation', 'ko', logX=True, logY=True)
    title = 'Coeff. of Var. abs(#Heads - #Tails)' + numTrialsString
    makePlot(xAxis, diffsCVs, title, 'Number of Flips', 'Coeff. of Var.', 'ko', logX=True)
    title = 'Coeff. of Var. Heads/Tails Ratio' + numTrialsString
    makePlot(xAxis, ratiosCVs, title, 'Number of Flips', 'Coeff. of Var.', 'ko', logX=True, logY=True)

flipPlot2(4, 20, 20)