Up

非心エフ分布のシミュレーション

非心エフ分布のシミュレーションを試みた。非心エフ分布に従う確率変数は、非心カイ2乗分布に従う確率変数をその自由度で割ったものを、別のカイ2乗分布に従う確率変数をその自由度で割ったもので割ることにより得ることができる。非心エフ分布のシミュレーションプログラムをリスト１に示す。

リスト１　非心エフ分布のシミュレーション

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

import seaborn as sb

n_smpls = 1_000_000

df1 = int(input('df1 = '))

df2 = int(input('df2 = '))

lmbd = float(input('noncentrality = '))

X1 = ss.ncx2.rvs(df1, lmbd, size = n_smpls)

X2 = ss.chi2.rvs(df2, size = n_smpls)

F = (X1/df1) / (X2/df2)

F = np.sort(F)

F = F[int(0.025*len(F)):int(0.975*len(F))] # 左右の裾2.5%ずつを除く

plt.hist(F, bins = 100, density = True, alpha = 0.5, label = 'Simulation')

x_coord = np.linspace(F[0], F[-1], 1000)

#

#　　　サンプルの５％除いた処理に合わせる

#

y_coord = ss.ncf.pdf(x_coord, df1, df2, lmbd) / 0.95

plt.plot(x_coord, y_coord, label = 'Model')

plt.title(f'df1 = {df1}, df2 = {df2}, noncentrality = {lmbd}')

plt.legend()

plt.show()

リスト１のスクリプトを実行すると、自由度と非心度の入力が求められる。

df1 = 3

df2 = 7

noncentrality = 2

非心度を入力すると、計算が始まり、計算後、シミュレーションによるサンプルのヒストグラムと非心エフ分布のグラフが表示される（図１）。

図１

サンプルのヒストグラムは、分布の左右の裾2.5％ずつが除かれていて、非心エフ分布のグラフはそれに合わせて調整されている。

非心エフ分布の説明を＜岡本安晴「データを使いこなすための統計入門」22世紀アート＞の付録　確率：E.9　非心エフ分布で行っている。

Up