26-8. 人文系の研究で卒業できなかった卒論のお話（後編）

import numpy as np
import numpy.random as nr
import scipy.stats as stats
import matplotlib.pyplot as plt

#乱数シード
nr.seed(0)

#パラメータ
num_student = 100
num_question = 500
p_border = 0.01


#高校１年時の様々なテストの大問スコア
ko1_score = (nr.randn(num_student, num_question) * 3 + 12).astype(int)
ko1_score = np.clip(ko1_score, 0, 20)
#確認
print(ko1_score.shape)
plt.hist(ko1_score.flatten(), bins=range(0, 21, 1))
plt.xlim(0, 20)
plt.show()

#入学大学（学部）の偏差値
hensachi = (nr.randn(num_student) * 9 + 50).astype(int)
hensachi = np.clip(hensachi, 25, 75)
#確認
print(hensachi.shape)
plt.hist(hensachi, bins=range(0, 101, 5))
plt.xlim(0, 100)
plt.show()

(100, 500)
(100,)

#大学偏差値が高い人のマスク、低い人のマスク
high_mask = (hensachi >= 60)
print(np.sum(high_mask))
low_mask = (hensachi <= 40)
print(np.sum(low_mask))

#大学偏差値が高い人の大問スコア、低い人の大問スコア
ko1_score_high = ko1_score[high_mask]
print(ko1_score_high.shape)
ko1_score_low = ko1_score[low_mask]
print(ko1_score_low.shape)

17
19
(17, 500)
(19, 500)

#ヒストグラム関数
def show(data1, data2, size=6):
    plt.hist(data1, bins=range(0, 21, 1), alpha=0.3, density=True, color='green', label='high hensachi daigaku')
    plt.hist(data2, bins=range(0, 21, 1), alpha=0.5, density=True, color='orange', label='low hensachi daigaku')
    
    #確率密度
    def gauss(x, mu, sigma):
        return np.exp(-(x - mu)**2/(2*sigma**2)) / (2*np.pi*sigma**2)**0.5
    
    #最尤推定でフィット曲線
    mu, sigma = np.mean(data1), np.std(data1)
    x_ = np.arange(0, 20, 0.001)
    y_ = gauss(x_, mu, sigma)
    plt.plot(x_, y_, color='green')
    
    mu, sigma = np.mean(data2), np.std(data2)
    x_ = np.arange(0, 20, 0.001)
    y_ = gauss(x_, mu, sigma)
    plt.plot(x_, y_, color='orange')
    
    plt.xlim(0, 20)
    plt.xlabel('ko1 daimon score')
    plt.legend()
    plt.show()

#全問題を調べる
for i in range(num_question):
    #対応のないt検定、等分散仮定
    t, p = stats.ttest_ind(ko1_score_high[:, i], ko1_score_low[:, i], equal_var=True)
    #有意差あり、highのほうが高平均なら表示
    if p < p_border and t > 0:
        print('i={}, p={}'.format(i, round(p, 4)))
        print('high mean={}, low mean={}'.format(round(np.mean(ko1_score_high[:, i]), 1), round(np.mean(ko1_score_low[:, i]), 1)))
        show(ko1_score_high[:, i], ko1_score_low[:, i])