<jupyter_start><jupyter_text>サンプルの分布とデータファイルの読み方前の節ではデータは全て Python で乱数を生成していた．この節では実際のデータをファイルから読み込む方法を説明し，そのデータを使って確率の世界の基本法則である**大数の法則**と**中心極限定理**を説明する．  ファイルの読み込みこの実習ではサンプルデータを読み込んで処理することが多いので，その練習をしておこう．屋久島の降水量 (気象庁で公開されているデータを加工) を読み込んでヒストグラムを表示してみる．まずは標準関数を使って書くと，<jupyter_code>from matplotlib import pyplot
%matplotlib inline
import re

rain_data = open("./dat/yakushima-rain-year.txt", "r")

x = []
y = []
sum = 0

for line in rain_data:
    if re.match(r"^\s+$", line):
        continue
    tmp = line.split()
    x.append(int(tmp[0]))
    y.append(float(tmp[1]))
    sum = sum + float(tmp[1])

rain_data.close()

mean = sum / len(y)

pyplot.hist(y)
pyplot.show()

pyplot.hlines([mean], 1938, 2018, "red")
pyplot.plot(x,y)
pyplot.show()


print("Average = %.1f" % mean)<jupyter_output><empty_output><jupyter_text>これを numpy で書くと下のようになる．そろそろ numpy のありがたみが分かってきたのではないだろうか．ファイルのオープンやクローズなどせずに`loadtxt()` だけで読み込みができるし，`unpack`オプションで行と列の転置も勝手にやってくれる．ついでに平均を求める `mean()` を使って平均値もすぐに求まる．<jupyter_code>import numpy as np

x,y = np.loadtxt("./dat/yakushima-rain-year.txt", unpack=True)

pyplot.hist(y)
pyplot.show()

pyplot.hlines([y.mean()], 1938, 2018, "red")
pyplot.plot(x,y)
pyplot.show()

print("Average = %.1f" % y.mean())<jupyter_output><empty_output><jupyter_text>ところで東京都の平年降水量は1528.8 mm (1981年～2010年 の30年統計)である．屋久島は「1ヶ月35日は雨が降る」と言われるが，おおげさではないことが分かる．  2次元配列のデータ./dat/yakushima-rain-month.txt に，月別の降水量データがある．年ごとのデータと違って2次元配列になりそうだが，これのヒストグラムを描くコードを書いてみよう． ヒント- 1列目は年なので，ここはいらない．`loadtxt()` で使う列を指定するには`usecols`オプションを使う．読み込む列のリストで指定するので，`range()`が使える．- 欲しいヒストグラムを得るには `hist()` に1次元配列でデータを渡す必要がある．2次元(多次元)配列を1次元配列にするには `flatten()`メソッドを使う<jupyter_code>data = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))

pyplot.hist(data.flatten())
pyplot.show()

print("Average = %.1f" % data.mean())<jupyter_output><empty_output><jupyter_text>用語の整理ここから先は用語をきちんと理解しておかないと混乱するので，用語を整理しておく．JISの統計用語 (JIS Z8101-1，ISO 3534-1に対応) に従うと，以下のように書ける．- **母集団 (population)** — 検討の対象となるアイテムの全体．- **サンプリング単位** — 母集団を分割した最小単位．- **標本，試料 (sample)** — 1つ以上のサンプリング単位からなる母集団の部分集合．- **観測値，測定値 (observed value)** — サンプルのある要素に関して観測された値．標本は母集団の統計的性質を反映したものであり，どのようにサンプリングされたかが重要である．一方で観測値はどのように取得したかを問題にしない．つまり統計的性質を気にしないという違いがある．さしあたり標本と観測値の違いは置いておいて，母集団と標本だけ考えていくことにする．"Statistics Hacks" では **スコア** と書くことで標本と観測値の違いをうまくあいまいにしている (ような気がする)．**観測値の集合**に対して数値的要約 (平均，標準偏差など) を求めるのを**記述統計 (descriptive statistics)** と呼ぶ．つまり，既知のスコアの性質を調べることである．一方で，標本から母集団の性質を推測するのが**推測統計 (inferential staitstics)** と呼ぶ．  基本法則統計をやるときに「とりあえず標本を増やせばいい」という考えの根拠になる2つの基本法則を確かめよう．- **大数の法則** — 母集団から無作為抽出した標本の平均は，標本のサイズを大きくすると母集団の平均 (母平均) に収束する- **中心極限定理** — 標本平均と母平均の誤差は，標本のサイズを大きくしたとき**正規分布**に収束する混同していたり同じものだと思っている人が結構いるが，別の法則である．また，統計の大原則のように思われがちだが，大数の法則が成立しない分布 (コーシー分布) や中心極限定理が正規分布に収束しない場合がないわけではない．ただ，そういう難しい話は考えないことにしよう．そんなことより忘れがちで大事はことは，中心極限定理は**母集団の分布が正規分布かどうかに関わらず**成立するということだ．上で使った屋久島の月別降水量を使って確かめよう．もう一度ヒストグラムを出してみると，<jupyter_code>data = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))

pyplot.hist(data.flatten())
pyplot.show()

print("Average = %.1f" % data.mean())<jupyter_output><empty_output><jupyter_text>どうも正規分布っぽくない．(正規分布がどういうものかの詳細は後で．)まずは大数の法則を確かめよう．ランダムサンプリングの回数を増やしていって，標本平均がどう変化するかをプロットしてみる．`numpy` を使えばわざわざ乱数を用意してあれこれしなくても簡単にランダムサンプリングができる．<jupyter_code>#---- パラメータ ----
# N_max : サンプリング回数の上限
N_max = 100
#------------------

x = []
y = []

population = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))

for N in range(1, N_max):
    sample = np.random.choice(population.flatten(), N, replace=False)
    x.append(N)
    y.append(sample.mean())
    
pyplot.xlabel("#samples")
pyplot.ylabel("sample mean")
pyplot.plot(x,y)
pyplot.hlines([population.mean()], 0, N_max, "red")
pyplot.show()<jupyter_output><empty_output><jupyter_text>大数の法則によれば，標本を増やすと標本平均は母集団の平均 (赤い線)に収束するはずである．$N=100$ だと収束しきっていない感じだが，Nmax を増やすと収束の様子が見える．Nmax を増やして結果を観察しよう．ちなみに母集団の要素数は972なので Nmax の上限は972 で，言うまでもなく N=972 のとき標本平均は母集団平均に一致する．次に中心極限定理を確かめる．母集団から N 個の標本を無作為抽出して，その標本群の平均と母平均との差を求める，という操作を Ndist 回試行してヒストグラムを描くプログラムである．(言葉で書くとややこしい)<jupyter_code>#---- パラメータ ----
# N : 標本の抽出数
N = 2
# N_dist : 試行回数
N_dist = 2000
#------------------

val = []

population = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))
# ave : 母平均
ave = population.mean()

for i in range(1, N_dist):
    sample = np.random.choice(population.flatten(), N, replace=False)
    val.append(sample.mean() - ave)
    
pyplot.hist(val)
pyplot.show()<jupyter_output><empty_output><jupyter_text>$N=1$だと元の分布形状が見えるし，最頻値が0からずれているが，$N$を大きくしていくとだんだん分布形状が左右対象で最頻値も0になり，正規分布に近づいていくことが分かる．つまり，標本の数を増やすと標本平均は高い確率で母集団の平均に一致するし，母集団の平均からのずれは正規分布に従う．これが中心極限定理である．  正規分布の前に説明を飛ばしていた正規分布について詳しく見ていこう．教科書なんかには，> 正規分布 (またはガウス分布) とは確率密度関数$f(x)$が平均$\mu$と分散$\sigma^2$を使って> $$f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$> で記述される分布で，$N(\mu, \sigma^2)$と書く．という感じで書かれている．(Wikipedia にはそう書かれている．)では平均 $\mu$ と標準偏差 $\sigma$ (2乗した$\sigma^2$が分散)は何者なのか？  平均 (mean)「平均」にもいろいろある，という話は置いておいて，算術平均 (相加平均) を考える．算術平均は言うまでもなく$$\mu = \frac{1}{N}\sum_{i=1}^N x_i = \frac{x_1 + x_2 + \cdots + x_N}{N}$$だが，この計算の意味は何だろうか？平均値$\mu$は，「全ての要素からの距離が最小になる点」になっている．距離とは何か？ を考え始めるとそれはまた難しいので二乗偏差の総和を最小化する，と言いかえてもいい．下のコードは，ある値と母集団の点との二乗偏差の和の合計をプロットするものである．赤の線が平均値 (348.4) を示す．<jupyter_code>import math

population = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))

x = []
y = []
ave = population.mean()

for a in range(200,500):
    d = 0
    for s in population.flatten():
        d += abs(s - a)**2
        
    x.append(a)
    y.append(math.sqrt(d))

pyplot.plot(x,y)
pyplot.xlabel("Value")
pyplot.ylabel("Distance")
pyplot.vlines([ave], 7e3, 9e3, "red")
pyplot.show()<jupyter_output><empty_output><jupyter_text>用語： average と mean「平均」はよく average と書くがこれは「中心傾向を表す代表値」というやや漠然としたもので，厳密に平均値を表すときは mean と書く．例えば RMS (二乗平均平方根) は Root Mean Square だ．中心傾向を表す値としては他にも中央値 (median) や最頻値 (mode) もある．ただ，正規分布であれば mean も median も mode も等しくなるので，あまり気にしない．  標準偏差標準偏差は「平均値との距離の平均」に相当する．  余談：定義が違う？「標準偏差」の定義は本によって違うことがある．なお，python の `numpy.std()` は $n$，その他の言語やExcel はだいたい $n-1$で計算するのでちょっと注意が必要かもしれない．ちなみに同じ python の `pandas.std()` は $n-1$で計算する．`statistics` では $n$ で計算する `pstdev()` と $n-1$ で計算する `stdev()` が用意されている．(オプション ddof をつければ `numpy.std()` を $n-1$で計算させることも，`pandas.std()` を $n$ で計算させることもできる．)  正規分布  中心極限定理の再確認中心極限定理の確認で表示したヒストグラムに正規分布のグラフを重ねて，どのぐらい正規分布かをもう少し詳しく見られるようにしてみよう．他にうまいやり方があるのかも知れないが，確率密度関数 `scipy.stats.norm.pdf()` を使ってグラフを描くために，ヒストグラムのビンの幅を固定したりごちゃごちゃやっている．$N$ を大きくしていくと，ヒストグラムと正規分布のグラフが重なっていくことが確認できるだろう．<jupyter_code>from scipy.stats import norm

#---- パラメータ ----
# N : 標本数
N = 2
# N_dist : 試行回数
N_dist = 2000
#------------------

# bin_width : ヒストグラムのビンの幅
bin_width = 20

val = np.array([])

population = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))
ave = population.mean()

for i in range(1, N_dist+1):
    sample = np.random.choice(population.flatten(), N, replace=False)
    val = np.append(val, sample.mean() - ave)

x = np.arange(-5 * val.std(), 5 * val.std(), 1)
pyplot.hist(val, bins=np.arange(min(val), max(val) + bin_width, bin_width))
pyplot.plot(x, norm.pdf(x, val.mean(), val.std()) * N_dist * bin_width)
pyplot.show()<jupyter_output><empty_output><jupyter_text>どのぐらい正規分布か — 尖度，歪度上の例ではヒストグラムと正規分布がどのぐらい似ているか，主観的に評価した．定量的に「どのぐらい正規分布しているか」を評価したければ，いくつかの指標を使うことができる．平均値$\mu$と分散$\sigma$で正規化した標準化確率変数$Z$$$Z=\frac{x - \mu}{\sigma}$$に対して，3次のモーメント (期待値)を**歪度 (skewness)**，4次のモーメント (から3を引いたもの)を**尖度 (kurtosis)**と呼ぶ．$$歪度 = E\left[Z^4\right]=E\left[\frac{(x - \mu)^4}{\sigma^4}\right] = \frac{E\left[(x-\mu)^4\right]}{\sigma^4},$$$$尖度 = E\left[Z^3\right] - 3=E\left[\frac{(x - \mu)^3}{\sigma^3}\right] -3 = \frac{E\left[(x-\mu)^3\right]}{\sigma^3} - 3,$$ただし$E(x)$は期待値である．この2つが0に近いほど正規分布に近いと言える．(尖度の $-3$ は正規分布で値を0にするための補正)分布の形状の話でなぜ期待値が出てくるのか？ と思うかも知れないが，平均が重心に対応するように期待値はモーメントに相当するので，期待値は分布の偏りに関係しているのである．歪度は分布の偏りを示し，正の値であれば左側に，負の値であれば右側に偏っていることを表す．尖度は正であれば分布がより平均値付近に集中し，負であれば分布の裾野が広いことを表す (分散の大小とは異なることに注意)．ちなみに歪度も尖度も定義できないことがある．例えばコーシー分布のように，期待値が発散する分布では定義ができない． Python で尖度，歪度は `scipy.stats` の中の `skew()`，`kurtosis()` で求めることができる．サンプル数$N$を増やしていったときに尖度と歪度がどう変化するかを調べてみよう．<jupyter_code>import scipy

#---- パラメータ ----
# N_max : 標本数の上限
N_max = 20
# N_dist : 試行回数
N_dist = 2000
#------------------

val = []
x = []
ys = []
yk = []

population = np.loadtxt("./dat/yakushima-rain-month.txt", usecols=range(1,13))
ave = population.mean()

for s in range(1, N_max+1):
    val = []
    for i in range(1, N_dist+1):
        sample = np.random.choice(population.flatten(), s, replace=False)
        val.append(sample.mean() - ave)

    x.append(s)
    ys.append(scipy.stats.skew(val))
    yk.append(scipy.stats.kurtosis(val))
    
pyplot.plot(x, ys, label="skewness")
pyplot.plot(x, yk, label="kurtosis")
pyplot.hlines([0], 0, Nmax, "red")
pyplot.hlines([1], 0, Nmax, "orange", linestyle="dashed")
pyplot.hlines([2], 0, Nmax, "blue", linestyle="dashed")
pyplot.legend()
pyplot.xlabel("#sample")
pyplot.ylabel("skewness, kurtosis")
pyplot.show()<jupyter_output><empty_output><jupyter_text>順調に尖度，歪度ともに0に向かっているので，中心極限定理は正しそうだ．「**尖度・歪度がどの程度の値であれば正規分布と見なしてよいか**」という基準については諸説あるようで，$-1 < 歪度 < 1$，$-2 < 尖度 < 2$ がよく言われているようだ．この例で言えば$N=2$でもう正規分布と見なしてよいということになる．ヒストグラムがけっこう歪んでいても正規分布として取り扱ってもあまり問題はないらしい．正規分布かどうかの判定には**正規性の検定**というのがあるのだが，もう少し別の知識が必要なのでこれについては後でやることにする．  中心極限定理を使う中心極限定理を使うと，標本群がどういう母集団から抽出されたかを推測することができる．例えば屋久島の降水量が鹿児島県内の降水量から見てどのぐらい特殊か，を評価することができる． `./dat/kagoshima-rain-year.txt` は1978年～2018年の41年間に屋久島を含む鹿児島県内の26地点で記録された年間降水量 (気象庁提供) である．とりあえずヒストグラムを見てみよう．青が屋久島，オレンジが鹿児島全体，青線が屋久島の平均，赤線が鹿児島全体の平均である．予想通り，屋久島はかなり雨が多いようだ．<jupyter_code>import scipy

yakushima = np.loadtxt("./dat/yakushima-rain-year.txt", usecols=[1])
kagoshima = np.loadtxt("./dat/kagoshima-rain-year.txt", usecols=range(1, 27))

y_ave = yakushima.mean()
k_ave = kagoshima.mean()

pyplot.xlim([0, 7000])
pyplot.xlabel("annual precipitation")
pyplot.ylabel("population")
pyplot.hist(yakushima)
pyplot.vlines(y_ave, 0, 15, "blue")
pyplot.vlines(k_ave, 0, 15, "red")
pyplot.title("Yakushima")
pyplot.show()

print("#sample: %d, mean: %d, sigma: %.2f, skew: %.3f, kurtosis: %.3f" % (yakushima.size, yakushima.mean(), yakushima.std(), scipy.stats.skew(yakushima), scipy.stats.kurtosis(yakushima)))

pyplot.xlim([0, 7000])
pyplot.xlabel("annual precipitation")
pyplot.ylabel("population")
pyplot.hist(kagoshima.flatten(), color="orange")
pyplot.vlines(y_ave, 0, 400, "blue")
pyplot.vlines(k_ave, 0, 400, "red")
pyplot.title("Kagoshima")
pyplot.show()

print("#sample: %d, mean: %d, sigma: %.2f, skew: %.3f, kurtosis: %.3f" % (kagoshima.size, kagoshima.mean(), kagoshima.std(), scipy.stats.skew(kagoshima.flatten()), scipy.stats.kurtosis(kagoshima.flatten())))<jupyter_output><empty_output><jupyter_text>ただ，これはひょっとしたら標本の取り方でたまたま偏っているだけかも知れない．(80個も標本を取ってそんなにずれるはずがない，と思うだろうがとりあえずそう疑っておく．) 尖度と歪度は十分に小さいので，屋久島の標本は正規分布していると見てよさそうだ．(標本の)平均値は 4180，(標本の)標準偏差は 786 となっている．分布の細かな形状はさておき，標本平均 (標本値のばらつきではなく標本**平均**のばらつきという点に注意) がサンプリングによってどのぐらいばらつくかは$$\frac{\sigma}{\sqrt{N}}$$で計算でき，これを**標本平均の標準誤差**と呼ぶ．屋久島の分布で計算すると，$785.5/\sqrt{81} = 87.3$ となる．つまり，「本当は屋久島の年間降水量の平均は2614 mmなのに標本の取り方でたまたま4180 mmになってしまった」という確率は，平均値が$\sigma\times(4180 - 2614)/87.3 = 18\sigma$ ずれる確率に等しい，ということになる．ちなみに正規分布で標本が$\pm18\sigma$ずれる確率はだいたい$10^{143}$分の1である． 注1標本平均の標準誤差の定義はより正確には，要素数$N$の母集団から$n$個の標本を抽出し，その標本の標準偏差が$\sigma$のとき，$$\sqrt{\frac{N - n}{N}}\frac{\sigma}{\sqrt{n}}$$である．$N$が十分大きいとすると $\sigma/\sqrt{n}$ となる． 注2上の例は標本が無作為抽出でないとかいろいろ突っこみどころがあるので注意されたし．  大数の法則標本の性質は母集団の性質と必ずしも一致するわけではない．サイコロの目が等しい確率で出るからと言って，6回振れば1から6の目が必ず1回ずつ出るわけではないし，1%の確率で当たりがでるくじは100回引けば必ず当たるわけでもない．母集団の性質と，標本の性質の差 (の比率) は$$ \frac{1}{\sqrt{samples}}$$で表される．つまり標本数の平方根で差は減っていく．これが大数の法則である．  練習1%の確率で大当りの出るくじを模擬するコードを書き，- 「1回目の大当りが出るまでの試行回数」の分布- 標本内の大当り確率が1%に十分に近づくには試行回数は何回ぐらい必要かなどを調べよ．<jupyter_code>import random

Ns = 10000
x = []

for t in range(Ns):
    for i in range(10000):
        if (random.random() < 0.01):
            break
    x.append(i)

pyplot.hist(x, bins=100)
pyplot.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17139261744966444, 'fraction_numerical': 0.026677852348993288, 'mean_word_length': 7.318911374738311, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '251421', 'license': 'MIT', 'n_tokens_mistral': 8335, 'n_tokens_neox': 8175, 'n_words': 607, 'path': '02_Distribution.ipynb', 'repo_name': 'atuchiya/python-statistics-tutorial', 'size': '11920'}