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

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

の写経をしたので記録を公開しておきます。

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

　 この記事では、19章をまとめておきます。

写経

　この章から関数と変数の名前をpep8に従うように変更してある。

以下使用モジュールとrcパラメータの設定。

import pylab
import scipy as sp
import scipy.stats
import random
import statistics as stat
import collections

pylab.rc('figure', dpi=200)

19.1 有意性のチェック

code 19.1

def plot_tdist(t_stat=-2.13165598142, num_bins=1000, num_samples=10000000, deg_of_freedom=198):
    t_dist = [sp.random.standard_t(deg_of_freedom) for i in range(num_samples)]
    pylab.hist(t_dist, bins=num_bins, weights=num_samples*[1.0/num_samples])
    pylab.axvline(t_stat, color='w')
    pylab.axvline(-t_stat, color='w')
    pylab.title(f'T-distribution with {deg_of_freedom} Degrees of Freedom')
    pylab.xlabel('T-statistic')
    pylab.ylabel('Probability')

plot_tdist()

code 19.2

メモ

　controlTimesとtreatmentTimesが急に出てきたけど、ファイルなどから読み込まれてはいないから使えない。なので、19.pyを参考にして生成した。 後で気が付いたが、code 19.3にcontrolTimesとtreatmentTimesを生成するコードが書いてあった。ふざけんな

random.seed(0)
treatment_times = [random.gauss(119.5, 5.0) for i in range(100)]
control_times = [random.gauss(120, 4.0) for i in range(100)]

def calc_tstat_pval():
    treatment_mean = stat.mean(treatment_times)
    control_mean = stat.mean(control_times)
    print(f'Treatment mean - control mean ={treatment_mean - control_mean} minutes')
    two_sample_test = sp.stats.ttest_ind(treatment_times, control_times, equal_var=False)
    print(f'The t-statistic from two-sample test is {two_sample_test[0]}')
    print(f'The p-value from two-sample test is {two_sample_test[1]}')

calc_tstat_pval()

Treatment mean - control mean =-0.9626159148482003 minutes
The t-statistic from two-sample test is -1.4677518436370114
The p-value from two-sample test is 0.1439256440649787

19.2 p値に注意せよ

メモ

　303ページに別の標本でp値を計算していたので、僕もやってみた。random.seed(1)を実行することで、別の乱数から生成された標本が得られるの。（別に、random.seedを使わなくても異なる乱数は得られるのだが、seedの理解の為にやった。）

random.seed(1)
treatment_times = [random.gauss(119.5, 5.0) for i in range(100)]
control_times = [random.gauss(120, 4.0) for i in range(100)]

def calc_tstat_pval():
    treatment_mean = stat.mean(treatment_times)
    control_mean = stat.mean(control_times)
    print(f'Treatment mean - control mean ={treatment_mean - control_mean} minutes')
    two_sample_test = sp.stats.ttest_ind(treatment_times, control_times, equal_var=False)
    print(f'The t-statistic from two-sample test is {two_sample_test[0]}')
    print(f'The p-value from two-sample test is {two_sample_test[1]}')

calc_tstat_pval()

Treatment mean - control mean =-1.4022836116031243 minutes
The t-statistic from two-sample test is -2.277359505985212
The p-value from two-sample test is 0.023866362341807106

figure 19.4

def plot_pvals(num_trial=10000, bins=100):
    pvals = []
    pval_small = 0
    for i in range(num_trial):
        treatment_times = [random.gauss(119.5, 5.0) for i in range(100)]
        control_times = [random.gauss(120, 4.0) for i in range(100)]
        pval = sp.stats.ttest_ind(treatment_times, control_times, equal_var=False)[1]
        pvals.append(pval)
        if pval < 0.05:
            pval_small += 1
    pval_median = stat.median(pvals)
    pval_mean = stat.mean(pvals)
    fraction_small = pval_small/num_trial
    pylab.hist(pvals, bins=bins, weights=num_trial*[1/num_trial])
    pylab.title('P-values from Actual Distribution')
    pylab.xlabel('p-value')
    pylab.ylabel('Probability')
    pylab.annotate(f'Mean = {round(pval_mean, 4)}, Median = {round(pval_median, 4)}\nFraction below 0.05 = {round(fraction_small, 4)}', xy=(0.1, 0.7), xycoords='axes fraction')

plot_pvals()

1標本検定追加

#random.seed(1)
treatment_times = [random.gauss(119.5, 5.0) for i in range(100)]
control_times = [random.gauss(120, 4.0) for i in range(100)]

def calc_tstat_pval():
    treatment_mean = stat.mean(treatment_times)
    control_mean = stat.mean(control_times)
    print(f'Treatment mean - control mean ={treatment_mean - control_mean} minutes')
    two_sample_test = sp.stats.ttest_ind(treatment_times, control_times, equal_var=False)
    print(f'The t-statistic from two-sample test is {two_sample_test[0]}')
    print(f'The p-value from two-sample test is {two_sample_test[1]}')
    one_sample_test = sp.stats.ttest_1samp(treatment_times, 120)
    print(f'The t-statistic from one-sample test is {one_sample_test[0]}')
    print(f'The p-value from one-sample test is {one_sample_test[1]}')

calc_tstat_pval()

Treatment mean - control mean =-0.3566556416159159 minutes
The t-statistic from two-sample test is -0.5549758199376972
The p-value from two-sample test is 0.5795487828148798
The t-statistic from one-sample test is -1.5022019326317575
The p-value from one-sample test is 0.13622803414860576

19.4 有意性の有無

def john1(num_games=1273, lyndsay_wins = 666):
    outcomes = lyndsay_wins*[1.0] + (num_games - lyndsay_wins)*[0.0]
    print(f'The p-value from a one-sample test is {sp.stats.ttest_1samp(outcomes, 0.5)[1]}')

john1()

The p-value from a one-sample test is 0.0982205871243577

code 19.4

def lyndsay(num_games=1273, lyndsay_wins=666, num_trials=1000):
    at_least = 0
    for t in range(num_trials):
        l_wins = 0
        for g in range(num_games):
            if random.random() < 0.5:
                l_wins += 1
        if l_wins >= lyndsay_wins:
            at_least += 1
    print(f'Probability of result at least this extreame by accident = {at_least/num_trials}')

lyndsay()

Probability of result at least this extreame by accident = 0.046

code 19.5

def john2(num_games=1273, lyndsay_wins=666, num_trials=1000):
    at_least = 0
    for t in range(num_trials):
        l_wins, j_wins = 0, 0
        for g in range(num_games):
            if random.random() < 0.5:
                l_wins += 1
            else:
                j_wins += 1
        if l_wins >= lyndsay_wins or j_wins >= lyndsay_wins:
            at_least += 1
    print(f'Probability of result at least this extreame by accident = {at_least/num_trials}')

john2()

Probability of result at least this extreame by accident = 0.088

figure 19.5

def plot_fig_19_5(mean_diff=0.5, stdev=5, num_trials=1000, start=1, stop=3000, step=50):
    mean1 = 0
    mean2 = mean1 + mean_diff
    mean_pvals=[]
    for i in range(start, stop, step):
        pvals = []
        for t in range(num_trials):
            pvals.append(sp.stats.ttest_ind([random.gauss(mean1, stdev) for j in range(i)],
                                            [random.gauss(mean2, stdev) for j in range(i)], equal_var=False)[1])
        mean_pvals.append(stat.mean(pvals))
    pylab.title(f'Gaussians with SD={stdev}, Means of {mean1} & {mean2}')
    pylab.plot(range(start, stop, step), mean_pvals, label='Mean p-value')
    pylab.xlabel('Sample size')
    pylab.ylabel('P-value')
    pylab.semilogy()
    pylab.axhline(0.05, linestyle='-', label='p = 0.05')
    pylab.axhline(0.01, linestyle='--', label='p = 0.01')
    pylab.legend(loc='best')

plot_fig_19_5()

19.6 多重仮説検定

bm_results2012_txt='bm_results2012.txt'

def get_BM_data(filename):
    data = {'name': [], 'gender': [], 'age':[], 'division': [], 'country': [], 'time': []}
    with open(filename) as f:
        for line in f.readlines():
            split = line.split(',')
            data['name'].append(split[0])
            data['gender'].append(split[1])
            data['age'].append(split[2])
            data['division'].append(split[3])
            data['country'].append(split[4])
            data['time'].append(float(split[5][:-1]))
    return data

data = get_BM_data(bm_results2012_txt)

code 19.6

countries_to_compare = ['BEL', 'BRA', 'FRA', 'JPN', 'ITA']
def fem_time_per_ctry(data, countries_to_compare):
    country_times = collections.defaultdict(list)
    for i in range(len(data['name'])):
        if data['country'][i] in countries_to_compare and data['gender'][i] == 'F':
            country_times[data['country'][i]].append(data['time'][i])
    return country_times

def comp_time(data, countries_to_compare):
    country_times = fem_time_per_ctry(data, countries_to_compare)
    for c1 in countries_to_compare:
        for c2 in countries_to_compare:
            if c1 < c2:
                pval = sp.stats.ttest_ind(country_times[c1], country_times[c2], equal_var=False)[1]
                if pval < 0.05:
                    print(f'{c1} and {c2} have significantly different means, p-value = {round(pval, 4)}')

comp_time(data, countries_to_compare)

ITA and JPN have significantly different means, p-value = 0.025

code 19.7

メモ

　このコードはどうやら、ガウス分布からサンプリングした値をpopulationリストに5000個保存して、populationリストからサンプリングしているようだ。どうせだったらガウス分布から直接サンプリングすればいいのにと思ったのでそうやるようなプログラムに書き直してみた。

　あと、テキストのsampleSize=30はsampleSize=200の間違いなのではないだろうかと思う。

def mult_hyp(num_hyps=20, sample_size=200):
    sample = collections.defaultdict(list)
    for i in range(num_hyps):
        sample[1].append([random.gauss(0, 1) for i in range(sample_size)])
        sample[2].append([random.gauss(0, 1) for i in range(sample_size)])
    num_sig = 0
    for i in range(num_hyps):
        if sp.stats.ttest_ind(sample[1][i], sample[2][i])[1] < 0.05:
            num_sig += 1
    print(f'Number of statistically significant (p < 0.05) result = {num_sig}')

mult_hyp()

Number of statistically significant (p < 0.05) result = 0

メモ

　mult_hyp関数の結果は1になったり0だったりと安定しないので、どんな感じの分布なのかを見てみるためのヒストグラムを書いてみた。

def mult_hyp(num_hyps=20, sample_size=200):
    sample = collections.defaultdict(list)
    for i in range(num_hyps):
        sample[1].append([random.gauss(0, 1) for i in range(sample_size)])
        sample[2].append([random.gauss(0, 1) for i in range(sample_size)])
    num_sig = 0
    for i in range(num_hyps):
        if sp.stats.ttest_ind(sample[1][i], sample[2][i])[1] < 0.05:
            num_sig += 1
    # print(f'Number of statistically significant (p < 0.05) result = {num_sig}')
    return num_sig

def plot_num_sig(num_trial=1000, num_hyps=20, sample_size=200):
    num_sigs = [mult_hyp(num_hyps, sample_size) for i in range(num_trial)]
    min_sig = min(num_sigs)
    max_sig = max(num_sigs) 
    bins = max_sig - min_sig + 1
    pylab.hist(num_sigs, range=(min_sig - 0.5, max_sig + 0.5), density=True, bins=bins)
    pylab.title('Probability of significants')
    pylab.xlabel('Significants')
    pylab.ylabel('Probability')

plot_num_sig()

　なるほど。稀に、20回のうちにp値が0.05を下回るとき6回あるときがあるのか。

　そもそも、正規分布から200個サンプリングした二組のt-検定量がどんな感じなんだ？と思って、最後にプロットしてみた。

def plot_pval(num_trials=100000, sample_size=200, sigma=1):
    t_stats = []
    for i in range(num_trials):
        sample = {}
        sample[1] = ([random.gauss(0, sigma) for i in range(sample_size)])
        sample[2] = ([random.gauss(0, sigma) for i in range(sample_size)])
        ttest = sp.stats.ttest_ind(sample[1], sample[2])
        t_stats.append(ttest[0])
    pylab.hist(t_stats, density=True, )
    pylab.title('Histogram of t-statistics')
    pylab.ylabel('probability')
    pylab.xlabel('the t-statistic')

plot_pval()

　どうやら、このt-検定の結果はだいたい正規分布に従うようだ。

おまけ

random.seedについて

　random.seedの意味が分からなかったので少し遊んでみた。 ドキュメントはこれrandom.seed(a=None, version=2)。 ここを読む限り、乱数の生成にはOSの時刻が使われているようだが、random.seedを実行することで、生成される乱数を制御できるようだ。実際にrandom.seedを実行するのとしないのを比較したらその違いがよく分かった。random.seed(0)を実行しないと各乱数の生成で得られるリストは異なるものになるが、実行すると毎回同じ乱数が得られる。

#random.seed(0)
[random.random() for i in range(4)]

[0.8159130965336595,
 0.1006075202160962,
 0.14635848891230385,
 0.6976706401912388]

#random.seed(0)
[random.random() for i in range(4)]

[0.5112747213686085,
 0.4049341374504143,
 0.7837985890347726,
 0.30331272607892745]

random.seed(0)
[random.random() for i in range(4)]

[0.8444218515250481,
 0.7579544029403025,
 0.420571580830845,
 0.25891675029296335]

random.seed(0)
[random.random() for i in range(4)]

[0.8444218515250481,
 0.7579544029403025,
 0.420571580830845,
 0.25891675029296335]

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