非心ティ分布のシミュレーション
非心ティ分布の説明を＜岡本安晴「データを使いこなすための統計入門」22世紀アート＞の付録 確率：E.8 非心ティ分布で行っている。
非心ティ分布に従う確率変数は、標準正規分布に従う確率変数に非心度を加えたものを、カイ２乗分布をその自由度で割ったものの平方根で割ることによって得られる。非心ティ分布のシミュレーションプログラムをリスト１のように作成してみた。
リスト１ 非心ティ分布のシミュレーション
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
n_smpls = 100_000
df = int(input('df = '))
delta = float(input('Noncentrality = '))
Z = ss.norm.rvs(0,
1, size = n_smpls) + delta
Y = ss.chi2.rvs(df,
size = n_smpls)
t = Z / ((Y/df)**0.5)
#
# tのサンプルから両裾のそれぞれ2.5%を除く
#
t = np.sort(t)
LB = int(len(t) * 0.025)
UB = int(len(t) * 0.975)
t = t[LB:UB]
x_coord = np.linspace(t[0],
t[-1], 1000)
#
# ｔのサンプルから５％除いた影響に合わせるため、0.95で割る
#
y_coord = ss.nct.pdf(x_coord,
df, delta) / 0.95
plt.hist(t, bins = 100, density = True, alpha
= 0.5, label = 'Simulation')
plt.plot(x_coord, y_coord, label = 'Model')
plt.title(f'df = {df}, NonCentrality = {delta}')
plt.legend()
plt.show()
リスト１のスクリプトを実行すると、自由度と非心度の入力が求められる。
df = 1
Noncentrality = 5
非心度が入力されると、計算が始まり、計算結果のグラフが表示される（図１）。
図１
シミュレーションによるサンプルのヒストグラムと確率分布モデルの曲線が表示されている。