Up

心理尺度の信頼性係数ω

ベイズ分析

このウェブサイトで信頼性係数オメガ（ω）のベイズ分析を扱っているが、Pystan 3でArvizが使えるようになったので更新版を用意した（2025.03.22）。Arvizを用いて、Rhatの計算、およびトレース図の作成を行った。ファイルは圧縮ファイルupdate2503files.zipにまとめた。

Minicondaのライブラリの更新は、使用している仮想環境においてコマンド

conda update –all

を実行すればよい。

Pystanの更新は、コマンド

pip install –U pystan

実行すればよい。アンインストールしてから再度インストールするという方法もある。

心理尺度の信頼性係数としてオメガω（McDonald, 1999）は、項目jの得点の因子モデル

に基づいて算出される（岡本、2014、第８章）。ここで、は測定対象の値（真値）、は項目jに固有の特性など誤差項としてまとめられるものである。は平均０、分散１に標準化されたものを考えるが、は素点である。

項目が項目あるとし、それらの合計点がテスト得点として算出されるとする。

である。このとき、オメガは次式で与えられる。

の平均が０、分散が１であることから式（３）は

と書ける。式番号（6.20a）はMcDonald (1999)におけるものである。

さらに、各が他の変数と独立であるとすると、式（6.20a）は

と書ける。

各項目の真値の係数（因子負荷量）が同じであるとき、アルファ係数の公式が導かれる。

ベイズ分析例

ベイズ分析のためのStanスクリプトを用いる例をリスト１に、リスト１によるMCMCサンプルの分析を行うスクリプトをリスト２に示した。ファイル類は、まとめて圧縮ファイルupdate2503files.zipとした。

リスト１のスクリプトは、WindowsにインストールされたUbuntu上で実行された。リスト２は、WindowsにインストールされたPythonで実行された。Windows上ではPyStan3は直接実行できないので、WindowsにWSL（Windows Subsystem for Linux）をインストールして、その上にUbuntuをインストールした。WindowsへのUbuntuのインストールはこのウェブサイトで説明している。PyStan3の簡単な説明のウェブサイトも用意している。WSLでは描画ができなかったので、リスト１で生成したMCMCサンプルの描画は、直接Windows上でPython 3.13により行った。リスト１のスクリプトファイルとリスト２のスクリプトファイルを同じフォルダに置いておくと、リスト１のUbuntuにおける実行に続いてリスト２をWindowsで実行することがスムーズに行うことができる。

リスト１のStanスクリプトは、式（１）

のモデルに対するものであるが、式（１）におけるTはリスト１ではFで表されている。Tはtrue valueを意味しているが、Fはfactor scoreであり、すなわち真値である。

パラメータFは、平均０，標準偏差１の標準正規分布に従うという仮定が置かれているが、事前分布での制約ではshrinkageが起こり、このためが過大評価されるという傾向が見られる。これを避けるために、Fの標準化を行っている。

また、MCMCサンプリングにおいてlabel switchingを避けるためにパラメータに非負の制約を置いている。

事前分布は、soft constraintを用いたが、Gelmanら（2021、p. 124）より緩いものを用いた。これは、Gelmanらでは、データ数がsparseの場合にも対応できるように設定されているが、尺度の信頼性を求めるときのデータはある程度の多さがあると考えられるからである。データ数が少なくてパラメータ推定が不安定な場合は、それに応じて事前分布を調整すればよい。PyStan 3では、事前分布として分散の大きいものを設定すると警告が出るので、一様分布に近いものとして分散の大きい分布を用いるときは、まず警告の出ない事前分布のパラメータを用意して、これを変換して分散の大きい分布としている（rescaling）。

リスト１はPyStan 3対応である。PyStan 3をWindowsで利用するときは、WindowsにUbuntuをインストールしてUbuntu上で利用するが、UbuntuのWindowsへのインストール、およびUbuntuでのPythonの準備は簡単なコマンドで出来る。このウェブサイトで説明しているので、必要ならば参照されたい。

リスト１のスクリプトを実行すると、入力データファイル名を聞いてくる。

(py312) *****/omegapystan3files$ python EstOmega_ps3.py

Input Excel data file (*.xlsx) =

入力データファイルはExcelファイルとして図１のように用意する。なお、図１のデータの再分析もこのウェブサイトの後半に示している。

図１

１行目に変数名を設定して、次の２行目からデータを設定する。ただし、１列目はデータの識別用値を設定する。１列目の値は、分析では無視される。

図１のファイルを入力データするときは、入力データファイル名を次のように設定する。

Input Excel data file (*.xlsx) = Data10PS3.xlsx

ファイル名を設定すると、計算が始まる。MCMCサンプリングが終了すると、出力ファイルResults.txtが保存されたことを表示してスクリプトの実行が終了する。

Mu[1]: Q1 = 50.773, Med = 51.782, Q3 = 52.784

Mu[2]: Q1 = 53.389, Med = 54.234, Q3 = 55.089

Mu[3]: Q1 = 49.044, Med = 50.187, Q3 = 51.332

Mu[4]: Q1 = 51.570, Med = 52.644, Q3 = 53.673

Mu[5]: Q1 = 51.922, Med = 52.893, Q3 = 53.820

Mu[6]: Q1 = 44.403, Med = 46.061, Q3 = 47.776

Mu[7]: Q1 = 51.574, Med = 52.507, Q3 = 53.406

Mu[8]: Q1 = 51.353, Med = 52.288, Q3 = 53.155

Mu[9]: Q1 = 42.242, Med = 43.764, Q3 = 45.282

Mu[10]: Q1 = 53.064, Med = 54.047, Q3 = 54.995

Results.txt was saved.

(py312) *****/omegapystan3files$

出力ファイルResults.txtは、エディタで開くことができる。

リスト１のスクリプトの実行中に保存された他のファイルは、次のリスト２のスクリプトの実行で読み込まれて分析される。リスト１とリスト２は同じフォルダ内に置かれているものとする。

リスト２のスクリプトを実行すると、まずトレース図が表示される（図２）。

図２

図２のWindowを閉じると、次にの事後分布のグラフが表示される（図２）。

図３

図３のWindowを閉じると、の事後分布が表示される（図４）。

図４

図４のWindowを閉じると、の事後分布のグラフが表示される（図５）。

図５

図５のWindowを閉じると、次に、変数ごとに各変数における真値の分散の各変数の分散に対する比

の事後分布のグラフが表示される（図６）。この値は、各変数に占める真値の割合と考えられる。

図６

図６のWindowを閉じると、ωの事後分布のグラフが表示される（図７）。

図７

グラフのタイトル部に、事後分布の4分位数が表示されている。点推定値としては、MAP推定値（モード）がよいという報告（岡本、2022）があるが、計算時間が掛かるので、ここではGelmanら（2021）に従って、中央値を点推定値とした。

図７のフォームを閉じると、図８の散布図が表示される。

図８

横軸に各の推定値（事後分布の中央値）、縦軸にテストスコアの和をとり、プロットしたものである。点がほぼ直線沿いに散らばっていることが判る。

図８のWindowを閉じると、スクリプトの実行終了となる。端末画面には、次のメッセージが表示されている。

Results_1.txt was saved.

出力ファイルResults_1.txtを開くと、以下のようになっている。

data...

[[ 1 79 46 68 65 89 38 58 70 49 65]

[ 2 48 71 61 50 43 27 48 54 34 61]

[ 3 38 42 38 45 55 43 47 39 65 54]

　　　　・

　　　　・

　　　　・

[48 36 46 49 54 31 75 43 40 30 49]

[49 44 24 56 51 58 34 42 52 67 48]

[50 71 75 64 84 58 25 67 65 47 56]]

Mu...

Mu[1]: Q1 = 50.785, Med = 51.810, Q3 = 52.811

Mu[2]: Q1 = 53.429, Med = 54.265, Q3 = 55.114

Mu[3]: Q1 = 49.031, Med = 50.209, Q3 = 51.384

Mu[4]: Q1 = 51.652, Med = 52.643, Q3 = 53.681

Mu[5]: Q1 = 51.983, Med = 52.912, Q3 = 53.871

Mu[6]: Q1 = 44.529, Med = 46.142, Q3 = 47.828

Mu[7]: Q1 = 51.595, Med = 52.498, Q3 = 53.336

Mu[8]: Q1 = 51.352, Med = 52.292, Q3 = 53.156

Mu[9]: Q1 = 42.201, Med = 43.731, Q3 = 45.296

Mu[10]: Q1 = 53.069, Med = 54.031, Q3 = 54.940

Sigma...

Sigma[1]: Q1 = 9.680 Med = 10.458, Q3 = 11.357

Sigma[2]: Q1 = 8.264 Med = 8.951, Q3 = 9.690

Sigma[3]: Q1 = 11.033 Med = 11.869, Q3 = 12.805

Sigma[4]: Q1 = 9.843 Med = 10.703, Q3 = 11.631

Sigma[5]: Q1 = 9.271 Med = 10.023, Q3 = 10.876

Sigma[6]: Q1 = 16.312 Med = 17.457, Q3 = 18.752

Sigma[7]: Q1 = 8.560 Med = 9.314, Q3 = 10.053

Sigma[8]: Q1 = 8.603 Med = 9.378, Q3 = 10.161

Sigma[9]: Q1 = 14.591 Med = 15.630, Q3 = 16.777

Sigma[10]: Q1 = 9.414 Med = 10.234, Q3 = 11.086

Lambda...

Lambda[1]: Q1 = 11.953 Med = 13.005, Q3 = 14.003 (Item_1)

Lambda[2]: Q1 = 11.723 Med = 12.621, Q3 = 13.448 (Item_2)

Lambda[3]: Q1 = 11.543 Med = 12.708, Q3 = 13.885 (Item_3)

Lambda[4]: Q1 = 15.081 Med = 16.122, Q3 = 17.148 (Item_4)

Lambda[5]: Q1 = 12.303 Med = 13.285, Q3 = 14.276 (Item_5)

Lambda[6]: Q1 = 0.578 Med = 1.263, Q3 = 2.212 (Item_6)

Lambda[7]: Q1 = 11.286 Med = 12.228, Q3 = 13.128 (Item_7)

Lambda[8]: Q1 = 14.821 Med = 15.733, Q3 = 16.678 (Item_8)

Lambda[9]: Q1 = 1.179 Med = 2.279, Q3 = 3.585 (Item_9)

Lambda[10]: Q1 = 14.781 Med = 15.841, Q3 = 16.864 (Item_10)

Lambda_j**/Var(X_j)...

Lambda[1]**2 / Var(X[1]): Q1 = 0.545, Med = 0.608, Q3 = 0.660

Lambda[2]**2 / Var(X[2]): Q1 = 0.611, Med = 0.667, Q3 = 0.711

Lambda[3]**2 / Var(X[3]): Q1 = 0.470, Med = 0.535, Q3 = 0.593

Lambda[4]**2 / Var(X[4]): Q1 = 0.645, Med = 0.695, Q3 = 0.739

Lambda[5]**2 / Var(X[5]): Q1 = 0.579, Med = 0.639, Q3 = 0.687

Lambda[6]**2 / Var(X[6]): Q1 = 0.001, Med = 0.005, Q3 = 0.016

Lambda[7]**2 / Var(X[7]): Q1 = 0.576, Med = 0.635, Q3 = 0.686

Lambda[8]**2 / Var(X[8]): Q1 = 0.692, Med = 0.738, Q3 = 0.780

Lambda[9]**2 / Var(X[9]): Q1 = 0.006, Med = 0.021, Q3 = 0.050

Lambda[10]**2 / Var(X[10]): Q1 = 0.655, Med = 0.707, Q3 = 0.749

omega...

Q1 = 0.897, Med = 0.904, Q3 = 0.910

まず、入力データが出力されている。

続いて、MCMCサンプリングの統計量が出力されている。

事後分布などのグラフで表示されていた情報が、統計量の数値データとして読むことができる。

図５（）、図６（）における他とは著しく異なる項目は、上の統計量から６番目と９番目の項目であることが解る。６番目と９番目を除いて再分析を行う。

再分析

図１のデータから、６番目と９番目の項目のデータを抜いたものを、図A1に示すように用意した。

図A1　　図１のデータから６番目と９番目のデータを抜いたもの

図A1のデータに対して、まずリスト１のスクリプトでMCMCサンプリングを行った。

(py312) *****/update2503files$ python EstOmega_ps3.py

Input Excel data file (*.xlsx) = Data10PS3_69.xlsx

リスト１のスクリプトによるMCMCサンプリングの後、MCMCサンプルの分析をリスト２のスクリプトで行った。

トレース図hは、図A2のようである。

図A2

の事後分布は、図A３のようである。

図A３

の事後分布は図A４のようである。

図A４

の事後分布は、図A５のようである。

図A５

の事後分布は図A６のようである。

図A６

ωの事後分布は図A７のようである。図７のものより右に寄っていることが分かる。

図A７

点の散布図は図A８のようである。図８のものより、点が狭い帯内にまとまっていることが分かる。

図A８

出力ファイルResults_1.txtを開くと、以下のようになっている。

data...

[[ 1 79 46 68 65 89 58 70 65]

[ 2 48 71 61 50 43 48 54 61]

[ 3 38 42 38 45 55 47 39 54]

　　　　・

　　　　・

　　　　・

[48 36 46 49 54 31 43 40 49]

[49 44 24 56 51 58 42 52 48]

[50 71 75 64 84 58 67 65 56]]

Mu...

Mu[1]: Q1 = 50.837, Med = 51.781, Q3 = 52.833

Mu[2]: Q1 = 53.409, Med = 54.240, Q3 = 55.084

Mu[3]: Q1 = 49.092, Med = 50.200, Q3 = 51.289

Mu[4]: Q1 = 51.647, Med = 52.661, Q3 = 53.699

Mu[5]: Q1 = 51.935, Med = 52.908, Q3 = 53.903

Mu[6]: Q1 = 51.618, Med = 52.480, Q3 = 53.363

Mu[7]: Q1 = 51.373, Med = 52.243, Q3 = 53.147

Mu[8]: Q1 = 53.064, Med = 54.042, Q3 = 55.040

Sigma...

Sigma[1]: Q1 = 9.669 Med = 10.431, Q3 = 11.271

Sigma[2]: Q1 = 8.236 Med = 8.952, Q3 = 9.683

Sigma[3]: Q1 = 11.016 Med = 11.874, Q3 = 12.810

Sigma[4]: Q1 = 9.845 Med = 10.723, Q3 = 11.647

Sigma[5]: Q1 = 9.253 Med = 10.007, Q3 = 10.889

Sigma[6]: Q1 = 8.568 Med = 9.263, Q3 = 10.051

Sigma[7]: Q1 = 8.611 Med = 9.388, Q3 = 10.189

Sigma[8]: Q1 = 9.471 Med = 10.298, Q3 = 11.157

Lambda...

Lambda[1]: Q1 = 11.930 Med = 13.029, Q3 = 14.067 (Item_1)

Lambda[2]: Q1 = 11.698 Med = 12.637, Q3 = 13.470 (Item_2)

Lambda[3]: Q1 = 11.541 Med = 12.670, Q3 = 13.855 (Item_3)

Lambda[4]: Q1 = 15.021 Med = 16.109, Q3 = 17.221 (Item_4)

Lambda[5]: Q1 = 12.230 Med = 13.254, Q3 = 14.207 (Item_5)

Lambda[6]: Q1 = 11.353 Med = 12.278, Q3 = 13.236 (Item_7)

Lambda[7]: Q1 = 14.741 Med = 15.713, Q3 = 16.623 (Item_8)

Lambda[8]: Q1 = 14.805 Med = 15.822, Q3 = 16.840 (Item_10)

Lambda_j**/Var(X_j)...

Lambda[1]**2 / Var(X[1]): Q1 = 0.547, Med = 0.609, Q3 = 0.662

Lambda[2]**2 / Var(X[2]): Q1 = 0.611, Med = 0.666, Q3 = 0.714

Lambda[3]**2 / Var(X[3]): Q1 = 0.470, Med = 0.532, Q3 = 0.592

Lambda[4]**2 / Var(X[4]): Q1 = 0.643, Med = 0.693, Q3 = 0.738

Lambda[5]**2 / Var(X[5]): Q1 = 0.577, Med = 0.634, Q3 = 0.687

Lambda[6]**2 / Var(X[6]): Q1 = 0.582, Med = 0.639, Q3 = 0.687

Lambda[7]**2 / Var(X[7]): Q1 = 0.691, Med = 0.738, Q3 = 0.776

Lambda[8]**2 / Var(X[8]): Q1 = 0.655, Med = 0.705, Q3 = 0.748

omega...

Q1 = 0.931, Med = 0.936, Q3 = 0.940

他と特に互いに異なる値の項目は認められない。ωの値も高くなっており

１０項目のときの

Q1 = 0.897, Med = 0.904, Q3 = 0.910

から、第６項目と第９項目を抜いた８項目でのωは

Q1 = 0.931, Med = 0.936, Q3 = 0.940

であり、四分位数の比較でも８項目に減らしたときのωの方が高い値であると言えよう。

リスト１　StanによるMCMCサンプリング

import numpy as np

import pandas as pd

import stan

import pickle

import arviz as az

StanCode = """

//

// Yasuharu Okamoto, 2023.08, 2025.03

//

data {

int N;

int M;

array[N,M] real X;

vector[M] i_mu;

vector[M] i_std;

}

parameters {

vector[N] vF;

vector[M] v_mu;

vector<lower=0.0>[M] v_Lmbd;

vector<lower=0.0>[M] v_sgm;

}

transformed parameters {

vector[N] F;

real sumF;

real ssF;

real meanF;

real sdF;

vector[M] mu;

vector[M] Lmbd;

vector[M] sgm;

array[N] vector[M] V;

//

// Standardized F values

//

sumF = 0.0;

ssF = 0.0;

for (i in 1:N) {

sumF += vF[i];

ssF += pow(vF[i], 2);

}

meanF = sumF / N;

sdF = pow((ssF - N * pow(meanF, 2)) / (N - 1), 0.5);

F = (vF - meanF) / sdF; // 標準化

//

// Rescaled parameters

//

mu = i_mu + 10.0 .* i_std .* v_mu;

Lmbd = 10.0 .* i_std .* v_Lmbd;

sgm = 0.001 + 10.0 .* i_std .* v_sgm;

//

// Regression model

//

for (i in 1:N)

V[i] = mu + Lmbd * F[i];

}

model {

vF ~ normal(0.0, 1.0);

v_mu ~ normal(0.0, 10.0);

v_Lmbd ~ exponential(0.1);

v_sgm ~ exponential(0.1);

for (i in 1:N)

X[i] ~ normal(V[i], sgm);

}

generated quantities {

real omg; // omega

vector[M] L2V;

real v1;

real v2;

v1 = 0.0;

v2 = 0.0;

L2V = pow(Lmbd, 2.0) ./ (pow(Lmbd, 2.0) + pow(sgm, 2.0));

for (j in 1:M) {

v1 += Lmbd[j];

v2 += pow(sgm[j], 2);

}

omg = pow(v1, 2) / (pow(v1, 2) + v2);

}

"""

pd.set_option('display.max_columns', 10000)

pd.set_option('display.max_rows', 10000)

fout = open('Results.txt', 'w')

#

# Excelデータファイルの読み込み

#

s = input("Input Excel data file (*.xlsx) = ")

xlsx_fl = pd.ExcelFile(s)

data_xlsx = pd.read_excel(xlsx_fl)

with open('data.pkl', 'wb') as f:

pickle.dump(data_xlsx, f)

var_names = [v for v in data_xlsx.keys()][1:]

print('var_names =', var_names)

data = data_xlsx.values

print('data...\n', data)

fout.write(f'Input data file = {s}\n')

fout.write(f'\ndata...\n{data}\n')

X = data[:,1:]

mus = np.mean(X, axis=0)

stds = np.std(X, axis=0)

Data = {'N': len(X), 'M': len(X[0]), 'X': X, 'i_mu': mus, 'i_std': stds}

sm = stan.build(StanCode, data = Data)

fit = sm.sample()

df = fit.to_frame() # DataFrame型への変換

with open('fit.pkl', 'wb') as f:

pickle.dump(df, f)

print()

fout.write('\n\nMu...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'mu.{j+1}'], [25, 50, 75])

print(f'Mu[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Mu[{j+1}]: Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

i_data = az.from_pystan(posterior=fit, posterior_model=sm)

with open('i_data.pkl', 'wb') as f:

pickle.dump(i_data, f) # InferenceData of MCMC samples is saved

print(az.summary(i_data)) # Ptint summary of MCMC sample

fout.write(f'\n{az.summary(i_data)}\n')

fout.close()

print('\nResults.txt was saved.\n')

リスト２　MCMCサンプルの分析

import numpy as np

import pandas as pd

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

fout = open('Results_1.txt', 'w')

with open('data.pkl', 'rb') as f: # 保存されたデータの読み込み

data_xlsx = pickle.load(f)

var_names = [v for v in data_xlsx.keys()][1:]

print('var_names =', var_names)

data = data_xlsx.values

print('data...\n', data)

fout.write(f'\ndata...\n{data}\n')

X = data[:,1:]

mus = np.mean(X, axis=0)

stds = np.std(X, axis=0)

# 保存されたMCMCサンプルの読み込み

with open('fit.pkl', 'rb') as f:

df = pickle.load(f)

with open('i_data.pkl', 'rb') as f:

i_data = pickle.load(f)

"""

MCMCサンプリングのトレース図

"""

az.plot_trace(i_data, var_names=['Lmbd', 'sgm', 'omg'], figsize=(10,7))

plt.tight_layout()

plt.show()

#

# パラメータmuの4分位数の計算

#

print()

fout.write('\n\nMu...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'mu.{j+1}'], [25, 50, 75])

print(f'Mu[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Mu[{j+1}]: Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

#

# muの事後分布

#

for j in range(len(X[0])):

sb.kdeplot(df[f'mu.{j+1}'], label=f'{var_names[j]}')

plt.title(r'$\mu$', fontsize = 16)

plt.xlabel('')

plt.legend(bbox_to_anchor=(1.0, 1.0))

plt.tight_layout()

plt.show()

#

# Sigmaの4分位数の計算

#

print()

fout.write('\nSigma...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'sgm.{j+1}'], [25, 50, 75])

print(fr'Sigma[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Sigma[{j+1}]: Q1 = {q1:.3f}' +

f' Med = {med:.3f}, Q3 = {q3:.3f}\n')

#

# Sigmaの事後分布

#

for j in range(len(X[0])):

sb.kdeplot(df[f'sgm.{j+1}'], label=f'{var_names[j]}')

plt.title(r'$\sigma$', fontsize = 16)

plt.xlabel('')

plt.legend(bbox_to_anchor=(1.0, 1.0))

plt.tight_layout()

plt.show()

#

# Lambdaの四分位数の計算

#

print()

fout.write('\nLambda...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'Lmbd.{j+1}'], [25, 50, 75])

print(fr'Lambda[{j+1}]:', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}', end = '')

print(f' ({var_names[j]})')

fout.write(f'Lambda[{j+1}]: Q1 = {q1:.3f}' +

f' Med = {med:.3f}, Q3 = {q3:.3f}' +

f' ({var_names[j]})\n')

#

# Lambdaの事後分布

#

for j in range(len(X[0])):

sb.kdeplot(df[f'Lmbd.{j+1}'], label=f'{var_names[j]}')

plt.title(r'$\lambda$', fontsize = 16)

plt.xlabel('')

plt.legend(bbox_to_anchor=(1.0, 1.0))

plt.tight_layout()

plt.savefig('FigLmbd.png')

plt.show()

#

# 真値の割合の計算

#

print()

fout.write('\nLambda_j**/Var(X_j)...\n')

for j in range(len(X[0])):

q1, med, q3 = np.percentile(df[f'L2V.{j+1}'], [25, 50, 75])

print(f'Lambda[{j+1}]**2 / Var(X[{j+1}]):', end = '')

print(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write(f'Lambda[{j+1}]**2 / Var(X[{j+1}]):' +

f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

for j in range(len(X[0])):

sb.kdeplot(df[f'L2V.{j+1}'], label=f'{var_names[j]}')

plt.title(r'$\lambda^{2}_{j}/Var(X_{j})$', fontsize = 16)

plt.xlabel('')

plt.legend(bbox_to_anchor=(1.0, 1.0))

plt.tight_layout()

plt.show()

#

# Omegaの四分位数の計算

#

print()

q1, med, q3 = np.percentile(df[f'omg'], [25, 50, 75])

print(f'omega: Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}')

fout.write('\nomega...\n')

fout.write(f' Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}\n')

#

# Omegaの事後分布

#

sb.kdeplot(df['omg'])

plt.title(f'Q1 = {q1:.3f}, Med = {med:.3f}, Q3 = {q3:.3f}', fontsize = 14)

plt.xlabel(r'$\omega$', fontsize = 16)

plt.tight_layout()

plt.show()

#

# （因子得点F、尺度得点和）の散布図

#

sumXs = np.sum(X, axis = 1)

Fs = []

for i in range(len(X)):

Fs.append(np.median(df[f'F.{i+1}']))

plt.plot(Fs, sumXs, 'o')

plt.xlabel(r'$\widehat{F}_{median}$', fontsize = 14)

plt.ylabel(r'$\sum{X_{j}}$', fontsize = 14)

plt.tight_layout()

plt.show()

fout.close()

print('\nResults_1.txt was saved.\n')

引用文献

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.

McDonald, R. P. (1999). Test Theory: A unified treatment. Psychology Press.

Home