心理尺度の信頼性係数オメガωのベイズ分析
本ウェブサイトのスクリプトは、CmdStanPyを利用している。CmdStanPyの準備の簡単な説明をこのウェブサイトに用意した。
心理尺度の信頼性係数として、オメガω（McDonald, 1999）は、項目jの得点
の因子モデル
に基づいて算出される（岡本、2014、第８章）。ここで、
は測定対象の因子の値、
は項目jに固有の特性など残差項（誤差項）としてまとめられるものである。は平均０、分散１に標準化されたものを考えるが、は素点である。
項目が
項目あるとし、それらの合計点
がテスト得点として算出されるとする。
である。このとき、オメガは次式で与えられる。
の平均が０、分散が１であることから式（３）は
と書ける。式番号（6.20a）はMcDonald (1999)におけるものである。
さらに、各が他の変数と独立であるとすると、式（6.20a）は
と書ける。
まず、与えられたデータの分析を行うが、不適切な項目が見つかったので、不適切な項目を除いた再分析も行った。
ベイズ分析のためのStanスクリプトを用いる例をリスト１に、リスト１のStanスクリプトを用いるPythonスクリプトをリスト２に示した。ファイルは、まとめて圧縮ファイルestomegafiles.zipとした。リスト２は、CmdStanPyを利用しているが、CmdStanPyの準備の簡単な説明は、このウェブサイトで行った
リスト２のスクリプト（EstOmega.py）は、リスト２のファイル、リスト１のファイル、および入力データファイルを同じフォルダに置き、カレントディレクトリをそのフォルダに移動してから、次のようなコマンドにより実行する。
(stan)
*****/estomegafiles$ python EstOmega.py
上のコマンドを実行すると、出力ファイル名と入力データファイル名の設定が求められる。
(stan) *****/estomegafiles$ python
EstOmega.py
Output file (*.txt) = Results.txt
Input Excel data file (*.xlsx) =
Data.xlsx
上の例では、出力ファイル名としてResults.txtが設定されている。出力ファイル名は、テキストファイル名（*.txt）であれば、任意である。
出力ファイル名を設定すると、入力データファイル名の設定が求められる。
入力データファイルはExcelファイルとして図A１のように用意する。
図A１
１行目に変数名を設定して、次の２行目からデータを設定する。ただし、１列目はデータの識別用値を設定する。１列目の値は、分析では無視される。
入力データファイル名を設定すると、ファイルが読み込まれ、リスト１のStanスクリプトのコンパイルが始まる。コンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、MCMCサンプリングの統計出力が行われ、その後、図A2のグラフが表示される。
図A２
の事後分布のグラフである。
図A２のWindowを閉じると、
の事後分布が表示される（図A３）。
図A３
図A３のWindowを閉じると、
の事後分布のグラフが表示される（図A４）。
図A４
図A４のWindowを閉じると、次に、変数ごとに各変数における因子
の分散の各変数の分散に対する比
の事後分布のグラフが表示される（図A５）。この値は、各変数に占める因子の割合と考えられる。
図A５
図A５のWindowを閉じると、ωの事後分布のグラフが表示される（図A６）。
図A６
グラフのタイトル部に、事後分布の四分位数が表示されている。点推定値は、Gelmanら（2021）に従えば中央値（median）を点推定値とすればよいので、
となる。
図A６のフォームを閉じると、図A７の散布図が表示される。
図A７
横軸に各
の推定値
（事後分布の中央値）、縦軸にテストスコアの和
をとり、プロットしたものである。点
がほぼ直線沿いに散らばっていることが判る。
図A７のWindowを閉じると、スクリプトの実行終了となる。端末画面には、次のメッセージが表示されている。
omega: Q1 = 0.898, Med = 0.904, Q3 = 0.911
Results.txt
was saved.
(stan)
*****/estomegafiles$
出力ファイルResults.txtを開くと、以下のようになっている。
Input data file = Data.xlsx
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]
[ 4 45 51 30 71 63 45 62 58 34 51]
[ 5 74 54 47 60 44 59 66 80 45 54]
[ 6 63 66 58 62 74 31 75 74 58 67]
[ 7 65 65 61 52 62 29 64 60 29 64]
[ 8 72 85 89 82 75 11 76 73 60 82]
[ 9 31 34 27 39 40 34 25 27 36 26]
[10 62 78 45 61 60 60 47 39 63 48]
[11 61 51 22 67 29 73 38 56 62 44]
[12 28 31 28 18 30 37 34 32 47 44]
[13 60 75 60 64 64 57 46 72 75 64]
[14 60 62 63 40 63 58 45 65 64 82]
[15 67 61 66 66 79 62 80 71 43 89]
[16 55 33 45 24 35 29 36 26 37 26]
[17 44 51 40 28 45 35 59 32 38 35]
[18 64 71 40 46 46 26 76 75 21 72]
[19 25 41 27 36 37 71 29 19 42 40]
[20 56 47 41 32 50 42 41 45 50 22]
[21 41 49 78 53 41 32 56 39 28 40]
[22 62 64 49 49 60 74 59 51 47 64]
[23 60 55 65 50 48 32 54 53 67 57]
[24 48 64 63 41 50 19 54 63 53 70]
[25 73 76 51 71 72 65 81 66 15 94]
[26 71 67 59 85 87 43 66 80 41 83]
[27 69 56 48 35 50 43 42 53 42 48]
[28 54 53 47 69 53 47 65 59 35 59]
[29 47 72 68 83 77 58 61 71 51 81]
[30 16 30 32 16 34 47 36 30 45 18]
[31 23 31 36 28 35 91 31 45 26 21]
[32 39 56 25 46 52 28 44 32 21 55]
[33 37 59 53 59 64 52 67 62 29 62]
[34 74 67 59 93 77 47 40 66 17 73]
[35 70 56 71 81 68 57 76 77 61 74]
[36 35 47 27 29 33 50 30 60 26 24]
[37 69 68 72 79 71 39 53 58 41 64]
[38 33 46 33 43 31 24 37 34 38 63]
[39 56 62 51 65 54 34 66 58 26 54]
[40 58 64 70 59 67 54 69 57 56 49]
[41 62 53 52 45 37 43 40 46 17 40]
[42 19 25 21 32 36 65 38 22 51 33]
[43 56 46 34 33 32 58 40 37 52 32]
[44 62 75 85 74 65 58 76 89 63 68]
[45 37 33 20 28 50 26 36 15 55 30]
[46 44 46 59 50 28 62 53 38 28 53]
[47 27 32 56 39 43 55 51 28 61 50]
[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]]
Mean
MCSE
StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ -1531.720000 0.301065 8.339540 ... 1521.21 250.589 1.000770
vF[1]
0.939066 0.006721 0.306042 ... 2643.87 692.703 1.000350
vF[2]
0.129547 0.006144 0.277600 ... 2710.55 682.873 1.001150
vF[3] -0.500504 0.006295 0.285406 ... 2301.82 692.745 1.001450
vF[4]
0.153873 0.005736 0.283963 ... 2812.10 823.800 1.000340
...
...
...
... ... ...
...
...
L2V[8] 0.730895 0.001083 0.067019 ... 2841.38 1438.970 0.999558
L2V[9] 0.034613 0.000526 0.040602 ... 1840.98 1161.000 0.999896
L2V[10] 0.697676 0.001112 0.071477 ... 2969.85 1498.250 1.000020
v1
115.490000 0.060979 4.727360 ... 3303.75 1984.820 0.999984
v2
1414.130000 1.444690 109.619000 ... 3156.40 1940.710 1.001600
[678 rows x 11 columns]
Mu...
Mu[1]: Q1 = 50.839, Med = 51.827, Q3 = 52.809
Mu[2]: Q1 = 53.410, Med = 54.230, Q3 = 55.075
Mu[3]: Q1 = 49.022, Med = 50.182, Q3 = 51.341
Mu[4]: Q1 = 51.669, Med = 52.650, Q3 = 53.684
Mu[5]: Q1 = 51.945, Med = 52.880, Q3 = 53.810
Mu[6]: Q1 = 44.363, Med = 46.127, Q3 = 47.851
Mu[7]: Q1 = 51.609, Med = 52.545, Q3 = 53.394
Mu[8]: Q1 = 51.359, Med = 52.254, Q3 = 53.166
Mu[9]: Q1 = 42.325, Med = 43.729, Q3 = 45.244
Mu[10]: Q1 = 53.096, Med = 54.045, Q3 = 55.008
Sigma...
Sigma[1]: Q1 = 9.662 Med = 10.413, Q3 = 11.267
Sigma[2]: Q1 = 8.234 Med = 8.965, Q3 = 9.722
Sigma[3]: Q1 = 10.973 Med = 11.850, Q3 = 12.780
Sigma[4]: Q1 = 9.872 Med = 10.732, Q3 = 11.640
Sigma[5]: Q1 = 9.246 Med = 10.013, Q3 = 10.845
Sigma[6]: Q1 = 16.298 Med = 17.502, Q3 = 18.792
Sigma[7]: Q1 = 8.617 Med = 9.285, Q3 = 10.003
Sigma[8]: Q1 = 8.574 Med = 9.355, Q3 = 10.189
Sigma[9]: Q1 = 14.553 Med = 15.615, Q3 = 16.776
Sigma[10]: Q1 = 9.413 Med = 10.278, Q3 = 11.176
Lambda...
Lambda[1]: Q1 = 11.977 Med = 12.989, Q3 = 13.988 (Item_1)
Lambda[2]: Q1 = 11.760 Med = 12.623, Q3 = 13.474 (Item_2)
Lambda[3]: Q1 = 11.578 Med = 12.731, Q3 = 13.875 (Item_3)
Lambda[4]: Q1 = 15.075 Med = 16.125, Q3 = 17.156 (Item_4)
Lambda[5]: Q1 = 12.270 Med = 13.290, Q3 = 14.232 (Item_5)
Lambda[6]: Q1 = 0.529 Med = 1.200, Q3 = 2.254 (Item_6)
Lambda[7]: Q1 = 11.370 Med = 12.247, Q3 = 13.189 (Item_7)
Lambda[8]: Q1 = 14.852 Med = 15.725, Q3 = 16.579 (Item_8)
Lambda[9]: Q1 = 1.102 Med = 2.249, Q3 = 3.640 (Item_9)
Lambda[10]: Q1 = 14.818 Med = 15.859, Q3 = 16.855 (Item_10)
Lambda_j**/Var(X_j)...
Lambda[1]**2 / Var(X[1]): Q1 = 0.549, Med = 0.608, Q3 = 0.662
Lambda[2]**2 / Var(X[2]): Q1 = 0.613, Med = 0.665, Q3 = 0.712
Lambda[3]**2 / Var(X[3]): Q1 = 0.473, Med = 0.537, Q3 = 0.592
Lambda[4]**2 / Var(X[4]): Q1 = 0.642, Med = 0.693, Q3 = 0.738
Lambda[5]**2 / Var(X[5]): Q1 = 0.581, 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.581, Med = 0.636, Q3 = 0.685
Lambda[8]**2 / Var(X[8]): Q1 = 0.693, Med = 0.739, Q3 = 0.778
Lambda[9]**2 / Var(X[9]): Q1 = 0.005, Med = 0.020, Q3 = 0.050
Lambda[10]**2 / Var(X[10]): Q1 = 0.655, Med = 0.706, Q3 = 0.750
Q1 = 0.898, Med =
0.904, Q3 = 0.911
まず、入力データが出力されている。
続いて、MCMCサンプリングの統計が出力されている。
事後分布などのグラフで表示されていた情報が、統計の数値データとして読むことができる。
図A４（）、図A５（
）における他とは著しく異なる項目は、上の統計量から６番目と９番目の項目であることが解る。６番目と９番目を除いて再分析を行う。
図A１のデータから、６番目と９番目の項目のデータを除いたものを、図B1に示すように用意した。
図B1 図A１のデータから６番目と９番目のデータを抜いたもの
図B1のデータに対して、リスト２のスクリプトを実行した。
(stan) *****/estomegafiles$ python
EstOmega.py
Output file (*.txt) = Results_69.txt
Input Excel data file (*.xlsx) =
Data_69.xlsx
出力ファイル名をResults_69.txtとした。上の分析の場合と同じ名前でもよいが、同じ名前にすると、ファイルは上書きされて、古いファイルの内容は消される。
入力データファイル名は、図B1のファイルである。
実行すると、の事後分布は、図B2のようである。
図B2
図B2のWindowを閉じると、の事後分布が表示される（図B3）。
図B３
図B3のWindowを閉じると、の事後分布のグラフが表示される（図B4）。
図B4
図B4のWindowを閉じると、の事後分布のグラフが表示される（図B5）。
図B5
図B5のWindowを閉じると、ωの事後分布のグラフが表示される（図B6）。図A６のものより右に寄っていることが分かる。
図B6
図B6のWindowを閉じると、点の散布図が表示される（図B7）。図A７のものより、点が狭い帯内にまとまっていることが分かる。
図B7
図B7のWindowを閉じると、スクリプトの実行終了である。
端末には以下のようなメッセージが表示されている。
omega: Q1 = 0.931, Med = 0.936, Q3 = 0.941
Results_69.txt was saved.
(stan) *****/estomegafiles$
出力ファイルResults_69.txtを開くと、以下のようになっている。
Input data file = Data_69.xlsx
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]
[ 4 45 51 30 71 63 62 58 51]
[ 5 74 54 47 60 44 66 80 54]
[ 6 63 66 58 62 74 75 74 67]
[ 7 65 65 61 52 62 64 60 64]
[ 8 72 85 89 82 75 76 73 82]
[ 9 31 34 27 39 40 25 27 26]
[10 62 78 45 61 60 47 39 48]
[11 61 51 22 67 29 38 56 44]
[12 28 31 28 18 30 34 32 44]
[13 60 75 60 64 64 46 72 64]
[14 60 62 63 40 63 45 65 82]
[15 67 61 66 66 79 80 71 89]
[16 55 33 45 24 35 36 26 26]
[17 44 51 40 28 45 59 32 35]
[18 64 71 40 46 46 76 75 72]
[19 25 41 27 36 37 29 19 40]
[20 56 47 41 32 50 41 45 22]
[21 41 49 78 53 41 56 39 40]
[22 62 64 49 49 60 59 51 64]
[23 60 55 65 50 48 54 53 57]
[24 48 64 63 41 50 54 63 70]
[25 73 76 51 71 72 81 66 94]
[26 71 67 59 85 87 66 80 83]
[27 69 56 48 35 50 42 53 48]
[28 54 53 47 69 53 65 59 59]
[29 47 72 68 83 77 61 71 81]
[30 16 30 32 16 34 36 30 18]
[31 23 31 36 28 35 31 45 21]
[32 39 56 25 46 52 44 32 55]
[33 37 59 53 59 64 67 62 62]
[34 74 67 59 93 77 40 66 73]
[35 70 56 71 81 68 76 77 74]
[36 35 47 27 29 33 30 60 24]
[37 69 68 72 79 71 53 58 64]
[38 33 46 33 43 31 37 34 63]
[39 56 62 51 65 54 66 58 54]
[40 58 64 70 59 67 69 57 49]
[41 62 53 52 45 37 40 46 40]
[42 19 25 21 32 36 38 22 33]
[43 56 46 34 33 32 40 37 32]
[44 62 75 85 74 65 76 89 68]
[45 37 33 20 28 50 36 15 30]
[46 44 46 59 50 28 53 38 53]
[47 27 32 56 39 43 51 28 50]
[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]]
Mean
MCSE
StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ -1188.020000 0.316133 8.132270 ... 1198.04 245.409 1.00504
vF[1] 0.931394 0.006707 0.303334 ... 2479.93 774.468 1.00027
vF[2] 0.133307 0.006507 0.283712 ... 2312.85 734.531 1.00097
vF[3] -0.514093 0.006231 0.280542 ... 2793.45 776.290 1.00060
vF[4] 0.160360 0.006567 0.287864 ... 2258.66 736.626 1.00074
...
...
...
... ... ...
...
...
L2V[6] 0.630956 0.001414 0.080581 ... 2934.79 1405.170 1.00063
L2V[7] 0.729854 0.001136 0.066082 ... 2898.76 1309.020 1.00135
L2V[8] 0.698674 0.001262 0.072172 ... 2638.18 1336.090 1.00028
v1
111.462000 0.060987 4.223700 ... 3378.10 1843.830 1.00211
v2
848.904000 1.146130 71.037400 ... 3049.17 1445.850 1.00066
[564 rows x 11 columns]
Mu...
Mu[1]: Q1 = 50.789, Med = 51.797, Q3 = 52.785
Mu[2]: Q1 = 53.426, Med = 54.263, Q3 = 55.100
Mu[3]: Q1 = 49.061, Med = 50.175, Q3 = 51.338
Mu[4]: Q1 = 51.633, Med = 52.634, Q3 = 53.634
Mu[5]: Q1 = 51.957, Med = 52.907, Q3 = 53.889
Mu[6]: Q1 = 51.613, Med = 52.475, Q3 = 53.338
Mu[7]: Q1 = 51.422, Med = 52.256, Q3 = 53.105
Mu[8]: Q1 = 53.094, Med = 54.016, Q3 = 54.970
Sigma...
Sigma[1]: Q1 = 9.633 Med = 10.400, Q3 = 11.236
Sigma[2]: Q1 = 8.242 Med = 8.957, Q3 = 9.694
Sigma[3]: Q1 = 11.022 Med = 11.863, Q3 = 12.800
Sigma[4]: Q1 = 9.847 Med = 10.686, Q3 = 11.612
Sigma[5]: Q1 = 9.274 Med = 10.024, Q3 = 10.868
Sigma[6]: Q1 = 8.568 Med = 9.242, Q3 = 9.980
Sigma[7]: Q1 = 8.626 Med = 9.399, Q3 = 10.249
Sigma[8]: Q1 = 9.368 Med = 10.211, Q3 = 11.103
Lambda...
Lambda[1]: Q1 = 12.015 Med = 12.991, Q3 = 14.034 (Item_1)
Lambda[2]: Q1 = 11.769 Med = 12.651, Q3 = 13.508 (Item_2)
Lambda[3]: Q1 = 11.545 Med = 12.702, Q3 = 13.881 (Item_3)
Lambda[4]: Q1 = 15.021 Med = 16.112, Q3 = 17.138 (Item_4)
Lambda[5]: Q1 = 12.249 Med = 13.271, Q3 = 14.253 (Item_5)
Lambda[6]: Q1 = 11.374 Med = 12.282, Q3 = 13.174 (Item_7)
Lambda[7]: Q1 = 14.786 Med = 15.757, Q3 = 16.682 (Item_8)
Lambda[8]: Q1 = 14.787 Med = 15.816, Q3 = 16.862 (Item_10)
Lambda_j**/Var(X_j)...
Lambda[1]**2 / Var(X[1]): Q1 = 0.551, Med = 0.610, Q3 = 0.662
Lambda[2]**2 / Var(X[2]): Q1 = 0.613, Med = 0.666, Q3 = 0.714
Lambda[3]**2 / Var(X[3]): Q1 = 0.470, Med = 0.534, Q3 = 0.591
Lambda[4]**2 / Var(X[4]): Q1 = 0.642, Med = 0.695, Q3 = 0.740
Lambda[5]**2 / Var(X[5]): Q1 = 0.576, Med = 0.636, Q3 = 0.687
Lambda[6]**2 / Var(X[6]): Q1 = 0.583, Med = 0.639, Q3 = 0.689
Lambda[7]**2 / Var(X[7]): Q1 = 0.691, Med = 0.737, Q3 = 0.777
Lambda[8]**2 / Var(X[8]): Q1 = 0.657, Med = 0.706, Q3 = 0.750
omega...
Q1 = 0.931, Med =
0.936, Q3 = 0.941
他と特に互いに異なる値の項目は認められない。ωの値も高くなっており
１０項目のときの
omega...
Q1 = 0.898, Med =
0.904, Q3 = 0.911
から、第６項目と第９項目を抜いた８項目でのωは
omega...
Q1 = 0.931, Med =
0.936, Q3 = 0.941
であり、四分位数の比較でも８項目に減らしたときのωの方が高い値であると言えよう。
リスト１ MCMCサンプリングのためのStanスクリプト（est_omega.stan）
//
//
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);
}
リスト２ リスト１のStanスクリプトを用いるPythonスクリプト（EstOmega.py）
from cmdstanpy import CmdStanModel
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
foutnm = input("Output file
(*.txt) = ")
fout = open(foutnm, 'w')
#
# Excelデータファイルの読み込み
#
s = input("Input Excel data file
(*.xlsx) = ")
xlsx_fl = pd.ExcelFile(s)
data_xlsx = pd.read_excel(xlsx_fl)
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}
model =
CmdStanModel(stan_file="est_omega.stan") # コンパイル
fit = model.sample(data=Data)
# サンプリング
print(fit.summary())
fout.write('\n{}'.format(fit.summary()))
df = fit.draws_pd()
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')
mus = np.mean(X, axis=0)
stds = np.std(X, axis=0)
#
# 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(foutnm, '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.
岡本安晴（2014）、心理学データ分析と測定．勁草書房