#Pythonでデータの分析を出来るようになりたい
ちなみに Python 初心者です。
Pythonは文法がすっきりしていて、ライブラリも豊富でパワフルだという評判を聞いてたんで、年末休みを利用してちょっとかじってみることにした。
#題材
何事も、勉強するなら具体的な題材・問題があった方がいい。
そこで最近読み始めた統計に関する本の内容をベースにすることにした。
統計学が最強の学問である[実践編]データ分析のための思想と方法 - 西内啓 著
第1章 統計学の実践は基本の見直しから始まる
05 なぜ、平均値は真実を捉えることができるのか?
この章に「裏が出る確率が2/3で表が出る確率が1/3、というコイン」を使った例がいくつかある。そこで Python を使って実際に「コイン」を何回も投げて(もちろんシミュレーションだけど)例と同じ結果になるか確認してみることにした。
(詳しい内容は本を買うか図書館で借りて調べるかしてください。P59~P64の部分です)
#Pythonのバージョン
使ったPython のバージョンは2.7。夏に Python のセミナーみたいなのを受けたときに言われるがままに入れた Spyderってのがインストールされてたんで、それを使った。
#コインを2回投げた時の確率
まずは、コインを2回投げるときの全ての組み合わせの確率についての例をPythonで実験してみた。
可能な組み合わせは:
裏・裏 (表0枚)
裏・表 (表1枚)
表・裏 (表1枚)
表・表 (表2枚)
##Python コード
from random import randint
from decimal import Decimal
from prettytable import PrettyTable
import numpy as np
def tossBiasedCoin():
""" Returns 0 or 1 with 0 having 2/3 chance """
return randint(0,2) % 2
# Make a 2x2 array
counts = [[0 for j in range(2)] for i in range(2)]
# Toss a coin many times to get counts
sampleCount = 500000
for num in range(sampleCount):
first = tossBiasedCoin()
second = tossBiasedCoin()
counts[first][second] += 1
# Conert all counts to perentage
TWOPLACES = Decimal(10) ** -2
for i in range(2):
for j in range(2):
value = counts[i][j]
counts[i][j] = (100 * Decimal(counts[i][j])/Decimal(sampleCount)).quantize(TWOPLACES)
print("Converted the value {} to percentage {}".format(value, counts[i][j]))
# Make summaries of number of heads.
keys = np.arange(3)
values = [counts[0][0],
counts[0][1]+counts[1][0],
counts[1][1]]
# Add row descriptions
counts[0].insert(0, '1st tail')
counts[1].insert(0, '1st head')
# Create table with column descriptions, add rows, then show it.
table = PrettyTable(["", "2nd tail", "2nd head"])
table.padding_width = 1
table.add_row(counts[0])
table.add_row(counts[1])
print table
# Draw a bar chart
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
rects = plt.bar(keys,
values,
0.5,
alpha=0.4,
align="center",
color='b')
plt.xlabel('Number of heads')
plt.ylabel('Probability (%)')
plt.title('Probabilities heads with a biased coin')
plt.xticks(keys, np.arange(3))
plt.tight_layout()
plt.show()
##実行結果
まずはP59の図表1-17とほぼ同じ表が表示される。4通りの組み合わせと、それぞれの確率。
(この表はPrettyTableを使った。後述)
「表・表」の組み合わせが出にくいのは、そもそも表が出る確率が1/3というコインを使ってるから。
つぎにp60の図表1-18とほぼ同じ棒グラフが表示される。
これは4通りのパターンを「表が何枚出るか」で3つにまとめたもの。「表が1枚」がパターンが(表・裏)と(裏・表)の2通りあったから、それは一つにまとめられている。
#コードの内容
##ライブラリのインポート
使用する関数などを実行環境に import する。
便利なライブラリは Python 標準ライブラリ にあるかもしれないが、膨大な数の third-party ライブラリが Python パッケージインデックス というレポジトリにある。
from random import randint
from decimal import Decimal
from prettytable import PrettyTable
import numpy as np
- 乱数生成にrandomモジュールの randint を使う。
- float型を下2桁に整えるのにdecimalモジュールの Decimal を使う。
- 表を作るにはprettytableモジュールの PrettyTable を使う。
- arange という関数が便利なのでnumpyモジュールをインポートするが、特にヘビーな使い方はしてない。
##Pythonでの関数定義
def tossBiasedCoin():
""" Returns 0 or 1 with 0 having 2/3 chance """
return randint(0,2) % 2
関数にするほどの内容じゃないけど、関数の定義の練習としてコインの表・裏(1か0)を返す関数を作った。0, 1, 2 のいずれかの値を乱数で生成して、値が偶数なら0を返し、それ以外なら1を return する。三つの値のうちふたつが偶数だから2/3の確率になる。
##2x2の変数を準備
出現頻度を記録するために2x2の sequence を使う。この場合は list 型になる。
# Make a 2x2 array
counts = [[0 for j in range(2)] for i in range(2)]
それぞれの変数は0に初期化しておく。
for 文 を、組み込み関数の一つである range関数 で作った即席の list の上でループ実行させている。2x2にするために、list を要素として持つ list を作る。
##コインを投げて結果を記録する
ここでは50万回投げるけど、そんなに投げなくてもいいかな(笑)
# Toss a coin many times to get counts
sampleCount = 500000
for num in range(sampleCount):
first = tossBiasedCoin()
second = tossBiasedCoin()
counts[first][second] += 1
結果は0か1なので、ちょうど2x2構造に対してのインデックスとして使用できる。インデックスしたマスの数を1づつ増やしていく。
ちなみにPythonではC言語でおなじみの ++ オペレータは無いもよう。Pythonのオペレータのリストはこちら。
##パーセントに変換する
頻度を投げた総数で割って、割合を出す。
# Convert all counts to perentage
TWOPLACES = Decimal(10) ** -2
for i in range(2):
for j in range(2):
value = counts[i][j]
counts[i][j] = (100 * Decimal(counts[i][j])/Decimal(sampleCount)).quantize(TWOPLACES)
print("Converted the value {} to percentage {}".format(value, counts[i][j]))
2x2 構造のそれぞれのマスを訪れることになるんで、iとjの二つのインデックスを使ってループ実行し、[i][j]
の形で任意のマスの値にアクセスする。値はin-placeで置き換えるけど、デバッグ用に変換前と変換後の値を表示している。
Decimal.quantize 関数に0.01
を渡すことで値を下2桁に丸める。
##棒グラフ用にデータを準備
棒グラフの棒は3つ。表が0枚、表が1枚、そして表が2枚。
# Make summaries of number of heads.
keys = np.arange(3)
values = [counts[0][0],
counts[0][1]+counts[1][0],
counts[1][1]]
表が1枚の頻度だけは、表がコイン1投目で出たか2投目で出たかに関係ないのでその両方の値を合わせている。
##表用にデータを準備
表の行となる list の左側に list.insert を使って「1投目が裏」と「1投目が表」という欄を加える。
# Add row descriptions
counts[0].insert(0, '1st tail')
counts[1].insert(0, '1st head')
##表を作る
PrettyTableというサード・パーティ・ライブラリ関数を使う。
# Create table with column descriptions, add rows, then show it.
table = PrettyTable(["", "2nd tail", "2nd head"])
table.padding_width = 1
table.add_row(counts[0])
table.add_row(counts[1])
print table
##棒グラフを作る
サード・パーティー製ライブラリ、matplotlibを使う。この matplotlib はそれだけで本が一冊書けそうなぐらい充実しているみたい(ギャラリー参照)
matplotlib の中の pyplot で棒グラフを書く。
# Draw a bar chart
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
rects = plt.bar(keys,
values,
0.5,
alpha=0.4,
align="center",
color='b')
plt.xlabel('Number of heads')
plt.ylabel('Probability (%)')
plt.title('Probabilities heads with a biased coin')
plt.xticks(keys, np.arange(3))
plt.tight_layout()
plt.show()
(その2)につづく