「統計解析のはなし―データに語らせるテクニック」１～４章

大村平「統計解析のはなし―データに語らせるテクニック」日科技連　の読後行間補充メモ

同書籍は、印刷数表と手計算による統計解析を紹介した良書。

本稿では、同書籍の設例について、プログラムによる演算処理をした場合の例を示す。

コードと出力例は、Python［Google Colaboratory］

---

１章　前口上 

---

p3　偏差値

平均値、標準偏差、Ａ～Ｆ氏の素点　→　Ａ～Ｆ氏の偏差値　の計算例（付録１［p289］）

import pandas as pd

def devi_score(ave, std, raw): # 偏差値の関数
  return 50 + 10*(raw-ave)/std

ave = 55 # 平均
std = 12 # 標準偏差
raw_list = [31,43,55,67,73,79] # 素点

devi_list = []
for raw in raw_list:
  ds = devi_score(ave, std, raw)
  devi_list.append(ds)

DF = pd.DataFrame([raw_list,devi_list],
                  index=['素点','偏差値'],
                  columns = ['A','B','C','D','E','F'])
DF

 		A		B		C		D		E		F
素点 		31.0		43.0		55.0		67.0		73.0		79.0
偏差値	30.0		40.0		50.0		60.0		65.0		70.0

---

２章　素人探偵ものがたり（推定のはなし　その１） 

---

p23　平均値20万円、標準偏差10万円の正規分布

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

mean = 20 # 平均
sigma = 10 # 標準偏差
X_list = np.linspace(0, 59, 60) #0から59まで60個のリスト
Y1 = norm.pdf(X_list,mean,sigma) # 確率密度関数のデータ
Y2 = np.linspace(0,0,60) # 0から0まで60個のリスト（y=0 のデータ）

# グラフの描画
plt.plot(X_list, Y1) # 正規分布
plt.plot(X_list, Y2, color = 'black') # x軸
plt.vlines(0,0,0.05, color = 'black') # y軸
plt.vlines(20,0,0.05) # 平均の縦軸

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=10),fc= 'lightgrey')
plt.fill_between(X_list, Y1, Y2, where=(X_list>=30),fc= 'white')

plt.title(f'Normal Distribution, mean={mean}, sigma={sigma}')
plt.show()

p25　平均値の区間推定

• １個のサンプル値\(x_1\)
• 母集団の標準偏差\(\sigma\)
• 信頼区間内に母平均を含む割合\(r\)
• \(x_1,\sigma,r\)　→ 母平均\(\mu\)を区間推定

import pandas as pd
from scipy.stats import norm

sample_score = float(input('サンプル値',))
sigma = float(input('母集団の標準偏差',))
conf_intervals = [68.3, 90, 95]
DF = pd.DataFrame()
for conf_int_element in conf_intervals:
  accum_probability = 1-(50-conf_int_element/2)/100
  upper = round(sample_score + norm.ppf(accum_probability)*sigma,1)
  lower = round(sample_score - norm.ppf(accum_probability)*sigma,1)
  seed = pd.DataFrame([round(conf_int_element,1),lower,upper])
  DF = pd.concat([DF,seed], axis=1)
DF.columns = ['ex1','ex2','ex3']
DF.index = ['信頼区間','下限','上限']
DF.T

サンプル値20
母集団の標準偏差10
信頼区間	下限	上限
ex1	68.3	10.0	30.0
ex2	90.0	3.6	36.4
ex3	95.0	0.4	39.6

p26～29　標本が２つになれば(正規分布の加法性)

import numpy as np
import math
import pandas as pd

np.random.seed(0) # 乱数
ave = 60 # 平均
std = 10 # 標準偏差
num = 200000 # サンプル数
scores1 = pd.Series(np.around(np.random.normal(ave, std, num))) # 集団１
scores2 = pd.Series(np.around(np.random.normal(ave, std, num))) # 集団２

# 新たな集団の形成
scores3s = [] # 和の集団
scores4s = [] # 平均の集団
i = 0
while i < num:
  s1 = float(scores1.sample(replace=True)) # 集団１よりランダム抽出
  s2 = float(scores2.sample(replace=True)) # 集団２よりランダム抽出
  scores3s.append(s1+s2) # 抽出した両要素の和
  scores4s.append((s1+s2)/2) # 抽出した両要素の平均
  i += 1
scores3 = pd.Series(scores3s) # 和の集団リスト
scores4 = pd.Series(scores4s) # 平均の集団リスト

# DF化
DF = pd.DataFrame([scores1, scores2, scores3, scores4]).T
DF.columns = ['集団１','集団２','和の集団','平均の集団']
Means = pd.DataFrame([DF.mean()],index=['平均値']) # 各集団の平均
Stds = pd.DataFrame([DF.std()],index=['標準偏差']) # 各集団の標準偏差
DFL = pd.concat([DF, Means, Stds]) # DFの結合

# 結果表示
print(f'（平均値）\n設定値 {ave}\n２倍値 {ave*2}\n')
print(f'（標準偏差）\n設定値 {std}\n√2倍 {round(std*math.sqrt(2),2)}\n1/√2倍 {round(std/math.sqrt(2),2)}\n')
pd.options.display.precision = 2 # DF表示桁数
DFL

（平均値）
設定値 60
２倍値 120

（標準偏差）
設定値 10
√2倍 14.14
1/√2倍 7.07

p29　標本が２つになれば(正規分布の加法性)グラフ

# コードは上記の続き
from matplotlib import pyplot as plt

bins = np.linspace(0, 200, 101) # 0～200を101等分したリスト
class_value = (bins[:-1] + bins[1:]) / 2  # 階級値

def graph(score, name): # グラフ描画用
  freq = score.value_counts(bins=bins, sort=False) # 度数分布表
  dist = pd.DataFrame({"階級値": class_value
                      ,"度数": freq})
  plt.bar(dist['階級値'],dist['度数'], width=0.5, label = name)

graph(scores1, 'Population_1') # 母集団１のグラフ
graph(scores3, 'Population_3') # 母集団３のグラフ

plt.ylim(0, 25000)
plt.xlabel('Score')
plt.ylabel('Degree')
plt.legend()
plt.show()

p30　標本２つからの母集団の平均値の区間推定

• ２個のサンプル値からの標本平均\(\bar{x}\)
• 母集団の標準偏差\(\sigma\)
• 信頼区間内に母平均を含む割合\(r\)
• \(\bar{x},\sigma,r\)　→ 母平均\(\mu\)を区間推定

import numpy as np
from scipy.stats import norm
import pandas as pd

sample_score1 = float(input('サンプル値1=',))
sample_score2 = float(input('サンプル値2=',))
sample_average = np.average([sample_score1,sample_score2])
sample_number =2 #2個のサンプル
m_sigma = float(input('母集団の標準偏差=',))
sigma = m_sigma/np.sqrt(sample_number)
conf_intervals = [68.3, 90, 95]
DF = pd.DataFrame()
for conf_int_element in conf_intervals:
  accum_probability = 1-(50-conf_int_element/2)/100
  upper = round(sample_average + norm.ppf(accum_probability)*sigma,1)
  lower = round(sample_average - norm.ppf(accum_probability)*sigma,1)
  seed = pd.DataFrame([round(conf_int_element,1),lower,upper])
  DF = pd.concat([DF,seed], axis=1)
DF.columns = ['ex1','ex2','ex3']
DF.index = ['信頼区間','下限','上限']
Dataset = DF.T
Dataset

サンプル値1=10
サンプル値2=30
母集団の標準偏差=10
信頼区間	下限	上限
ex1	68.3	12.9	27.1
ex2	90.0	8.4	31.6
ex3	95.0	6.1	33.9

p42　標準正規分布

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

mean = 0 # 平均
sigma = 1 # 標準偏差
X_list = np.linspace(-3, 2.9, 60) #0から59まで60個のリスト
Y1 = norm.pdf(X_list,mean,sigma) # 確率密度関数のデータ
Y2 = np.linspace(0,0,60) # 0から0まで60個のリスト（y=0 のデータ）

# グラフの描画
plt.plot(X_list, Y1) # 正規分布
plt.plot(X_list, Y2, color = 'black') # x軸
plt.vlines(0,0,0.5, color = 'black') # y軸
plt.title('Standard Normal Distribution') # 表題

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=-sigma),fc= 'lightgrey')
plt.fill_between(X_list, Y1, Y2, where=(X_list>=sigma),fc= 'white')
plt.show()

p45　ｔ分布

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

X_list = np.linspace(-3, 3, 100) #-3から3で100個のリスト
f_inf = float('inf') # 無限大

for dof in [1,5,f_inf]: # 自由度のリスト
  Y_list = stats.t.pdf(X_list,dof) # t分布の確率密度関数のリスト
  plt.plot(X_list, Y_list, label=f"Φ={dof}") # グラフを描画

plt.hlines(0,-3,3, color = 'black') # x軸
plt.vlines(0,0,0.5, color = 'black') # y軸
plt.title('T Distribution') # 表題
plt.legend() # 凡例
plt.show()

p46　両裾の面積（確率）→ｔの値

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

dof = 3 # 自由度
prob = 0.1 # 両裾の面積（確率）
t_value = stats.t.ppf(prob/2, dof) # 片側確率のｔ値

X_list = np.linspace(-3, 3, 100) #-3から3で100個のリスト
Y1 = stats.t.pdf(X_list,dof) # t分布の確率密度関数のリスト
Y2 = np.linspace(0,0,100) # 0から0まで100個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list >= -t_value),fc= 'lightgrey')
plt.fill_between(X_list, Y1, Y2, where=(X_list <= t_value),fc= 'lightgrey')

plt.hlines(0,-3,3, color = 'black') # x軸
plt.vlines(0,0,0.4, color = 'black') # y軸
plt.title(f'TD, Both_prob={prob*100}% → t={-round(t_value,4)}') # 表題
plt.legend() # 凡例
plt.show()

p46　両裾の面積（確率）→ｔの値の一覧表

import numpy as np
from scipy import stats
import pandas as pd

f_inf = float('inf') # 無限大
horizon = [0.1,0.05,0.01] # 表の横軸(面積ｐ)
vertical = [1,2,3,f_inf] # 表の縦軸（自由度Φ）

t_list = []
for dof in vertical:
  tvalues = []
  for prob in horizon:
    t_value = round(-stats.t.ppf(prob/2, dof),3) # 片側確率のｔ値
    tvalues.append(t_value) # 各行を形成
  t_list.append(tvalues) # 各行を連結

# 階層付きDF化
mult_index = pd.MultiIndex.from_product([["Degrees of Freedom"], vertical]) #縦軸の階層化
mult_column = pd.MultiIndex.from_product([["Probability"], horizon])# 横軸の階層化
DF = pd.DataFrame(t_list, index =mult_index, columns= mult_column)
DF

---

３章　名探偵ものがたり（推定のはなし　その２） 

---

p59　χ二乗分布（自由度１）

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

dof = 1 # 自由度
prob = 0.1 # 両裾の面積（確率）
x_length = 8
plot_n =120

X_list = np.linspace(0, x_length, plot_n) #0から8で120個のリスト
Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで120個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画

plt.xlim(0, ) # 描画範囲を正の数に限定
plt.ylim(0, )

plt.title(r'$\chi^2$'  f'distribution') # 表題
plt.legend() # 凡例
plt.show()

p59　χ二乗分布（自由度１）の上側５％部分

# import 文は上記設例

dof = 1 # 自由度
prob = 0.1 # 両裾の面積（確率）
chi2_lower = stats.chi2.ppf(prob/2, dof) # 下側確率のχ二乗値
chi2_upper = stats.chi2.ppf(1-prob/2, dof) # 上側確率のχ二乗値
x_length = 6
plot_n =100

X_list = np.linspace(0, x_length, plot_n) #0から6で100個のリスト
Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで100個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画
plt.ylim(0, 0.5) # 描画範囲を限定

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=chi2_upper),fc= 'lightgrey')
plt.fill_between(X_list, Y1, Y2, where=(X_list<=chi2_lower),fc= 'lightgrey') # too small to display

plt.xlim(0, ) # 描画範囲を正の数に限定
plt.ylim(0, )

plt.title(r'$\chi^2$'  f', Both_P={prob*100} % →'  r'$\chi^2 upper$' f'={round(chi2_upper,3)}') # 表題
plt.legend() # 凡例
plt.show()

p59　χ二乗分布（自由度１）の下側５％部分

# dof, prob, chi2_lower, chi2_upper は上記設例数値を使用

x_length = 0.05
plot_n =1000

X_list = np.linspace(0, x_length, plot_n) #0から0.05で1000個のリスト
Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで1000個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画
plt.ylim(0, 20) # 描画範囲を限定

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=chi2_upper),fc= 'lightgrey') # outside the drawing area
plt.fill_between(X_list, Y1, Y2, where=(X_list<=chi2_lower),fc= 'lightgrey')

plt.xlim(0, ) # 描画範囲を正の数に限定
plt.ylim(0, )

plt.title(r'$\chi^2$'  f', Both_P={prob*100} % →'  r'$\chi^2 lower$' f'={round(chi2_lower,4)}') # 表題
plt.legend() # 凡例
plt.show()

p62　χ二乗分布（自由度３）の上側・下側５％部分

# import 文は上記設例

dof = 3 # 自由度
prob = 0.1 # 両裾の面積（確率）
chi2_lower = stats.chi2.ppf(prob/2, dof) # 下側確率のχ二乗値
chi2_upper = stats.chi2.ppf(1-prob/2, dof) # 上側確率のχ二乗値
x_length = 10
plot_n =100

X_list = np.linspace(0, x_length, plot_n) #0から10で100個のリスト
Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで100個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画
plt.ylim(0, max(Y1)*1.1) # 描画範囲を限定

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=chi2_upper),fc= 'grey')
plt.fill_between(X_list, Y1, Y2, where=(X_list<=chi2_lower),fc= 'grey')

plt.xlim(0, ) # 描画範囲を正の数に限定
plt.ylim(0, )

plt.title(f'BP={prob*100} % →'  r'$\chi^2$' f'=low_{round(chi2_lower,3)}, up_{round(chi2_upper,3)}') # 表題
plt.xlabel(r'$\chi^2$')
plt.legend() # 凡例
plt.show()

p65　χ二乗分布（自由度１～６）

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

x_length = 12
plot_n =120

for dof in [1,2,3,4,5,6]:
  X_list = np.linspace(0, x_length, plot_n) #0から12で120個のリスト
  Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
  Y2 = np.linspace(0,0,plot_n) # 0から0まで120個のリスト（y=0 のデータ）
  plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画

plt.xlim(0, ) # 描画範囲を正の数に限定
plt.ylim(0, )

plt.title(r'$\chi^2$'  f'distribution') # 表題
plt.legend() # 凡例
plt.show()

p65　χ二乗値の一覧表（自由度１～６）

import numpy as np
from scipy import stats
import pandas as pd

horizon = [0.95,0.05,0.025] # 表の横軸(面積ｐ)上側確率
vertical = [1,2,3,4,5,6] # 表の縦軸（自由度Φ）

chi_list = []
for dof in vertical:
  chivalues = []
  for prob in horizon:
    chi_value = round(stats.chi2.ppf(1-prob, dof),3) # 上側確率のχ二乗値
    chivalues.append(chi_value) # 各行を形成
  chi_list.append(chivalues) # 各行を連結

# 階層付きDF化
mult_index = pd.MultiIndex.from_product([["Degrees of Freedom"], vertical]) #縦軸の階層化
mult_column = pd.MultiIndex.from_product([["Probability"], horizon])# 横軸の階層化
DF = pd.DataFrame(chi_list, index =mult_index, columns= mult_column)
DF

      Probability                
            0.950   0.050   0.025
DOF 1       0.004   3.841   5.024
    2       0.103   5.991   7.378
    3       0.352   7.815   9.348
    4       0.711   9.488  11.143
    5       1.145  11.070  12.833
    6       1.635  12.592  14.449

p73　F分布のグラフ（下側５％、上側５％）

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

dfn = 9 # 分子の自由度
dfd = 6 # 分母の自由度

prob = 0.1 # 両裾の面積（確率）
F_lower = stats.f.ppf(prob/2, dfn, dfd) # 下側確率のF値
F_upper = stats.f.ppf(1-prob/2, dfn, dfd) # 上側確率のF値
x_length = 6
plot_n =100

X_list = np.linspace(0, x_length, plot_n) #0から x_lengthまで plot_n個のリスト
Y1 = stats.f.pdf(X_list,dfn, dfd) # F分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで plot_n個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"dfn={dfn} \ndfd={dfd}") # グラフを描画

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=F_upper),fc= 'lightgrey')
plt.fill_between(X_list, Y1, Y2, where=(X_list<=F_lower),fc= 'lightgrey')

plt.xlim(0, ) # 描画範囲を正の数に限定
plt.ylim(0, )

plt.title(f'F_Dist, Both_P={prob*100} % → F_low ={round(F_lower,2)}, F_up ={round(F_upper,2)}') # 表題
plt.xlabel('F')
plt.legend() # 凡例
plt.show()

p73　F値の一覧表（上側５％）

import numpy as np
from scipy import stats
import pandas as pd

prob = 0.1
dfns = np.linspace(1,10,10)  # 横軸(分子の自由度)
dfds = np.linspace(1,10,10)  # 縦軸(分母の自由度)

F_list = []
for dfd in dfds:
  Fvalues = []
  for dfn in dfns:
    F_value = stats.f.ppf(1-prob/2,dfn, dfd) # 上側確率のF値
    Fvalues.append(round(F_value,2)) # 各行を形成
  F_list.append(Fvalues) # 各行を連結

# 階層付きDF化
mult_index = pd.MultiIndex.from_product([["DoD"], dfds]) #縦軸の階層化
mult_column = pd.MultiIndex.from_product([["Degrees of Numerator"], dfns])# 横軸の階層化
DF = pd.DataFrame(F_list, index =mult_index, columns= mult_column)
DF

• 分子の自由度10、分母の自由度9について書籍73頁の表では「2.14」とあるが、正しくは「3.14」である。

p75～78　時計２つの遅れ進みデータ　→ばらつき比を区間推定

import math
from scipy.stats import f
import statistics

tick1 = [69,70,68,72,67,73,71]
tick2 = [57,45,49,50,44,55,43,51,55,51]


dfn = len(tick2)-1
dfd = len(tick1)-1
Nper = float(input('信頼区間[%]',))
pvalue = (100-Nper)/2/100
upper = f.ppf(1-pvalue, dfn, dfd)
lower = f.ppf(pvalue, dfn, dfd)
V1 = statistics.variance(tick1)
V2 = statistics.variance(tick2)
print(Nper,'%信頼区間：',round(1/math.sqrt(upper/(V2/V1)),3),'～',round(1/math.sqrt(lower/(V2/V1)),2))

信頼区間[%]95
95.0 %信頼区間： 0.956 ～ 4.67

p81　二項分布を正規分布で代用する

import numpy as np
import math
from scipy.stats import norm
from scipy.stats import binom
import matplotlib.pyplot as plt

# 二項分布
n_attempts = 100  #試行回数
prob = 0.2 #成功確率

x_list = range(n_attempts)    # 確率変数 x のリスト
prob_d = binom(n_attempts, prob) # 確率分布
y_list = prob_d.pmf(x_list)    # y のリスト
plt.bar(x_list, y_list, color = 'grey', label = 'Binomial D')

# 正規分布
mean = n_attempts * prob # 平均
sigma = math.sqrt(n_attempts * prob *(1-prob)) # 標準偏差
number = n_attempts

X_list = np.linspace(0, number-1, number) #0から number-1 まで numbre個のリスト
Y1 = norm.pdf(X_list,mean,sigma) # 確率密度関数のデータ
Y2 = np.linspace(0,0,number) # 0から0まで number個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, color = 'red', label = 'Normal D')

# 表示設定
plt.xlim(0,50) # 描画範囲を限定
plt.ylim(0,) 

plt.title(f'Binomial/Normal Distribution, attempts={n_attempts}, prob ={prob}')
plt.legend() # 凡例
plt.show()

p79～83　標本の割合　→母集団の割合を区間推定

import math
from scipy.stats import norm

n = int(input('標本となる学生人数：',))
r = int(input('標本中の横浜ファンの人数：',))

p_h = r/n
mu_h = n*p_h
sigma_h =math.sqrt(n*p_h*(1-p_h))

print('ｐ推定値：',p_h)
print('μ推定値：',mu_h)
print('σ推定値：',round(sigma_h,2))

# 割合pの区間推定
confidence_interval = float(input('信頼区間[％]：',))
accum_probability = 1-(50-confidence_interval/2)/100
upper = mu_h + norm.ppf(accum_probability)*sigma_h
lower = mu_h - norm.ppf(accum_probability)*sigma_h
print(f'割合ｐの区間推定：{round(lower,2)}～{round(upper,2)}％')

標本となる学生人数：100
標本中の横浜ファンの人数：30
ｐ推定値： 0.3
μ推定値： 30.0
σ推定値： 4.58
信頼区間[％]：95
割合ｐの区間推定：21.02～38.98％

p87～90　標本の割合　→母集団の割合をF分布で区間推定

from scipy.stats import f

n = int(input('標本数：',))
r = int(input('標本中の対象数：',))
phi1 = 2*(r+1)
phi2 = 2*(n-r)
phi1d = 2*(n-r+1)
phi2d = 2*r

confidence_interval = float(input('信頼区間[％]：',))
alpha_half = (100-confidence_interval)/2/100

F_a = f.ppf(1-alpha_half, phi1, phi2)
F_b = f.ppf(1-alpha_half, phi1d, phi2d)

P_up = phi1*F_a/(phi2+phi1*F_a)
P_low = phi2d/(phi2d+phi1d*F_b)
print('母集団中の割合ｐの',confidence_interval,'％信頼区間:　',round(P_low*100,2),'～',round(P_up*100,2),'％')

標本数：10
標本中の対象数：2
信頼区間[％]：95
母集団中の割合ｐの 95.0 ％信頼区間:　 2.52 ～ 55.61 ％

p91　確率紙

import numpy as np
from scipy.stats import f
import matplotlib.pyplot as plt

# 標本数 sn、信頼区間 ci
sn = int(input('標本数：',))
ci = float(input('信頼区間 [％]：',))

# 関数：　標本数 sn、標本比率 sp、信頼区間 ci　→母集団比率の信頼区間(Ｆ分布)
def CI_by_F(sn,sp,ci):
  stn = sn*sp # 標本中の対象数
  phi1 = 2*(stn+1)
  phi2 = 2*(sn-stn)
  phi1d = 2*(sn-stn+1)
  phi2d = 2*stn
  alpha_half = (100-ci)/2/100 # 信頼区間ci →F分布確率

  F_a = f.ppf(1-alpha_half, phi1, phi2)
  F_b = f.ppf(1-alpha_half, phi1d, phi2d)
  P_up = phi1*F_a/(phi2+phi1*F_a) # 上限
  P_low = phi2d/(phi2d+phi1d*F_b) # 下限
  return P_low, P_up

# 確率紙のグラフ
SP = np.arange(0, 1, 0.001) # 標本中の比率の列
MPU = np.array([]) # 母集団中の比率上限の列
MPL = np.array([]) # 母集団中の比率下限の列

for sp in SP:
  mpl = CI_by_F(sn,sp,ci)[0] # 下限値を関数計算
  MPL = np.append(MPL, mpl)
  mpu = CI_by_F(sn,sp,ci)[1] # 上限値を関数計算
  MPU = np.append(MPU, mpu)

plt.plot(SP, MPU, color='red')
plt.plot(SP, MPL, color='blue')
plt.grid(color='gray')
plt.xlabel('Sample proportion')
plt.ylabel('Population proportion')

標本数：10
信頼区間 [％]：90

---

４章　名行司ものがたり（検定のはなし） 

---

p102～105　観客の反応と性別差

import math

# 関数：　# 男性ｘ人のうちｍ人、女性ｙ人のうちｎ人が失神する確率
def faint(x,m,y,n):
   return math.comb(x,m)*math.comb(y,n)/math.comb(x+y, m+n)

# 女性が６人含まれる確率
x,m,y,n = [8,2,7,6] 
try1 = faint(x,m,y,n)

# 女性が７人含まれる確率
x,m,y,n = [8,1,7,7]
try2 = faint(x,m,y,n)

# 合計の確率
prob_faint = try1+ try2
f'{round(try1*100,2)}％＋{round(try2*100,2)}％＝{round(prob_faint*100,2)}％'

3.05％＋0.12％＝3.17％

p106～107　観客の反応と性別差（多人数版）

# 男性454人女性584人のうち失神者65人。同失神者の中に女性が44人以上含まれる確率
import math
import pandas as pd

allocation =[454,21,584,44]

# 関数：　# 男性ｘ人のうちｍ人、女性ｙ人のうちｎ人が失神する確率
def faint(x,m,y,n):
   return math.comb(x,m)*math.comb(y,n)/math.comb(x+y, m+n)

addmax = allocation[2]-allocation[3] # 女性総数が女性失神者の増加上限
decmax = allocation[1]-0 # ０名が男性失神者の減少下限
limit = min(decmax,addmax) # 失神者の増減移動の限界

DF = pd.DataFrame(allocation)

for m in range(limit):　# 女性44人から65人の全場合をDF化
  allocation[1] = allocation[1]-1
  allocation[3] = allocation[3]+1
  new_allo = pd.DataFrame(allocation)
  DF = pd.concat([DF,new_allo],axis=1)
DF.index = ['Male','M_faint','Female','F_faint']
DF.columns=range(DF.shape[1])
print(DF)

# 各場合の確率
Probability = []
for col in range(DF.shape[1]): # 列数だけ繰り返す
  prob = faint(DF[col]['Male'],DF[col]['M_faint'],DF[col]['Female'],DF[col]['F_faint'])
  Probability.append(prob)
Prob = pd.DataFrame(Probability).T
DF = pd.concat([DF, Prob])
DF.index = ['Male','M_faint','Female','F_faint','Probability']
sum = DF.loc['Probability'].sum()
print(f'確率：{round(sum*100,2)}％')

p109～114　食い違いの検定

import pandas as pd
import scipy.stats

# 基礎データ
DF = pd.DataFrame([
    [200,8],
    [150,22],
    [100,15]
    ],
                  index = ['Male', 'Female', 'Child'],
                  columns = ['headcount','value'])

headcount_sum = DF['headcount'].sum() # 縦合計値
value_sum = DF['value'].sum() # 縦合計値

# DF の形成
DF['expected'] = DF['headcount'] * value_sum / headcount_sum # 期待値
DF['difference'] = DF['value'] - DF['expected'] # 偏差
DF['square'] = DF['difference']**2 # 偏差の累乗
DF['Chi_square'] = DF['square']/DF['expected'] # 期待値で割る

# # χ二乗値、自由度　→確率
Chi_2 = DF['Chi_square'].sum() # χ二乗値
dof = DF.shape[0]-1 # 自由度
prob = 1 - stats.chi2.cdf(Chi_2, dof) # 確率

# 参考（確率、自由度　→χ二乗値）
prob_005 = 0.05 # 確率５％
chi_005 = stats.chi2.ppf(1-prob_005, dof) # 確率５％の χ二乗値

print(f'本件：χ二乗値:{round(Chi_2,2)} の確率:{round(prob*100,3)} ％')
print(f'参考：χ二乗値:{round(chi_005,2)} の確率:{prob_005*100}％')
DF

p111　適合度検定の関数（chisquare）を利用

import pandas as pd
from scipy.stats import chisquare

DF = pd.DataFrame([
    [200,8],
    [150,22],
    [100,15]
    ],
                  index = ['Male', 'Female', 'Child'],
                  columns = ['headcount','value'])

headcount_sum = DF['headcount'].sum() # 縦合計値
value_sum = DF['value'].sum() # 縦合計値
DF['expected'] = DF['headcount'] * value_sum / headcount_sum # 期待値
sta, pvalue = chisquare( DF['value'], DF['expected']) # 観察データ、期待値データ
print(f'確率：{round(pvalue*100,3)}％')
DF

p114　表4.2設例

import pandas as pd
import math
from scipy import stats

# ----------------------------------------- 設定値の入力
m_h = int(input('男性：',))
m_f = int(input('男性の失神者：',))
f_h = int(input('女性：',))
f_f = int(input('女性の失神者：',))

# ----------------------------------------- χ二乗検定
DFX = pd.DataFrame([
    [m_h,m_f],
    [f_h,f_f]
    ],
                  index = ['Male', 'Female'],
                  columns = ['headcount','value'])

headcount_sum = DFX['headcount'].sum()
value_sum = DFX['value'].sum()
DFX['expected'] = round(DFX['headcount'] * value_sum / headcount_sum,1) # 有効数字３桁
DFX['difference'] = DFX['value'] - DFX['expected']
DFX['square'] = round(DFX['difference']**2,1) # 有効数字３桁
DFX['Chi_square'] = round(DFX['square']/DFX['expected'],2) # 有効数字３桁

# # χ二乗値、自由度　→確率
Chi_2 = round(DFX['Chi_square'].sum(),2) # 有効数字３桁
dof = DFX.shape[0]-1 # 自由度
kprob = 1-stats.chi2.cdf(Chi_2, dof) # 確率

# 参考（確率、自由度　→χ二乗値）
prob_005 = 0.05 # 確率５％
chi_005 = stats.chi2.ppf(1-prob_005, dof) # 確率５％の χ二乗値

# ----------------------------------------- 出力
print(f'本件：χ二乗値：{Chi_2}の確率：{round(kprob*100,2)}％')
print(f'参考：χ二乗値：{round(chi_005,2)}の確率：{prob_005*100}％')
DFX

p116　設例（コードは上述 p114 と同じ）

p116　適合度検定の関数（chisquare）を利用

DFXの［value］[expected]列は上述 p114 と同じ

sta, pvalue = chisquare( DFX['value'], DFX['expected']) # 観察データ、期待値データ
print(f'確率：{round(pvalue*100,2)}％')

確率：10.29％

p118　数値０～９の乱数表　の偏り

import numpy as np
import pandas as pd
import scipy.stats as stats
from scipy.stats import chisquare

# 基礎データ
tenN  = np.linspace(0,9,10)
values = pd.Series([9,13,12,8,10,8,10,11,8,11]).T

# DFの形成
DF = pd.DataFrame([tenN,values]).T
DF.columns = ['Number','value']
Number_sum = DF['Number'].sum() # 縦合計値
value_sum = DF['value'].sum() # 縦合計値
DF['expected'] = value_sum / len(tenN) # 期待値
DF['difference'] = DF['value'] - DF['expected'] # 偏差
DF['square'] = DF['difference']**2 # 偏差の累乗
DF['Chi_square'] = DF['square']/DF['expected'] # 期待値で割る

# χ二乗値、自由度　→確率
Chi_2 = DF['Chi_square'].sum() # χ二乗値
dof = DF.shape[0]-1 # 自由度
prob = 1 - stats.chi2.cdf(Chi_2, dof) # 確率

# chisquare で検算：
sta, pvalue = chisquare( DF['value'], DF['expected']) # 観察データ、期待値データ

# 参考（確率、自由度　→χ二乗値）
prob_005 = 0.05 # 確率５％
chi_005 = stats.chi2.ppf(1-prob_005, dof) # 確率５％の χ二乗値

# 出力
print(f'本件：χ二乗値:{round(Chi_2,2)} の確率:{round(prob*100,3)} ％')
print(f'本件：chisquare での確率：{round(pvalue*100,3)}％')
print(f'参考：χ二乗値:{round(chi_005,2)} の確率:{prob_005*100}％')
DF

p120　両側５％危険率の検定（数値０～９の乱数表）グラフ

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

dof = 9 # 自由度
prob = 0.05 # 両裾の面積（確率）
magni = 1.3 # グラフ対描画平面の拡大率

chi2_lower = stats.chi2.ppf(prob/2, dof) # 下側確率のχ二乗値
chi2_upper = stats.chi2.ppf(1-prob/2, dof) # 上側確率のχ二乗値
x_length = chi2_upper*magni
plot_n =100

X_list = np.linspace(0, x_length, plot_n) #0から10で100個のリスト
Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで100個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=chi2_upper),fc= 'grey')
plt.fill_between(X_list, Y1, Y2, where=(X_list<=chi2_lower),fc= 'grey')

plt.xlim(0, chi2_upper*magni) # 描画範囲を正の数に限定
plt.ylim(0, max(Y1)*magni) # 描画範囲を限定

plt.title(f'Both Prob ={prob*100} % →  '  r'$\chi^2$' f'=low_{round(chi2_lower,3)}, up_{round(chi2_upper,2)}') # 表題
plt.xlabel(r'$\chi^2$')
plt.legend() # 凡例
plt.show()

p120　片側５％危険率の検定（数値０～９の乱数表）グラフ

# import 文は、上記と同様
# chi2_lower, chi2_upperは、上記設例（両側検定の上側確率）を使用

dof = 9 # 自由度
oneside_prob = 0.05 # 片側の面積（確率）
magni = 1.3 # 縦拡大率

chi2_lower_s = stats.chi2.ppf(oneside_prob, dof) # 下側確率のχ二乗値
chi2_upper_s = stats.chi2.ppf(1-oneside_prob, dof) # 上側確率のχ二乗値
x_length = chi2_upper*magni # 描画の横範囲を両側検定にあわせる
plot_n =100

X_list = np.linspace(0, x_length, plot_n) #0から10で100個のリスト
Y1 = stats.chi2.pdf(X_list,dof) # χ二乗分布の確率密度関数のリスト
Y2 = np.linspace(0,0,plot_n) # 0から0まで100個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画

# 塗りつぶし
plt.fill_between(X_list, Y1, Y2, where=(X_list>=chi2_upper_s),fc= 'blue')
plt.fill_between(X_list, Y1, Y2, where=(X_list<=chi2_lower_s),fc= 'orange')
plt.fill_between(X_list, Y1, Y2, where=(X_list>=chi2_upper),fc= 'grey') # 両側検定の描画
plt.fill_between(X_list, Y1, Y2, where=(X_list<=chi2_lower),fc= 'grey') # 両側検定の描画

plt.xlim(0, chi2_upper*magni) # 描画の横範囲を両側検定にあわせる
plt.ylim(0, max(Y1)*magni) # 描画の縦範囲を限定

plt.title(f'One_Sided_Prob ={oneside_prob*100} % →  '  r'$\chi^2$' f'=low_{round(chi2_lower_s,2)}, up_{round(chi2_upper_s,2)}') # 表題
plt.xlabel(r'$\chi^2$')
plt.legend() # 凡例
plt.show()

p120　平均値を検定する

gram = pd.Series([94,97,98,99,102,98]) # 標本値
m_average = 100 # 母平均（仮説）
side_prob = 0.05 # 片側確率（確率）

s_average = gram.mean() # 標本平均
s_s = gram.std(ddof=0) # 標本標準偏差
numbers = len(gram) # 要素数
dof = numbers -1 # 自由度
t_value = -stats.t.ppf(side_prob, dof) # 片側確率のｔ値
myu_up = s_average + t_value * (s_s/math.sqrt(numbers-1)) # ｔ値対応のｘ上限

def t_func(sa,ma,ss,num): # ｔ値を求める関数
  return (sa - ma)/(ss/math.sqrt(num-1))

tva = -t_func(s_average,m_average,s_s,numbers) # ｔ値
prob = stats.t.sf(tva, df = numbers-1) # ｔ値に対応する上側確率

# ｔ検定の結果表示
if prob>= 0.05:
  print(f't ： {round(tva,2)} ＜ {round(t_value,3)}')
else:
    print(f't ： {round(tva,2)} ＞ {round(t_value,3)}')

if prob>= 0.05:
  print(f'weight ： {round(m_average)} ＜ {round(myu_up,3)}')
else:
  print(f'weight ： {round(m_average)} ＞ {round(myu_up,3)}')

if prob>= 0.05:
  print(f'side_prob：{round(prob*100,1)}％ ＞ {side_prob*100}％,\n then the hypothesis [μ＝{m_average}g] is Not rejected.')
else:
  print(f'side_prob：{round(prob*100,1)}％ ＜ {side_prob*100}％,\n then the hypothesis [μ＝{m_average}g] is Rejected.')

# グラフ表記
X_list = np.linspace(95, 101, 100) #95から101で100個のリスト
Y1 = stats.t.pdf(X_list,dof,s_average) # t分布の確率密度関数のリスト
Y2 = np.linspace(0,0,100) # 0から0まで100個のリスト（y=0 のデータ）
plt.plot(X_list, Y1, label=f"Φ={dof}") # グラフを描画
plt.fill_between(X_list, Y1, Y2, where=(X_list>=myu_up),fc= 'lightgrey') # 塗りつぶし
plt.xlim(95,101) # 描画領域設定
plt.ylim(0,) # 描画領域設定
plt.vlines(s_average,0,stats.t.pdf(0,dof), color = 'blue') # 標本平均の場所に縦線
plt.vlines(m_average,0,stats.t.pdf(m_average-s_average,dof), color = 'red') # 母平均（仮説値）の場所に縦線
plt.title(f'TD, Upper_Side_Prob={side_prob*100}% → weight_upper={round(myu_up,3)}') # 表題
plt.legend() # 凡例
plt.show()

t ： 1.88 ＜ 2.015
weight ： 100 ＜ 100.145
side_prob：6.0％ ＞ 5.0％,
 then the hypothesis [μ＝100g] is Not rejected.
 
 

p125　ばらつきを検定する

import numpy as np
import statistics as stat
from scipy import stats

savings = np.array([10,10,30,30]) # 標本
sigma = 30 # 母標準偏差（仮定）
prob = 0.05 # 危険率

dof = len(savings)-1 # 自由度
ss = stat.pstdev(savings) # 標本の標準偏差
smean = stat.mean(savings) # 標本の平均
sum_squ = sum((s - smean)**2 for s in savings) # 標本の平方和
chi2_h = sum_squ/sigma**2 # 帰無仮説下でのχ二乗値
chi2 = stats.chi2.ppf(prob, dof) # 下側確率のχ二乗値

print(f'標本{savings}の標準偏差は{ss}。\n帰無仮説（母標準偏差{sigma}）を危険率{prob*100}％で検定。')
print(f'帰無仮説下でのχ二乗値：{round(chi2_h,2)}\n下側確率{prob*100}％のχ二乗値：{round(chi2,2)}')

if chi2_h > chi2:
  print('帰無仮説は維持。')
else:
  print('帰無仮説を棄却。')

標本[10 10 30 30]の標準偏差は10.0。
帰無仮説（母標準偏差30）を危険率5.0％で検定。
帰無仮説下でのχ二乗値：0.44
下側確率5.0％のχ二乗値：0.35
帰無仮説は維持。

p127　ばらつきの比率を検定する

import pandas as pd

tick1 = [69,70,68,72,67,73,71]
tick2 = [57,45,49,50,44,55,43,51,55,51]

prob = 0.05 # 危険率

def statistic(list):
  length = len(list) # 標本数
  smean = stat.mean(list) # 標本の平均
  ss = round(stat.pstdev(list),2) # 標本の標準偏差
  ubv = round(length/(length-1) * ss**2,1)  # 標本の不偏分散
  return length, smean, ss, ubv

# DataFrame の形成
DF = pd.DataFrame([statistic(tick1), statistic(tick2)],
                  index = ['tick1','tick2'],
                  columns = ['Numbers','Mean','Standard Dev','Unbiased Var']).T

F_value_h = DF['tick2']['Unbiased Var']/DF['tick1']['Unbiased Var'] # 帰無仮説下でのF値
F_value = stats.f.ppf(0.95, DF['tick2']['Numbers']-1,  DF['tick1']['Numbers']-1) # 上側確率のF値

print(f'帰無仮説（２つの時計にばらつき無し）を危険率{prob*100}％で片側検定。')
print(f'帰無仮説下でのF値：{round(F_value_h,2)} \n上側確率{prob*100}％のF値：{round(F_value,2)}')

if F_value_h < F_value:
  print('帰無仮説は維持。')
else:
  print('帰無仮説を棄却。')

DF

帰無仮説（２つの時計にばらつき無し）を危険率5.0％で片側検定。
帰無仮説下でのF値：5.0 
上側確率5.0％のF値：4.1
帰無仮説を棄却。
				tick1	tick2
Numbers			7.0		10.0
Mean			70.0 	50.0
Standard Dev 		2.0 		4.6
Unbiased Var 		4.7 		23.5

p130　パーセンテージを検定する

import math

mtv = 230 # 該番組を見たTV台数
ntv = 618 # TV台数
pratings = 0.3 # 母集団の視聴率（仮定）
sratings = round(mtv/ntv,3) # 標本の視聴率
prob = 0.05 # 危険率

kh = (sratings - pratings) / math.sqrt(sratings*(1-sratings)/ntv) # 帰無仮説下でのz値
kt = stats.norm.ppf(1-prob) # 上側確率のz値


print(f'標本の視聴率：{sratings*100}％。\n帰無仮説（母集団の視聴率{pratings*100}％）を危険率{prob*100}％で検定。')
print(f'帰無仮説下でのz値：{round(kh,3)}\n上側確率{prob*100}％のz値：{round(kt,3)}')

if kh < kt:
  print('帰無仮説は維持。')
else:
  print('帰無仮説を棄却。')

標本の視聴率：37.2％。
帰無仮説（母集団の視聴率30.0％）を危険率5.0％で検定。
帰無仮説下でのz値：3.703
上側確率5.0％のz値：1.645
帰無仮説を棄却。