Up

心理尺度の信頼性係数ω

ベイズ分析

 

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

 

 

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

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

 

 

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

 

 

の平均が0、分散が1であることから式(3)は

 

 

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

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

 

 

と書ける。

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

 

 

 

ベイズ分析例

 

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

 

Pystan3との組み合わせでArvizが使えるようになったので、更新版をこのウェブサイトに用意した(2025.03.22)。MCMCサンプルのRhat値とトレース図が得られるようにした。

 

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

リスト1のStanスクリプトは、式(1)

 

 

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

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

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

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

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

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

 

 

 

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

Input Excel data file (*.xlsx) =

 

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

図1

 

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

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

 

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は、エディタで開くことができる。

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

 

リスト2のスクリプトを実行すると、まずの事後分布のグラフが表示される(図2)。

図2

 

図2のWindowを閉じると、の事後分布が表示される(図3)。

図3

 

図3のWindowを閉じると、の事後分布のグラフが表示される(図4)。

図4

 

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

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

図5

 

図5のWindowを閉じると、ωの事後分布のグラフが表示される(図6)。

図6

 

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

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

図7

 

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

図7の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.830,  Med = 51.831,  Q3 = 52.865

Mu[2]:  Q1 = 53.353,  Med = 54.220,  Q3 = 55.137

Mu[3]:  Q1 = 49.028,  Med = 50.159,  Q3 = 51.319

Mu[4]:  Q1 = 51.609,  Med = 52.646,  Q3 = 53.687

Mu[5]:  Q1 = 51.915,  Med = 52.871,  Q3 = 53.848

Mu[6]:  Q1 = 44.512,  Med = 46.133,  Q3 = 47.731

Mu[7]:  Q1 = 51.613,  Med = 52.482,  Q3 = 53.394

Mu[8]:  Q1 = 51.407,  Med = 52.248,  Q3 = 53.139

Mu[9]:  Q1 = 42.293,  Med = 43.748,  Q3 = 45.241

Mu[10]:  Q1 = 53.102,  Med = 54.068,  Q3 = 55.051

 

Sigma...

Sigma[1]:  Q1 = 9.642  Med = 10.410,  Q3 = 11.247

Sigma[2]:  Q1 = 8.284  Med = 8.978,  Q3 = 9.735

Sigma[3]:  Q1 = 10.983  Med = 11.842,  Q3 = 12.747

Sigma[4]:  Q1 = 9.919  Med = 10.687,  Q3 = 11.617

Sigma[5]:  Q1 = 9.202  Med = 9.970,  Q3 = 10.813

Sigma[6]:  Q1 = 16.324  Med = 17.401,  Q3 = 18.750

Sigma[7]:  Q1 = 8.567  Med = 9.282,  Q3 = 10.048

Sigma[8]:  Q1 = 8.575  Med = 9.357,  Q3 = 10.224

Sigma[9]:  Q1 = 14.564  Med = 15.605,  Q3 = 16.794

Sigma[10]:  Q1 = 9.443  Med = 10.207,  Q3 = 11.154

 

Lambda...

Lambda[1]:  Q1 = 12.019  Med = 13.024,  Q3 = 14.026  (Item_1)

Lambda[2]:  Q1 = 11.692  Med = 12.602,  Q3 = 13.466  (Item_2)

Lambda[3]:  Q1 = 11.536  Med = 12.743,  Q3 = 13.908  (Item_3)

Lambda[4]:  Q1 = 15.121  Med = 16.120,  Q3 = 17.180  (Item_4)

Lambda[5]:  Q1 = 12.310  Med = 13.275,  Q3 = 14.242  (Item_5)

Lambda[6]:  Q1 = 0.537  Med = 1.211,  Q3 = 2.216  (Item_6)

Lambda[7]:  Q1 = 11.363  Med = 12.291,  Q3 = 13.136  (Item_7)

Lambda[8]:  Q1 = 14.754  Med = 15.739,  Q3 = 16.642  (Item_8)

Lambda[9]:  Q1 = 1.191  Med = 2.267,  Q3 = 3.587  (Item_9)

Lambda[10]:  Q1 = 14.841  Med = 15.834,  Q3 = 16.835  (Item_10)

 

Lambda_j**/Var(X_j)...

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

Lambda[2]**2 / Var(X[2]):  Q1 = 0.610,  Med = 0.662,  Q3 = 0.710

Lambda[3]**2 / Var(X[3]):  Q1 = 0.474,  Med = 0.538,  Q3 = 0.594

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

Lambda[5]**2 / Var(X[5]):  Q1 = 0.584,  Med = 0.640,  Q3 = 0.688

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

Lambda[7]**2 / Var(X[7]):  Q1 = 0.580,  Med = 0.636,  Q3 = 0.688

Lambda[8]**2 / Var(X[8]):  Q1 = 0.691,  Med = 0.739,  Q3 = 0.778

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

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

 

omega...

         Q1 = 0.898,  Med = 0.905,  Q3 = 0.911

 

 

 

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

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

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

図4()、図5()における他とは著しく異なる項目は、上の統計量から6番目と9番目の項目であることが解る。6番目と9番目を除いて再分析を行う。

 

 

 

再分析

 

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

A1  図1のデータから6番目と9番目のデータを抜いたもの

 

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

 

 

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

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

 

 

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

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

A2

 

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

A3

 

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

A4

 

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

A5

 

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

A6

 

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

A7

 

出力ファイル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.764,  Med = 51.760,  Q3 = 52.738

Mu[2]:  Q1 = 53.410,  Med = 54.258,  Q3 = 55.118

Mu[3]:  Q1 = 49.034,  Med = 50.143,  Q3 = 51.341

Mu[4]:  Q1 = 51.568,  Med = 52.626,  Q3 = 53.648

Mu[5]:  Q1 = 51.974,  Med = 52.900,  Q3 = 53.838

Mu[6]:  Q1 = 51.618,  Med = 52.488,  Q3 = 53.403

Mu[7]:  Q1 = 51.372,  Med = 52.248,  Q3 = 53.125

Mu[8]:  Q1 = 53.082,  Med = 54.033,  Q3 = 54.994

 

Sigma...

Sigma[1]:  Q1 = 9.636  Med = 10.431,  Q3 = 11.266

Sigma[2]:  Q1 = 8.233  Med = 8.907,  Q3 = 9.690

Sigma[3]:  Q1 = 11.030  Med = 11.894,  Q3 = 12.825

Sigma[4]:  Q1 = 9.887  Med = 10.713,  Q3 = 11.659

Sigma[5]:  Q1 = 9.228  Med = 9.968,  Q3 = 10.830

Sigma[6]:  Q1 = 8.638  Med = 9.297,  Q3 = 10.042

Sigma[7]:  Q1 = 8.619  Med = 9.381,  Q3 = 10.245

Sigma[8]:  Q1 = 9.387  Med = 10.177,  Q3 = 11.057

 

Lambda...

Lambda[1]:  Q1 = 11.991  Med = 13.022,  Q3 = 14.071  (Item_1)

Lambda[2]:  Q1 = 11.698  Med = 12.601,  Q3 = 13.507  (Item_2)

Lambda[3]:  Q1 = 11.491  Med = 12.644,  Q3 = 13.867  (Item_3)

Lambda[4]:  Q1 = 15.017  Med = 16.109,  Q3 = 17.154  (Item_4)

Lambda[5]:  Q1 = 12.293  Med = 13.268,  Q3 = 14.273  (Item_5)

Lambda[6]:  Q1 = 11.346  Med = 12.275,  Q3 = 13.163  (Item_7)

Lambda[7]:  Q1 = 14.737  Med = 15.733,  Q3 = 16.612  (Item_8)

Lambda[8]:  Q1 = 14.813  Med = 15.798,  Q3 = 16.847  (Item_10)

 

Lambda_j**/Var(X_j)...

Lambda[1]**2 / Var(X[1]):  Q1 = 0.551,  Med = 0.610,  Q3 = 0.663

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

Lambda[3]**2 / Var(X[3]):  Q1 = 0.468,  Med = 0.530,  Q3 = 0.589

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

Lambda[5]**2 / Var(X[5]):  Q1 = 0.584,  Med = 0.640,  Q3 = 0.690

Lambda[6]**2 / Var(X[6]):  Q1 = 0.577,  Med = 0.637,  Q3 = 0.685

Lambda[7]**2 / Var(X[7]):  Q1 = 0.689,  Med = 0.737,  Q3 = 0.778

Lambda[8]**2 / Var(X[8]):  Q1 = 0.658,  Med = 0.706,  Q3 = 0.750

 

omega...

         Q1 = 0.931,  Med = 0.936,  Q3 = 0.941

 

 

 

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

10項目のときの

 

         Q1 = 0.898,  Med = 0.905,  Q3 = 0.911

 

 

から、第6項目と第9項目を抜いた8項目でのωは

 

         Q1 = 0.931,  Med = 0.936,  Q3 = 0.941

 

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

 

 

 

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

 

import numpy as np

import pandas as pd

import stan

import pickle

 

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')

 

             

fout.close()

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

 

 

 

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

 

import numpy as np

import pandas as pd

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

 

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)

 

with open('fit.pkl', 'rb') as f:     #  保存されたMCMCサンプルの読み込み

    df = pickle.load(f)

#

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

#

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()

#

#        Sigma4分位数の計算

#

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