2018年1月5日金曜日

【今日の読書】 R によるやさしい統計学

書名: R によるやさしい統計学
著者: 山田 剛史、杉澤 武俊、 村井 潤一郎
出版: オーム社
ISBN: 978-4-274-06710-5

題名の通り、R 言語を用いて統計学の基礎を説明している本です。

統計については基本的な検定 (t 分布を用いた検定、カイ二乗検定、無相関検定、t 検定) や分散分析あたりまでが対象です。

R については Windows と MacOS でインストールから解説していますが、いかんせん刊行が 2008 年と古いので現在の環境とは乖離しているかもしれません。(僕は R ではなく Python で実装をすすめたので R の事は不明です)



2017年12月19日火曜日

「R によるやさしい統計学」を Python で実習する 第 4 章

第 4 章 母集団と標本

2、3 章でやったのは、母集団であれ標本であれ手元にあるデータの集合を統計量を使って記述するという事。 この章では、大きな集団 (母集団) から取り出した幾つかの標本を元に、母集団の性質について推測する推測統計。

あとこの辺まで読んでびっくりしたのが、この本の中で「心理統計学の基礎 統合的理解のために」が何度も取り上げられている事。

In [17]:
#import statistics as stat
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# 17 歳の日本人男性全体 (母集団) の標本として 10 個抽出
# その平均で母数の平均を推測する
h = np.array([165.2, 175.9, 161.7, 174.2, 172.1, 163.3, 170.9, 170.6, 168.4, 171.3])
h_ave = np.mean(h)
print(h_ave)
169.36

推定量と 推定値

母集団の値 (母数) の推定の為に利用した、標本統計量を「推定量」と言う。上記の場合なら「算術平均値」を母数の推定量に使った、という事。

また実際に標本データを使って計算した値を「推定値」と言う。上記で言えば h_ave = 169.36 が推定値。

プログラミング的に表現するなら推定量は理論的に算出出来る「関数」で、推定値は「リテラル」。

In [18]:
# numpy で不偏分散を計算
h_var = np.var(h, ddof=1)
print(h_var)

# デフォルトだと標本分散になる
print(np.var(h))
21.6671111111
19.5004
In [19]:
# 乱数生成 numpy.random で出来る
# https://qiita.com/yubais/items/bf9ce0a8fefdcc0b0c97
from numpy.random import *
# 1 から 6 の整数の乱数を 1 つ生成
print(randint(1,6+1))
list(randint(1,6+1,20))
4
Out[19]:
[2, 3, 4, 3, 1, 6, 3, 6, 6, 2, 6, 3, 2, 4, 1, 1, 5, 1, 5, 5]
In [20]:
from collections import Counter
counter = Counter(randint(1,6+1,6))
for word, cnt in counter.most_common():
    print(word,cnt)

print("\n")
counter100 = Counter(randint(1,6+1,1000000))
for word,cnt in counter100.most_common():
    print(word,cnt)
5 2
1 1
2 1
4 1
6 1


2 167127
4 166761
6 166677
3 166649
1 166457
5 166329
In [25]:
plt.bar(["M", "F"], [2/3, 1/3])
plt.show()
In [37]:
# 正規分布の確率密度関数
f = lambda x: (np.exp(-x**2/2)) / np.sqrt(2*np.pi)
n = 1000
p = [f(x) for x in np.linspace(-4,4,n)]
plt.scatter(np.linspace(-4,4,n),p)
plt.show()
In [45]:
# 正規母集団から単純無作為抽出
# 平均 50、標準偏差 10、標本数 5
sample = normal(50, 10, 5)
print(sample)
plt.hist(sample,bins=6)
plt.show()
[ 43.18817731  45.93021025  48.61258077  51.41517152  22.69788665]
In [47]:
# 標本数 1万、中心極限定理が観察できる
sample2 = normal(50, 10, 10000)
plt.hist(sample2, bins=16)
plt.show()
In [51]:
# 母集団から単純無作為抽出で標本を複数回抽出しても、当然抽出されるものは同じではない
# したがって、抽出する度にその標本値 (たとえば平均値) は異なる
sample3 = normal(50, 10, 10)
sample4 = normal(50, 10, 10)
print(np.mean(sample3), sample3)
print(np.mean(sample4), sample4)
45.0366051714 [ 54.10431473  47.91864812  55.47318835  46.74851373  51.75296109
  57.33983397  36.11324872  35.29190763  24.09140537  41.53203001]
61.4958695533 [ 62.52188233  40.08330706  53.17507861  53.59026655  62.85203939
  59.70500566  64.5804723   70.17621897  74.42216631  73.85225834]
In [62]:
# 標本値の分布を求める
# 平均 50、標準偏差 10の正規母集団から標本数 10 の単純無作為抽出を 1 万回行なう
# それぞれの回の平均値を計算して、その平均値を分布を求める
s_ave = np.empty(10000)
for i in range(10000):
    s = normal(50, 10, 10)  # N=10でサンプリング
    s_ave[i] = np.mean(s)   # 標本平均を計算
plt.hist(s_ave, bins=16)
plt.show()
In [70]:
s_5 = abs(s_ave - np.array([50]))<5
cnt = Counter(s_5)
for word, cnt in cnt.most_common():
    print(word, cnt)
print(np.mean(s_ave))
print(np.var(s_ave))
True 8858
False 1142
50.0017853539
10.1963531462
In [72]:
# 標本分散と普遍分散の標本分布
c = 10000
s_var = np.empty(c)
s_uvar = np.empty(c)

for i in range(c):
    s = normal(50, 10, 10)
    s_var[i] = np.var(s)
    s_uvar[i] = np.var(s, ddof=1)
print(np.mean(s_var), np.mean(s_uvar))
89.9455442826 99.9394936473
In [74]:
print(np.std(s_var), np.std(s_uvar))
42.3424359641 47.0471510713
In [85]:
plt.hist(s_var,range=(0,500),bins=50)
plt.show()
plt.hist(s_uvar,range=(0,500),bins=50)
plt.show()
In [93]:
s_var100 = s_var>200
cnt = Counter(s_var100)
for w, c in cnt.most_common():
    print(w,c)
s_uvar100 = s_uvar>200
cnt2 = Counter(s_uvar100)
for w, c in cnt2.most_common():
    print(w,c)
False 9799
True 201
False 9639
True 361
In [96]:
print(np.mean(np.sqrt(s_uvar)))
9.72558172541
In [108]:
# 中央値の標本分布
s_ave = np.empty(c)
s_med = np.empty(c)
for i in range(c):
    s = normal(50, 10, 10)
    s_ave[i] = np.mean(s)
    s_med[i] = np.median(s)
print(np.mean(s_ave), np.mean(s_med))
49.7934955374 49.8567878207
In [111]:
print(np.std(s_ave), np.std(s_med))
3.12068896852 3.7639414451
In [113]:
plt.hist(s_ave, bins=16)
plt.show()
plt.hist(s_med, bins=16)
plt.show()

2017年12月15日金曜日

「R によるやさしい統計学」を Python で実習する 第 3 章

引き続き第 3 章。jupyter notebookをhtmlにする方法が分かったので、それをコピペした。

BloggerにJupyter Notebookを貼り付ける方法

第3章 2つの変数の記述統計

相関(量的変数同士の関係)や連関(質的変数同士の関係)について

In [20]:
import statistics as stat
import numpy as np
import matplotlib.pyplot as plt

stat_test1 = np.array([6,10,6,10,5,3,5,9,3,3,11,6,1,9,7,5,8,7,7,9])
stat_test2 = np.array([10,13,8,15,8,6,9,10,7,3,18,14,18,11,12,5,7,12,7,7])
psych_test = np.array([13, 14, 7, 12, 10, 6, 8, 15, 4, 14, 9, 6, 10, 12, 5, 12, 8, 8, 12, 15])

# https://pythondatascience.plavox.info/matplotlib/%E6%95%A3%E5%B8%83%E5%9B%B3
plt.scatter(stat_test1, stat_test2)
# 他の図とスケールを合わせる為にx,y軸の最小最大値を指定しておく
plt.ylim(0, 20)
plt.xlim(0, 20)
plt.show()
In [16]:
plt.scatter(psych_test, stat_test1)
plt.ylim(0, 20)
plt.xlim(0, 20)
plt.show()
In [17]:
plt.scatter(psych_test, stat_test2)
plt.ylim(0, 20)
plt.xlim(0, 20)
plt.show()
In [72]:
# 共分散
# https://docs.python.jp/3/library/statistics.html
cov12 = np.mean((stat_test1 - np.mean(stat_test1)) * (stat_test2 - np.mean(stat_test2)))
print(cov12)

# 不偏共分散
uncov12 = np.sum((stat_test1 - np.mean(stat_test1)) * (stat_test2 - np.mean(stat_test2))) / (20 - 1)
print(uncov12)
# なんか書籍と計算結果が違う
3.55
3.73684210526
In [77]:
# 相関係数
# https://qiita.com/takaki@github/items/247ada674b594dd8fdce
cor12 = cov12 / (np.std(stat_test1) * np.std(stat_test2))
print(cor12)
cor12np = np.corrcoef(stat_test1, stat_test2)
print(cor12np[1,0])
0.333212171275
0.333212171275
In [90]:
# クロス集計
# http://docs.pyq.jp/python/pydata/pandas/pivot_table.html
import pandas as pd
df = pd.DataFrame({'math' : np.array(["h", "h","l","l","h","h","h","h","h","l","l","h","l","h","h","l","h","h","h","h"]),
                   'stat' : np.array(["l","l","l","l","h","h","h","h","h","h","l","l","l","h","l","h","h","h","h","h"])})
df.pivot_table(index='math', columns='stat', aggfunc=[len])
Out[90]:
len
stat h l
math
h 10 4
l 2 4
In [97]:
# ファイ係数
# 内包表記で変換
# http://www.lifewithpython.com/2014/09/python-list-comprehension-and-generator-expression-and-dict-comprehension.html
mathd = [0 if x == 'h' else 1 for x in df['math']]
statd = [0 if x == 'h' else 1 for x in df['stat']]
print(mathd,statd)
print(np.corrcoef(mathd,statd)[1,0])
[0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0] [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0]
0.35634832255