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

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

の写経をしたので記録を公開しておきます。 この記事では、15章のcode 15.9より後ろをまとめておきます。

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

vals = []
for i in range(1000):
    num1 = random.choice(range(0, 101))
    num2 = random.choice(range(0, 101))
    vals.append(num1 + num2)
pylab.hist(vals, bins = 10)
pylab.xlabel('Number of Occurrences')

    Text(0.5, 0, 'Number of Occurrences')

code 15.10

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

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

def labelPlot(numFlips, numTrials, mean, sd):
    pylab.title(str(numTrials) + ' trials of ' + str(numFlips) + ' flips each')
    pylab.xlabel('Fraction of Heads')
    pylab.ylabel('Number of Trials')
    pylab.annotate('Mean = ' + str(round(mean, 4)) + '\nSD = ' + str(round(sd, 4)), size='x-large', xycoords='axes fraction ', xy=(0.67, 0.5))

def makePlots(numFlips1, numFlips2, numTrials):
    val1, mean1, sd1 = flipSim(numFlips1, numTrials)
    pylab.hist(val1, bins=20)
    xmin, xmax = pylab.xlim()
    labelPlot(numFlips1, numTrials, mean1, sd1)
    pylab.figure()
    val2, mean2, sd2 = flipSim(numFlips2, numTrials)
    pylab.hist(val2, bins=20)
    pylab.xlim(xmin, xmax)
    labelPlot(numFlips2, numTrials, mean2, sd2)

makePlots(100, 1000, 100000)

print(sp.integrate.quad(abs, 0, 5)[0])

12.5

以下のコードのsp.stats.norm.pdfは scipy.stats.norm を参照した。テキストのgaussianはなんだかわからない。

print(sp.integrate.quad(sp.stats.norm.pdf, -2, 2, (0, 1))[0])

0.9544997361036417

code 15.11

def checkEmprical(numTrials):
    for t in range(numTrials):
        mu = random.randint(-10, 10)
        sigma = random.randint(1, 10)
        print('For mu =', mu, 'and sigma =', sigma)
        for numStd in (1, 2, 3):
            area = sp.integrate.quad(sp.stats.norm.pdf, mu - numStd*sigma, mu + numStd*sigma, (mu, sigma))[0]
            print('    Fraction within', numStd, 'std =', round(area, 4))

checkEmprical(3)

For mu = 6 and sigma = 4
    Fraction within 1 std = 0.6827
    Fraction within 2 std = 0.9545
    Fraction within 3 std = 0.9973
For mu = 5 and sigma = 10
    Fraction within 1 std = 0.6827
    Fraction within 2 std = 0.9545
    Fraction within 3 std = 0.9973
For mu = -10 and sigma = 1
    Fraction within 1 std = 0.6827
    Fraction within 2 std = 0.9545
    Fraction within 3 std = 0.9973

code 15.12

def showErrorBars(minExp, maxExp, numTrials):
    means, sds, xVals = [], [], []
    for exp in range(minExp, maxExp + 1):
        xVals.append(2**exp)
    for numFlips in xVals:
        fracHeads, mean, sd = flipSim(numFlips, numTrials)
        means.append(mean)
        sds.append(sd)
    pylab.errorbar(xVals, means, yerr=1.96*pylab.array(sds))
    pylab.semilogx()
    pylab.title('Mean Fraction of Heads (' + str(numTrials) + 'trials)')
    pylab.xlabel('Number of flips per trial')
    pylab.ylabel('Fraction of heads & 95% confidence ') # テキストだとインデントが入っていないが、おそらく、この行も関数内だろう。

showErrorBars(3, 10, 100)

指練習(224ページ)

def diceThrowPlot(min=2, max=100):
    """
    minは2以上の整数
    maxはmin以上の整数
    """
    probabilities = []
    xAxis = []
    for t in range(min, max + 1):
        probabilities.append(math.factorial(t)/(math.factorial(2)*math.factorial(t-2)))
        xAxis.append(t)
    pylab.title('Finger practice')
    pylab.xlabel('Throw number')
    pylab.ylabel('Probability of two times success')
    pylab.plot(xAxis, probabilities)

diceThrowPlot()

code 15.13

def clear(n, p, steps):
    numRemaining = [n]
    for t in range(steps):
        numRemaining.append(n*(1-p)**t)
    pylab.plot(numRemaining)
    pylab.xlabel('Time')
    pylab.ylabel('Molecules Remaining')
    pylab.title('Crearance of Drug')

clear(1000, 0.01, 1000)

def clear(n, p, steps):
    numRemaining = [n]
    for t in range(steps):
        numRemaining.append(n*(1-p)**t)
    pylab.plot(numRemaining)
    pylab.xlabel('Time')
    pylab.ylabel('Molecules Remaining')
    pylab.semilogy()
    pylab.title('Crearance of Drug')

clear(1000, 0.01, 1000)

ここは、random.expovariateで遊んでみた。確かに、random.expovariateは指数分布に従う確率変数になっているようだ。

random.expovariate(10)

0.010019297929931845

def plotExpDist(length=10000, lambd=1):
    samples = []
    for i in range(length):
        samples.append(random.expovariate(lambd))
    pylab.hist(samples)

plotExpDist()

code 15.14

def successfulStarts(successProb, numTrials):
    """
    eventProb：1回の試行で成功する確率を表す浮動小数点 
    numTrials：正の整数
    各実験において、成功するまでに必要な試行の回数を出力する。
    """
    # 上のdocstringのeventProbはsuccessProbのことのはず。
    triesBeforeSuccess = []
    for t in range(numTrials):
        consecFailures = 0
        while random.random() > successProb:
            consecFailures += 1
        triesBeforeSuccess.append(consecFailures)
    return triesBeforeSuccess

successProb  = 0.5 # テキストではprobOfSuccessだったが、successfulStats関数を作成したときの引数名であるsuccessProbにそろえた。
numTrials = 5000
distribution = successfulStarts(successProb, numTrials)
pylab.hist(distribution, bins = 14)
pylab.xlabel('Tries Before Success')
pylab.ylabel('Number of Occurrences Out of ' + str(numTrials))
pylab.title('Probability of Starting Each Try = ' + str(successProb))

    Text(0.5, 1.0, 'Probability of Starting Each Try = 0.5')

code 15.15

def collisionProb(n, k):
    prob = 1.0
    for i in range(1, k):
        prob = prob * ((n-i)/n)
    return 1 - prob

def simInsertions(numIndices, numInsertions):
    choices = range(numIndices)
    used = []
    for i in range(numInsertions):
        hashVal = random.choice(choices)
        if hashVal in used:
            return 1
        else:
            used.append(hashVal)
    return 0

def findProb(numIndices, numInsertions, numTrials):
    collisions = 0
    for t in range(numTrials):
        collisions += simInsertions(numIndices, numInsertions)
    return collisions/numTrials

print('Actual probability of a collision =', collisionProb(1000, 50))
print('Est. probability of a collision =', findProb(1000, 50, 10000))
# Est.はEstimated（推計の）の略。findProbはシミュレーション（推計）の結果を返す関数だった。
print('Actual probability of a collision =', collisionProb(1000, 200))
print('Est. probability of a collision =', findProb(1000, 200, 10000))

Actual probability of a collision = 0.7122686568799875
Est. probability of a collision = 0.7199
Actual probability of a collision = 0.9999999994781328
Est. probability of a collision = 1.0

code 15.16

def playSeries(numGames, teamProb):
    numWon = 0
    for game in range(numGames):
        if random.random() <= teamProb:
            numWon += 1
    return (numWon > (numGames//2))

def fractionWon(numGames, teamProb, seriesLen): #引数の名前がわかりにくかったので、playSeries関数のものと揃えた。
    won = 0
    for series in range(seriesLen): 
        if playSeries(numGames, teamProb):
            won += 1
    return won/(float(seriesLen))

def simSeries(seriesLen):
    teamProb = 0.5
    fracsWon, probs = [], []
    while teamProb <= 1.0:
        fracsWon.append(fractionWon(7, teamProb, seriesLen))
        probs.append(teamProb)
        teamProb += 0.01
    pylab.axhline(0.95) #Draw line at 95%
    pylab.plot(probs, fracsWon, 'k', linewidth=5)
    pylab.xlabel('Probability of Winning a Game')
    pylab.ylabel('Probability of Winning a Series')
    pylab.title(str(seriesLen) + ' Seben-Game Series')

simSeries(400)

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