平均値の条件差検定のシミュレーション
データを母集団からのランダムなサンプリングにより得られたと考えるとき、得られるデータとその統計量はサンプリングごとに変動する。変動の様子を見るシミュレーションプログラムをリスト１に示すように作成してみた。プログラムファイルは、圧縮ファイルsimu2popfile.zipにまとめたので、ダウンロード展開すれば自由に利用できる。ただし、これはWindows用である。CentOS Stream 9などのOS用は、このウェブサイトの後半で説明している。
２つのデータセットx1およびx2がそれぞれ平均
および
、標準偏差
および
の正規分布からの
個および
個のサンプルであるとする。
このときのデータのサンプリングを次のスクリプトで行った。
x1 = ss.norm.rvs(mu1, sgm1,
n1)
x2 = ss.norm.rvs(mu2, sgm2,
n2)
このデータに対して、平均値と標準偏差を次のスクリプトで求めた。
m1_hat = np.mean(x1)
s1_hat = np.std(x1)
m2_hat = np.mean(x2)
s2_hat = np.std(x2)
また、ｔ検定の値を、次のスクリプトで求めた。
r = ss.ttest_ind(x1,
x2)
# Student
rw = ss.ttest_ind(x1, x2,
equal_var = False) # Welch
サンプルされたデータと統計量、検定結果などの情報は、ライブラリPyGameを利用してシミュレーションのウィンドウ画面上に表示した。PyGameのインストールは、コマンド
pip install
pygame
の実行で簡単に行うことができた。
リスト１では、表示ウィンドウを前面に出すためにモジュールwin32guiを用いているが、これはpywin32のインストールによりimportできるようになる。pywin32のインストールは、まず次のコマンドを実行する。
pip install
pywin32
続いて、次のコマンドをPythonのインストールされているフォルダに移って実行する。
C:\Pythonのインストールされているフォルダ\python.exe
Scripts\pywin32_postinstall.py -install
これで、再起動後、次のようにwin32guiがimportできる。
import
win32gui
ただし、Anacondaの場合は、上のpipによるインストールで完了である。\python.exe Scripts\pywin32_postinstall.py -installによるインストール作業は不要である。
リスト１のプログラムを実行すると、まず、２つの母集団の平均と標準偏差、およびそれぞれのサンプルサイズの入力が求められる（図１）。
図１
図１の状態で、求められたすべての値を設定すると、シミュレーションが始まるが、その前に図２のウィンドウが表示される。
図２
これは、表示ウィンドウにキーボード入力が向けられるようにするためである。
図２のウィンドウの画面は、２回クリックすると消されて、シミュレーションが始まり、計算結果が図３のように表示される。
図３
最上段に、標準ｔ検定の結果とWelchのｔ検定の結果が表示されている。分散が等しいかどうか不明の場合は、Welchのｔ検定を用いる方がよいと思われる。
次の段に、サンプルx1の母集団パラメータと統計量が表示され、x1の分布が緑の小円で描かれている。
その下の段には、サンプルx2の母集団パラメータと統計量、およびx2の分布が青の小円で描かれている。
最下段に、キーの選択が示されている。Ｎキーを押すと、次のサンプリングが行われて、結果が表示される。
下図に、Nキーを押して、さらに４回のシミュレーションを行った例を示す。
Eキー押すと、プログラムの実行終了である。
なお、Pythonによるデータ分析の入門書＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を上梓した。
リスト１ ２つの正規分布からのサンプリング（simu_2pop.py）
import pygame as pg
import scipy.stats as ss
import numpy as np
from win32gui import SetWindowPos,
SetFocus
import win32con
mu1 = float(input('mu_1 = '))
mu2 = float(input('mu_2 = '))
sgm1 = float(input('sgm_1 = '))
sgm2 = float(input('sgm_2 = '))
n1 = int(input('n_1 = '))
n2 = int(input('n_2 = '))
wt = int(3000 / (n1 + n2)) + 1
pg.init()
pg.display.set_caption('Two
Independent Sampling')
screen =
pg.display.set_mode(size=(800, 600))
SetWindowPos(pg.display.get_wm_info()['window'],
-1, 0, 0, 0, 0,
win32con.SWP_ASYNCWINDOWPOS | win32con.SWP_NOMOVE |
win32con.SWP_NOSIZE | win32con.SWP_SHOWWINDOW)
screen.fill((255,255,255))
pg.display.flip()
my_font = pg.font.SysFont(None, 32)
my_font_L = pg.font.SysFont(None, 40)
img_str_L = my_font_L.render('Next
-> key N, End ->
key E', True,
(0,0,0), (255,255,255))
rect_str_L = img_str_L.get_rect()
my_font_L1 = pg.font.SysFont(None, 60)
img_str_mouse =
my_font_L1.render('Click here, twice', True,
(0,0,0), (255,255,255))
rect_str_mouse =
img_str_mouse.get_rect()
screen.blit(img_str_mouse, (400 -
rect_str_mouse.width//2,
300 - rect_str_mouse.height//2))
pg.display.flip()
ck1 = True
while ck1:
for e in pg.event.get():
if
e.type == pg.MOUSEBUTTONDOWN:
ck1 = False
ck = True
while ck:
x1 = ss.norm.rvs(mu1, sgm1,
n1)
x2 = ss.norm.rvs(mu2, sgm2,
n2)
min_x = np.min(x1) if
np.min(x1) < np.min(x2) else np.min(x2)
max_x = np.max(x1) if
np.max(x1) > np.max(x2) else np.max(x2)
def scale_x(v):
return int(30 + 740 * (v - min_x) / (max_x - min_x))
screen.fill((255,255,255))
for i in range(n1):
pg.draw.circle(screen, (0,200,0), (scale_x(x1[i]), 250), 10, 3)
pg.display.flip()
pg.time.wait(wt)
m1_hat = np.mean(x1)
s1_hat = np.std(x1)
str1 = f'mu1 = {mu1}, sgm1 =
{sgm1},' + \
f' m1_hat = {m1_hat:.2f},
s1_hat = {s1_hat:.2f}'
img_str1 =
my_font.render(str1, True, (0,0,0), (255,255,255))
screen.blit(img_str1, (50,
150))
for i in range(n2):
pg.draw.circle(screen, (0,0,200), (scale_x(x2[i]), 400), 10, 3)
pg.display.flip()
pg.time.wait(wt)
m2_hat = np.mean(x2)
s2_hat = np.std(x2)
str2 = f'mu2 = {mu2}, sgm2 =
{sgm2},' + \
f' m2_hat = {m2_hat:.2f},
s2_hat = {s2_hat:.2f}'
img_str2 =
my_font.render(str2, True, (0,0,0), (255,255,255))
screen.blit(img_str2, (50,
330))
r = ss.ttest_ind(x1,
x2)
# Student
rw = ss.ttest_ind(x1, x2,
equal_var = False) # Welch
strg = f't =
{r.statistic:.3f}, p = {r.pvalue:.5f}, df = {r.df:.0f}'
strgw = f"t' = {rw.statistic:.3f}, p' =
{rw.pvalue:.5f}, df' = {rw.df:.3f}"
img_strg =
my_font.render(strg, True, (0,0,0), (255,255,255))
screen.blit(img_strg, (400 -
(img_strg.get_rect().width//2), 10))
img_strgw =
my_font.render(strgw, True, (0,0,0), (255,255,255))
screen.blit(img_strgw, (400
- (img_strgw.get_rect().width//2), 40))
screen.blit(img_str_L, (400
- (rect_str_L.width//2),
600 - (rect_str_L.height*2)) )
pg.display.flip()
pg.event.set_grab(True)
ck1 = True
while ck1:
for
e in pg.event.get():
if e.type == pg.KEYDOWN:
if e.key == pg.K_n:
ck1 = False
elif e.key == pg.K_e:
ck = False
ck1
= False
pg.quit()
Linux（Ubuntu、CentOS
Stream 9など）の場合
２つのデータセットx1およびx2がそれぞれ平均および、標準偏差およびの正規分布からの個および個のサンプルであるとする。
このときのデータのサンプリングを次のスクリプトで行った。
x1 = ss.norm.rvs(mu1, sgm1,
n1)
x2 = ss.norm.rvs(mu2, sgm2,
n2)
このデータに対して、平均値と標準偏差を次のスクリプトで求めた。
m1_hat = np.mean(x1)
s1_hat = np.std(x1)
m2_hat = np.mean(x2)
s2_hat = np.std(x2)
また、ｔ検定の値を、次のスクリプトで求めた。
r = ss.ttest_ind(x1,
x2)
# Student
rw = ss.ttest_ind(x1, x2, equal_var
= False) # Welch
サンプルされたデータと統計量、検定結果などの情報は、ライブラリPyGameを利用してシミュレーションのウィンドウ画面上に表示した。PyGameのインストールは、コマンド
pip install
pygame
の実行で簡単に行うことができた。
シミュレーションはリスト２のスクリプトで行ったが、スクリプトファイルは圧縮ファイルsimu2pop_cos.zipにまとめた。
リスト２のプログラムを実行すると、まず、２つの母集団の平均と標準偏差、およびそれぞれのサンプルサイズの入力が求められる。
mu_1 = 70
mu_2 = 80
sgm_1 = 10
sgm_2 = 12
n_1 = 50
n_2 = 55
これらの値の設定後、シミュレーションが始まる。シミュレーションの結果は、下図のように表示される。
最上段に、標準ｔ検定の結果とWelchのｔ検定の結果が表示されている。分散が等しいかどうか不明の場合は、Welchのｔ検定を用いる方がよいと思われる。
次の段に、サンプルx1の母集団パラメータと統計量が表示され、x1の分布が緑の小円で描かれている。
その下の段には、サンプルx2の母集団パラメータと統計量、およびx2の分布が青の小円で描かれている。
最下段に、キーの選択が示されている。Ｎキーを押すと、次のサンプリングが行われて、結果が表示される。
Eキー押すと、プログラムの実行終了である。
リスト２のスクリプトは、CentOS Stream 9のPython 3.12.3/Anacondaで実行したが、Ubuntu 24/WSL/Windows 10のPython 3.12.3/Anacondaでも実行できた。
なお、Pythonによるデータ分析の入門書＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を上梓した。
リスト２ Linux（Ubuntu、CentOS Stream
9など）用サンプルスクリプト
import pygame as pg
import scipy.stats as ss
import numpy as np
mu1 = float(input('mu_1 = '))
mu2 = float(input('mu_2 = '))
sgm1 = float(input('sgm_1 = '))
sgm2 = float(input('sgm_2 = '))
n1 = int(input('n_1 = '))
n2 = int(input('n_2 = '))
wt = int(3000 / (n1 + n2)) + 1
pg.init()
pg.display.set_caption('Two
Independent Sampling')
screen =
pg.display.set_mode(size=(800, 600))
screen.fill((255,255,255))
pg.display.flip()
my_font = pg.font.SysFont(None, 32)
my_font_L = pg.font.SysFont(None, 40)
img_str_L = my_font_L.render('Next
-> key N, End ->
key E', True,
(0,0,0), (255,255,255))
rect_str_L = img_str_L.get_rect()
ck = True
while ck:
x1 = ss.norm.rvs(mu1, sgm1,
n1)
x2 = ss.norm.rvs(mu2, sgm2, n2)
min_x = np.min(x1) if
np.min(x1) < np.min(x2) else np.min(x2)
max_x = np.max(x1) if
np.max(x1) > np.max(x2) else np.max(x2)
def scale_x(v):
return int(30 + 740 * (v - min_x) / (max_x - min_x))
screen.fill((255,255,255))
for i in range(n1):
pg.draw.circle(screen, (0,200,0), (scale_x(x1[i]), 250), 10, 3)
pg.display.flip()
pg.time.wait(wt)
m1_hat = np.mean(x1)
s1_hat = np.std(x1)
str1 = f'mu1 = {mu1}, sgm1 =
{sgm1},' + \
f' m1_hat = {m1_hat:.2f},
s1_hat = {s1_hat:.2f}'
img_str1 =
my_font.render(str1, True, (0,0,0), (255,255,255))
screen.blit(img_str1, (50,
150))
for i in range(n2):
pg.draw.circle(screen, (0,0,200), (scale_x(x2[i]), 400), 10, 3)
pg.display.flip()
pg.time.wait(wt)
m2_hat = np.mean(x2)
s2_hat = np.std(x2)
str2 = f'mu2 = {mu2}, sgm2 =
{sgm2},' + \
f' m2_hat = {m2_hat:.2f},
s2_hat = {s2_hat:.2f}'
img_str2 = my_font.render(str2,
True, (0,0,0), (255,255,255))
screen.blit(img_str2, (50,
330))
r = ss.ttest_ind(x1,
x2)
# Student
rw = ss.ttest_ind(x1, x2,
equal_var = False) # Welch
strg = f't =
{r.statistic:.3f}, p = {r.pvalue:.5f}, df = {r.df:.0f}'
strgw = f"t' = {rw.statistic:.3f}, p' =
{rw.pvalue:.5f}, df' = {rw.df:.3f}"
img_strg =
my_font.render(strg, True, (0,0,0), (255,255,255))
screen.blit(img_strg, (400 -
(img_strg.get_rect().width//2), 10))
img_strgw =
my_font.render(strgw, True, (0,0,0), (255,255,255))
screen.blit(img_strgw, (400
- (img_strgw.get_rect().width//2), 40))
screen.blit(img_str_L, (400
- (rect_str_L.width//2),
600 - (rect_str_L.height*2)) )
pg.display.flip()
pg.event.set_grab(True)
ck1 = True
while ck1:
for
e in pg.event.get():
if e.type == pg.KEYDOWN:
if e.key == pg.K_n:
ck1 = False
elif e.key == pg.K_e:
ck = False
ck1 = False
pg.quit()