Blogs

Pythonプログラミングイントロダクション(17章) まで読みました。 17章では中心極限定理について説明されています。

この章の目標の一つは「一つの標本から標本平均の信頼度を出す」ことにあります。 そのための道具として中心極限定理を利用します。

中心極限定理

中心極限定理についての引用は下記になります。

- 同じ母集団から抽出された標本サイズが十分に大きいと、標本の平均(標本平均)はおおよそ正規分布に従う
- この正規分布の平均は母集団の平均に非常に近い。
- 標本平均の分散は母集団の分散を標本サイズで割ったものに非常に近い

中心極限定理により標本平均は正規分布に従います。

また、標本平均の分散は母集団の分散を標本サイズで割ったものに非常に近いという事は、標本平均の標準偏差は母集団の標準偏差を標本サイズの平方根で割ったものに非常に近いということになります。 これを式にすると、標本平均の標準偏差を求める式は下記のようになります。

母集団の標準偏差を標本サイズの平方根で割ったものは標準誤差とも言われます。 そしてこのとき必要な母集団の標準偏差は、ある一定数の標本サイズがあれば標本の標準偏差を用いて良いとされているそうです。

まとめると、下記の手順で中心極限定理が正しいかを見ていきます。

1. 標本の値から標準偏差を求め、母集団の標準偏差とする。
2. 標本平均の標準偏差を求める。
3. この標準偏差から95%信頼区間を導きます。
4. 実際の母集団の平均値と標本平均の誤差を求める

アニメのレーテイングデータで検証する

kaggleのサイトでアニメのレーティングに関するデータを利用させていただきました。 本書と同様の手順でこのデータから標本サイズ200のデータを元に信頼区間を導き、 実際の母集団から求めた平均との絶対誤差を元に信頼区間の値が妥当を見ていきます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

import random
import pylab
import csv

def variance(X):
    """分散"""
    mean = sum(X)/len(X)
    tot = 0.0
    for x in X:
        tot += (x - mean)**2
    return tot/len(X)


def stdDev(X):
    """標準偏差(分散の平方根)"""
    return variance(X)**0.5

def read_anime_score(path):
    """データを読み込みスコアのリストを返す"""
    scores = []
    with open(path, newline='') as f:
        reader = csv.reader(f)
        # skip header
        next(reader)
        for row in reader:
            scores.append(float(row[15]))

    return scores

data = read_anime_score('./anime_cleaned.csv')

# 母平均
popMean = sum(data) / len(data)
# 標本サイズ
sampleSize = 200
numBad = 0
for t in range(10000):
    sample = random.sample(data, sampleSize)
    # 標本平均
    sampleMean = sum(sample) / sampleSize
    # 標準誤差
    se = stdDev(sample) / sampleSize**0.5

    # 母平均と標本平均の誤差から信頼区間外のデータ数を求める
    if abs(popMean - sampleMean) > 1.96 * se:
        numBad += 1
print('Fraction outside 95% confidence interval =', numBad / 10000)

: Fraction outside 95% confidence interval = 0.05

概ね5%程度の結果が求められます。 中心極限定理で定義されている通り、標本サイズが十分大きいと標本平均が正規分布に従い、母集団の平均に非常に近いと言えそうです。

まとめ

本章は話の流れがわかりにくかったので最初は何を言っているのか正直わかりませんでしたが、 まとめてみると最終的に標本平均から信頼区間を求めたいのだとわかりました。

母集団の標準偏差が標本の標準偏差で代用できる点は気持ちの悪いところですが、本書でも下記のように説明されています。

多くの統計学者が言うには、母集団の分布がおおよそ正規分布であれば標本サイズが30から40で十分に大きい。
標本サイズが小さい場合にはt分布と呼ばれるものを使って区間計算したほうが良い。

t分布はこの後出てくるようなので一旦読み進めていくとします。

Pythonプログラミングイントロダクション(16章) まで読みました。 16章ではモンテカルロ・シミュレーションについての説明になります。

モンテカルロ・シミュレーションって大そうな名前がついていますが、この章までにやってきたランダムネスを利用したシミュレーションのことをモンテカルロ・シミュレーションといいます。 この章でもサイコロゲームを題材にランダムネスを利用したモンテカルロシミュレーションの実装がいくつか例示されています。

ランダムネスで円周率を求めるプログラミングなんかは痺れました。

ジャンケンのシミュレーション

ここではジャンケンのモンテカルロ・シミュレーションをPythonで実装して理解を深めていきたいと思います。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

import pylab
import random


def janken():
    JANKEN_PATTERN = ['G', 'C', 'P']
    while True:
        a_t = random.choice(JANKEN_PATTERN)
        b_t = random.choice(JANKEN_PATTERN)
        if a_t != b_t:
            break
    return janken_judge(a_t, b_t)


def janken_kai():
    JANKEN_WIN_PATTERN = ['G', 'G', 'C', 'P']
    JANKEN_LOSE_PATTERN = ['G', 'C', 'C', 'P']
    while True:
        a_t = random.choice(JANKEN_WIN_PATTERN)
        b_t = random.choice(JANKEN_LOSE_PATTERN)
        if a_t != b_t:
            break
    return janken_judge(a_t, b_t)


def janken_judge(a_t, b_t):
    """ a_tとb_tは違う値の想定。
    aが勝っていたらTrue, 負けていたらFalseを返す """

    win_pattern = [('G', 'C'), ('C', 'P'), ('P', 'G')]
    return (a_t, b_t) in win_pattern


def show_plot(trials, target):
    a_rate = []
    a_win = 0
    for k in range(1, trials + 1):
        if target():
            a_win += 1

        a_rate.append(a_win / k * 100)
    pylab.plot(a_rate)
    pylab.ylim(0, 100)
    pylab.xlim(1, trials)
    pylab.xlabel('times')
    pylab.ylabel('rate')
    pylab.title('Junken')
    pylab.show()


def get_stdDev(trials):
    a_rate = []
    for i in range(0, 100):
        a_win = 0
        for k in range(0, trials):
            if janken():
                a_win += 1
        a_rate.append(a_win / k * 100)
    return stdDev(a_rate)


def variance(X):
    """分散"""
    mean = sum(X)/len(X)
    tot = 0.0
    for x in X:
        tot += (x - mean)**2
    return tot/len(X)


def stdDev(X):
    """標準偏差(分散の平方根)"""
    return variance(X)**0.5


show_plot(1000, janken)
show_plot(1000, janken_kai)
print('10 trials stdDev: ', get_stdDev(10))
print('100 trials stdDev: ', get_stdDev(100))
print('1000 trials stdDev: ', get_stdDev(1000))

janken関数では、a,bそれぞれのグー(G)、チョキ(C)、パー(P)の中からランダムに選択した上で勝ち負けを判定しています。 あいこの場合はもう一度ランダムに選択しなおします。

完全にランダムにグー、チョキ、パーを選ぶ場合は試行回数が増えるについれて勝率は5割に近づいていきます。 15章で学んだ大数の法則の通りですね。

人間がジャンケンをする場合は完全にランダムとはいかず、グー、チョキ、パーのどれかに選択が偏ることがあるでしょう。 ランダムではありますが、aさんがグーを出しがち、bさんがチョキを出しがちな場合はaさんの勝率が上がっていきます。そんなjanken_kaiを作ってみました。 この場合はaさんの勝率が5割より上回っていますね。

続いて標準偏差を出してみました。 試行回数を増やすにつれ、標準偏差が減っていきます。

: 10 trials stdDev:  17.91612832955234
: 100 trials stdDev:  5.337158102067871
: 1000 trials stdDev:  1.6880153857382443

まとめ

これまで学んだ統計の知識を織り交ぜつつ、ジャンケンを題材にモンテカルロ・シミュレーションをやってみました。 ここ何章かは思ったより統計学の内容が多いことに若干戸惑いもありましたが、内容も楽しいですしトレンド的にも良い題材なのかと思います。 Pythonを学ぶついでに統計学も学べるお得な一冊だと思います。

Pythonプログラミングイントロダクション(15章)まで読みました。 15章では確率統計を題材に、正規分布や関連する話題として信頼区間について述べられています。

様々な確率統計の題材について、Pythonでの実装とpylabを用いたグラフ表示付きで説明されているので、写経をしていくことで確率統計の知識が深まります。 その中でも確率分布についてみていこうと思います。

様々な確率分布

この章では様々な確率分布を例にPythonでの実装が例示されています。 この章で学ぶ確率分布の一部を列挙してみました。

2項分布の練習問題

上記にあげた確率分布の中でも、2項分布については練習問題がついていたのでみていこうと思います。 練習問題は下記になります。

k回の公平なサイコロ投げでちょうど3が2回出る確率を計算する関数を実装せよ。
この関数を用いて、kを2から100まで変化させたときの確率をプロットせよ。

まず、サイコロを投げて3が出る確率は 1/6 になります。 k回サイコロを投げて2回しか3がでず、残りは3以外の数字が出るということになります。 k回サイコロを投げて3が2回出る組み合わせを考える必要があります。

k = 2 のときは

k = 3 のときは、1,2回目が3の時、1,3回目が3の時、2,3回目が3の時の3通りの確率を足し合わせます。

このような2項分布と呼ばれる確率分布の式は以下で表されます。

Pythonでの実装は下記のようになります。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import math
import pylab

def binomical_distribution(k, p, t):
    c = math.factorial(k) / (math.factorial(t) * math.factorial(k - t))
    return c * p**t * (1 - p)**(k - t)


def show_plot(p, t):
    x, y = [], []
    for k in range(t, 101):
	x.append(k)
	y.append(binomical_distribution(k, p, t))

    pylab.plot(x)
    pylab.plot(y)
    pylab.ylim(0, 0.5)
    pylab.xlim(2, 100)
    pylab.xlabel('k')
    pylab.ylabel('p')
    pylab.title('Binomical Distribution')
    pylab.show()

show_plot(1 / 6, 2)

まとめ

2項分布の指練習について、実際にPythonで実装してみました。 この章で学んだ信頼区間についてはこの後でも深掘りされていくようなので楽しみです。

「統計学が最強の学問である(数学編)」を読み終わりました。機械学習を学ぶための数学の知識を叩き込んで、機械学習の入り口まで連れて行くぞという著者の意志を強く感じる本でした。

記号論理学から始まって、連立方程式、線形代数、微分積分まで、幅広い知識の中から機械学習に必要なエッセンスを詰め込みながら、機械学習・ディープラーニングの入り口まで連れて行ってくれる一冊になっています。

内容

この本の流れを簡単にまとめると下記のようになります。

1. 簡単な確率、記号論理学、集合から最終的には迷惑メールフィルタでも用いられていたナイーブベイズ
2. 1次、2次の関数を学び、ちょうど良い数を求めるための平方完成
3. 二項定理から二項分布、対数、ネイピア数について
4. たくさんの数を効率的に表記するためにΣ、ベクトル、行列
5. ちょうど良いところを探すための微分、確率密度関数を理解するための積分
6. 線形代数と微積分を同時に使うことで機械学習、ディープラーニングの入り口へ！

やはりクライマックスはこれまでの知識を積み上げて機械学習の入り口が見えてくる6章です。特に行列、偏微分のありがたみがわかった最小二乗法の解法はしびれました。

行列、偏微分を用いた最小二乗法とnumpy

2章では訪問回数と契約件数のデータを元に比例関係の式を最小二乗法で2章で学んだ最小二乗法は平方完成により答えをめました。 6章ではこれまで学んできた知識、つまり微分、行列を使った解法を解説してくれています。

下記が2章の例示されていたデータになります。

上記の表で表されるモデルが訪問回数をx, 契約件数をyとした場合に直線で表されるとします。 この場合にモデルを表す式は下記のようになります。

εが誤差でこの誤差を最小にするa, bを求めます。 データは3つありますので、それぞれx, yに代入すると3つの式が出てくるわけですが、これを行列で表現すると3つの式ではなく一つの行列式で表すことができます。

\[ y = X\beta + \varepsilon \]

X, y, β,εはそれぞれ下記のような行列になります。

εの二乗和が最小になるようなa, とbを求めるのが最小二乗法です。 行列の場合、誤差の二乗の和を求める方法としてはεの転置行列とεの内積で表現することができます。

\[ \varepsilon^{ \mathrm { T } }\varepsilon \]

この二乗和が最小になる時のa,bを求めるためにβで偏微分してやります。

\[ \frac{ \partial \varepsilon^{ \mathrm { T } }\varepsilon } { \partial \beta } \]

ここからは、この章にたどり着くまでに得た行列や偏微分の知識を駆使して進めて行くのですが、 詳しい解法は是非本書で確認してみてください。結果として下記の式が得られます。

\[ \beta = (X^{ \mathrm { T } }X)^{-1}X^{ \mathrm { T } }y \]

\[ (X^{ \mathrm { T } }X)\beta = X^{ \mathrm { T } }y \]

この式に今回の例を代入すると、下記のようになります。

ここまでくるとPythonで機械学習をやる上で定番のnumpyを使い、逆行列を用いて傾きと切片を求めることができます。 この本で学んだ数学の知識が、自分がこれまで学んできた機械学習の定番、numpyと出会った瞬間でした。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

import numpy

x = numpy.array([[1, 1], [1, 2], [1, 3]])
y = numpy.array([10, 40, 82])

left = numpy.dot(x.T , x)
right = numpy.dot(x.T,  y)
inv_left = numpy.linalg.inv(left)
beta = numpy.dot(inv_left, right)
print(beta)

結果は下記のようになります。

1

[-28.  36.]

データ(今回の場合x, y)が増えていっても、この式に当てはめれば容易に最小二乗法を用いたモデルを導くことができます。

まとめ

自分自身はこの本を読む前に機械学習の本を読み始めてしまっていましたが、これから学習する方には是非ともお勧めしたい一冊です。 この本で得た知識を携えて、なんとなく通り過ぎているこれまでの知識をもう一度ブラッシュアップしたいと思います。

引き続きPythonプログラミングイントロダクション(12章)をみていきます。 この章ではナップザック問題を題材に貪欲アルゴリズムについて説明がされていますが、貪欲アルゴリズムはハフマン符号という符号化にも利用されています。 そして、ハフマン符号はJPEGやZIPで利用されている符号化です。

JPEGやZIPのような馴染みのある要素技術として利用されているハフマン符号とは何か、貪欲アルゴリズムの使いどころはどこなのかをみていきます。

ハフマン符号

ある文章の情報を表現するときにそれぞれの文字に2進数を割り当てたとしたとします。 文章中によく出る文字には短い2進数の割り当てを行い、あまり出てこない文字には長い2進数の割り当てを行うと、全体として短い2進数で文章を表現できるかもしれません。 これがハフマン符号です。

簡単な例として、a,b,c,d,eの文字が文章内に出てくる頻度を表にしました。 頻度の一番多い b には短い2進数 11 が割り当てられ、頻度の一番少ない a には長い2進数 000 が割り当てられています。

Pythonによる実装

この割り当てを求めるためのPythonによる実装は下記になります。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

class Node():

    def __init__(self, node, freq):
        self.freq = freq
        self.node_list = []
        self.node_list.append(node)

    def extend(self, node):
        self.node_list.extend(node.node_list)

    def __lt__(self, other):
        return self.freq < other.freq

def huffman(nodes):
    from heapq import heappush, heappop
    huffmanQueue = []
    ans = {}

    # 優先度つきキューにノードを格納する
    for node in nodes:
        heappush(huffmanQueue, (node.freq, node))

    while len(huffmanQueue) > 1:
        (leftFreq, leftNode) = heappop(huffmanQueue)
        (rightFreq, rightNode) = heappop(huffmanQueue)
        # 一番頻度が少ないノードに'0'を追加
        for s in leftNode.node_list:
            ans[s] = '0' + ans.get(s, '')

        # 二番目に頻度が少ないノードに'1'を追加
        for s in rightNode.node_list:
    	    ans[s] = '1' + ans.get(s, '')

        totalFreq = leftFreq + rightFreq
        leftNode.extend(rightNode)
        heappush(huffmanQueue, (totalFreq, leftNode))
    return ans

print(huffman([Node('a', 5), Node('b', 35), Node('c', 25), Node('d', 15), Node('e', 30)]))

Nodeクラスは頻度と文字を管理するクラスです。文字はそのノード配下にある全ての文字を管理できるようにリストで配下の文字を管理しています。 まず出現頻度が一番低い文字と二番目に低い文字を探し出します。この時、前回も利用した優先度付きキューを利用して出現頻度が低い順にpopしています。 最終的に割り当てられる二進数文字はディクショナリで管理しています。一番目に出現頻度が低い文字には'0’を、二番目に頻度が少ないノードには'1’を追加していきます。

この2つのノードについて頻度を足し合わせた数を求め、ノードも結合した物を再度優先度付きキューに格納します。 改めてこの優先度付きキューについて出現頻度が低いものから順に同様の処理を行います。このように、頻度が低いノードをひたすら処理していくことで、割り当て文字を求めていきます。

実行結果は下記のようになります。

: {'a': '000', 'd': '001', 'c': '01', 'e': '10', 'b': '11'}

図で表すと下記のようになります。最初に a , d が処理され、 0, 1 が割り振られます。2つの出現頻度を足し合わせた 20 と二つのノードを結合させたものを優先度付きキューに再度格納します。 続いてこの 20 の出現頻度のノード ( a と d が結合したもの ) と c について処理を行います。 a , d には 0 が追加で割り振られ、 c には 1 が割り振られます。あとは同様に頻度を足し合わせ、ノードを結合していきながら繰り返し処理を行なっています。

まとめ

貪欲アルゴリズムを題材にハフマン符号について取り上げました。 そこまで複雑ではないアルゴリズムが、普段利用しているファイル形式の内部で利用されているということでテンションの上がる内容でした。