Up

マグニチュード推定法

ベイズ分析(CmdStanPy 対応)

 

CmdStanPyの準備の簡単な説明のウェブサイトを用意した。

 

プロセティックな感覚(量的な感覚)の場合、物理刺激量と感覚量の間の関係を表す式として、べき法則と呼ばれている次式

 

 

あるいは、これをより一般化した次式

 

 

がある。ただし、式番号は<岡本安晴(2014)「心理学データ分析と測定」勁草書房>におけるものである。

心理物理学実験において、第試行で提示された物理刺激の物理量をで、その感覚量のマグニチュード推定量(Magnitude Estimation: ME)をで表すとき、以下のモデルを設定する(Stevens, 1975; p. 169)。

モデル(10.2.1)に対しては、

 

 

モデル(10.2.2)に対しては、

 

 

ここで、は平均、分散の正規分布を表す。

上記の確率モデルをデータに適用して、ベイズ分析によりパラメータ値を求める。ベイズ分析については、<岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版>の題3部、第9章〜第15章、で説明している。

 

まず、一般化されたモデル(10.2.4)による分析について説明する。その後、モデル(10.2.4)においてと置いた場合のモデル(10.2.3)による分析を示す。

なお、ファイルはmefiles.zipにまとめている。自由に利用されたい。

 

 

 

一般化されたモデル(10.2.4)の場合

 

 

モデル(10.2.4)に基づいてベイズ分析を行うためのStanスクリプトをリストG1のように用意した。

 

モデル(10.2.2)では、3つのパラメータの間に強い連関が予想される。パラメータが小さくなったときは、パラメータの値を大きく、あるいはパラメータの値を負の方向に大きくすればモデル(10.2.2)とデータの関係の維持が期待される。このパラメータ間の連関は推定の不安定性として現れる。推定の不安定さを抑える方法として、事前分布による制約が考えられる。例えば、パラメータが絶対閾値を反映するように事前分布を設定することが考えられる。リストGでは、パラメータの範囲を、0以上、提示刺激値の最小値以下としている。パラメータの推定が不安定なときに事前分布による制約を用いる方法は、Wichmann & Hill(2001により推奨されている。,

 

リストG1Stanスクリプトを用いるPythonスクリプトをリストG2に示す。リストG1およびリストG2のファイルなどは、圧縮ファイルmefiles.zipにまとめた。ダウンロード展開して自由に利用されたい。

 

リストG1のファイル(MEc.stan)、リストG2のファイル(ME_cBysn.py)と入力データファイルを同じフォルダ内に置き、そのフォルダをカレントディレクトリとする。CmdStanPyのインストールされた環境で次のコマンドを実行すればよい。

 

(stan) *****/mefiles$ python ME_cBysn.py

 

上のコマンドを実行すると、入力データファイル名と出力ファイル名の設定が求められる。次の例では、入力データファイルとしてData.xlsxが、出力ファイルとしてResults.txtが設定されている。

 

(stan) *****/mefiles$ python ME_cBysn.py

Input Data File (*.xlsx) = Data.xlsx

Output Text File (*.txt) = Results.txt

 

出力ファイルは適当なテキストファイル名を設定すればよい。

入力データファイルは、Excelファイルとして図G1の形式で用意する。

 

アプリケーション, テーブル

AI 生成コンテンツは誤りを含む可能性があります。

テーブル

AI 生成コンテンツは誤りを含む可能性があります。

G1

 

1行目は、データの変数名を入力する。データは2行目から入力する。1列目はデータの識別用で、図1では試行が入力されている。この第1列目は、データの識別用であり、数値でも文字でも何でもよい。また、識別用であるが、同じ値を繰り返して用いてもよい。データの分析においては無視される。2列目に提示刺激の物理量、3列目にマグニチュード推定値を入力する。

出力ファイル名を設定すると、入力データが読み込まれ、リストG1Stanスクリプトのコンパイルが始まる。このコンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが行われ、サンプリングが終了すると、トレース図が表示される(図G2)。

 

G2

 

G2のフォームを閉じると、パラメータの3次元散布図が表示される(図G3)。

 

G3

 

G3の3D散布図は、図上をマウスでドラッグすると回転することができる(図G4)。

 

G4

 

G4Windowを閉じると、図G5のグラフが表示される。

 

G5

 

パラメータabcの事後分布のグラフである。

G5Windowを閉じると、図G6のグラフが表示される。

 

G6

 

パラメータ対の散布図である。

G6Windowを閉じると、図G7のグラフが表示される。

 

G7

 

左側のグラフでは、データの点と、事後分布の中央値をパラメータの推定値として式(10.2.2)のグラフが描かれている。Gelmanら(2021は、事後分布の代表値として、平均値より安定した中央値を勧めている。

右側のグラフでは、MCMCサンプリングの結果から500個のサンプル値をランダムに選び、各サンプル値に対して点線でグラフ

 

 

が描かれている。第1四分位数Q!と第3四分位数Q3に対応する値が水色とライムの破線で描かれている。

中央の赤線は、各パラメータの中央値を推定値としてグラフ(10.2.2)を描いたものである。

 

G7のフォームを閉じると、スクリプトの実行終了である。

 

 

 

リストG1 モデル(10.2.4)のStanスクリプト(MEc.stan

 

data {

    int N;

    array[N] real Sts;

    array[N] real MEs;

}

transformed data {

    array[N] real LogMEs;

    real minSt;

    for (i in 1:N) {

        LogMEs[i] = log(MEs[i]);

    }

    minSt = Sts[1];

    for (i in 2:N) {

        if (minSt > Sts[i]) {

            minSt = Sts[i];

        }

    }

}

parameters {

    real<lower = 0.0> a;

    real<lower = 0.0, upper = minSt> c;

    real<lower = 0.0> b;

    real<lower = 0.0> sgm;

}

model {

    a ~ exponential(0.001);

    b ~ exponential(0.001);

    c ~ uniform(0.0, minSt);

    sgm ~ exponential(0.001);

    for (i in 1:N) {

        LogMEs[i] ~ normal(log(a * pow(Sts[i] - c, b)), sgm);   //   Eq. (10.2.4)

    }

}

 

 

 

リストG2 リストG1を利用するPythonスクリプト(ME_cBysn.py

 

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import arviz as az

import matplotlib as mplt

import matplotlib.pyplot as plt

import seaborn as sb

 

finnm = input('Input Data File (*.xlsx) = ')

 

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl,  header = None)

data = data_xlsx.values

 

foutnm = input('Output Text File (*.txt) = ')

fout = open(foutnm, 'w')

fout.write('Input Data File = {}\n'.format(finnm))

 

N = len(data) - 1

Sts = [v for v in data[1:, 1]]

MEs = [v for v in data[1:, 2]]

print('N =', N)

print(len(Sts))

print(len(MEs))

 

for st, me in zip(Sts, MEs):

    print(st, me)

 

Data = {'N':N, 'Sts':Sts, 'MEs':MEs}

 

model = CmdStanModel(stan_file="MEc.stan")    #   コンパイル

fit = model.sample(data=Data)                 #   サンプリング

 

print(fit.summary())

fout.write('\n{}'.format(fit.summary()))

 

inf_data = az.from_cmdstanpy(fit)           #   Arviz(az)用のデータに変換

az.plot_trace(inf_data, figsize=(10,7))     #   トレースプロット

plt.tight_layout()

plt.show()

 

dframe = fit.draws_pd()                     #   Pandas DataFrameへ変換

 

fig = plt.figure(figsize = (7,7))

ax = fig.add_subplot(111, projection = '3d')

ax.scatter(dframe['a'], dframe['b'], dframe['c'])

ax.set_xlabel('a', fontsize = 16)

ax.set_ylabel('b', fontsize = 16)

ax.set_zlabel('c', fontsize = 16)

plt.title('Scatter plot of (a,b,c)', fontsize = 20)

plt.show()

 

plt.figure(figsize=(10,8))

 

plt.subplot(2,4,(1,2))

 

mean_a = dframe['a'].mean()

P025_a, P50_a, P975_a = np.percentile(dframe['a'], [2.5, 50, 97.5])

 

plt.title('Posterior Distribution of a', fontsize = 20)

sb.kdeplot(dframe['a'], clip = (0.0, max(dframe['a'])),

           label ='Mean = {0:.2f}\nMed = {1:.2f}\n95%CI = [{2:.2f}, {3:.2f}]'.

           format(mean_a, P50_a, P025_a, P975_a))

plt.xlabel('a', fontsize = 18)

plt.ylabel('Density', fontsize = 14)

plt.legend(fontsize = 12)

plt.tight_layout()

 

plt.subplot(2,4,(3,4))

 

mean_b = dframe['b'].mean()

P025_b, P50_b, P975_b = np.percentile(dframe['b'], [2.5, 50, 97.5])

 

plt.title('Posterior Distribution of b', fontsize = 20)

sb.kdeplot(dframe['b'], label = 'Mean = {0:.2f}\nMed = {1:.2f}\n95%CI = [{2:.2f}, {3:.2f}]'.

           format(mean_b, P50_b, P025_b, P975_b))

plt.xlabel('b', fontsize = 18)

plt.ylabel('Density', fontsize = 14)

plt.legend(fontsize = 12)

plt.tight_layout()

 

plt.subplot(2,4,(6,7))

 

mean_c = dframe['c'].mean()

P025_c, P50_c, P975_c = np.percentile(dframe['c'], [2.5, 50, 97.5])

 

plt.title('Posterior Distribution of c', fontsize = 20)

sb.kdeplot(dframe['c'], label ='Mean = {0:.2f}\nMed = {1:.2f}\n95%CI = [{2:.2f}, {3:.2f}]'.

           format(mean_c, P50_c, P025_c, P975_c))

plt.xlabel('c', fontsize = 18)

plt.ylabel('Density', fontsize = 14)

plt.legend(fontsize = 12)

plt.tight_layout()

plt.show()

 

plt.figure(figsize=(10, 8))

 

plt.subplot(2,4,(1,2))

 

plt.title('Scattergram of (a,b)', fontsize = 20)

plt.plot(dframe['a'], dframe['b'], 'bo')

plt.xlabel('a', fontsize = 18)

plt.ylabel('b', fontsize = 18)

plt.tight_layout()

 

plt.subplot(2,4,(3,4))

plt.title('Scattergram of (a,c)', fontsize = 20)

plt.plot(dframe['a'], dframe['c'], 'bo')

plt.xlabel('a', fontsize = 18)

plt.ylabel('c', fontsize = 18)

plt.tight_layout()

 

plt.subplot(2,4,(6,7))

plt.title('Scattergram of (b,c)', fontsize = 20)

plt.plot(dframe['b'], dframe['c'], 'bo')

plt.xlabel('b', fontsize = 18)

plt.ylabel('c', fontsize = 18)

plt.tight_layout()

plt.show()

 

min_St = min(Sts)

max_St = max(Sts)

St_Values = np.arange(min_St, max_St, (max_St - min_St) / 1000.0)

 

def PowerLaw(x, a, b, c):

    return a * ((x - c) ** b)

 

plt.figure(figsize=(12,5))

 

plt.subplot(1,2,1)

plt.plot(Sts, MEs, 'bo', alpha = 0.5, label = 'Data')

plt.plot(St_Values, [PowerLaw(v, P50_a, P50_b, P50_c) for v in St_Values], 'k-',

         label = r'$\psi=a(\varphi-c)^b$' + ('\na = {0:.3f}'.format(P50_a)) +

                     ('\nb = {0:.3f}'.format(P50_b)) +

                      ('\nc = {0:.3f}'.format(P50_c)))

plt.xlabel(r'$\varphi$', fontsize = 18)

plt.ylabel(r'$\psi$', fontsize = 18)

plt.legend(fontsize = 12)

plt.tight_layout()

 

plt.subplot(1,2,2)

 

St_grid = np.linspace(min_St, max_St, 100)

tn_samples = len(dframe)

n_samples = 500

idx_samples = np.random.choice(range(tn_samples), n_samples, replace = False)

curves_smpl = []

for s in idx_samples:

    va = dframe['a'][s]

    vb = dframe['b'][s]

    vc = dframe['c'][s]

    phis = PowerLaw(St_grid, va, vb, vc)

    plt.plot(St_grid, phis, ls = ':', lw = 1, alpha = 0.5)

    curves_smpl.append(phis)

curves_smpl = np.array(curves_smpl)

c_s_q1, c_s_q3 = np.percentile(curves_smpl, [25, 75], axis = 0)

plt.plot(St_grid, c_s_q1, c = 'aqua', ls = '--', label = 'Q1')

plt.plot(St_grid, c_s_q3, c = 'lime', ls = '--', label = 'Q3')

plt.plot(St_Values, [PowerLaw(v, P50_a, P50_b, P50_c) for v in St_Values],

         'r-', lw=3, label = r'$\psi=a(\varphi-c)^b$' + ('\na = {0:.3f}'.format\

                        (P50_a)) + ('\nb = {0:.3f}'.format(P50_b)) +

                      ('\nc = {0:.3f}'.format(P50_c)))

plt.plot(Sts, MEs, 'o', c = 'b', alpha = 0.5, label = 'Data')

plt.xlabel(r'$\varphi$', fontsize = 18)

plt.ylabel(r'$\psi$', fontsize = 18)

plt.legend(fontsize = 12)

plt.tight_layout()

plt.show()

 

fout.close()

print()

print(foutnm + ' was saved.')

print()

 

 

 

モデルの場合

 

確率モデル

 

 

StanスクリプトをリストS1のように用意した。

リストS1を用いるPythonスクリプトは、リストS2のように用意した。使用法は、リストG2と同じである。スクリプトファイルなどは、まとめて圧縮ファイルmefiles.zipとした。ダウンロード展開すれば自由に利用できる。

リストS1のファイル、リストS2のファイル、および入力データファイルを同じフォルダに置いてカレントディレクトリをそのフォルダに移動して、CmdStanPyのインストールされた環境において、次のようなコマンドを実行すればよい。CmdStanPyの準備の簡単な説明は、このウェブサイトで行った。

 

(stan) *****/mefiles$ python MEBysn.py

 

上のコマンドを実行すると、入力データファイル名の設定が求められる。

 

(stan) *****/mefiles$ python MEBysn.py

Input Data File (*.xlsx) = Data.xlsx

Output Text File (*.txt) = Results.txt

 

入力データファイルは、Excelファイルとして図Gに示す様式で準備する。入力データファイル名を入力すると、次に出力ファイル名の設定が求められる。出力ファイルは、テキストファイル(ファイル拡張子が「.txt」)であれば任意の名前でよい。

 

出力ファイル名を設定すると、データが読み込まれ、Stanスクリプトがコンパイルされる。コンパイルには多少時間が掛かる。コンパイル後、MCMCサンプリングが実行され、サンプリング、MCMCサンプリングのトレース図(図S1)が表示される。

 

S1

 

S1Windowを閉じると、図S2のグラフが表示される。

 

S2

 

パラメータabの事後分布と、abの散布図である。

S2Windowを閉じると、図S3のグラフが表示される。

 

S3

 

左側のグラフは、事後分布の中央値を点推定値として次式モデル

 

 

を描いたものである。中央値を点推定値とすることは、Gelmanら(2021が平均値より安定しているとして勧めている。

右のグラフは、MCMCサンプリングの結果から500個のサンプル値をランダムに選び、各サンプル値に対して点線でグラフを描いたものである。第1四分位数Q!と第3四分位数Q3に対応する値が水色とライムの破線で描かれている。中央の赤線は、各パラメータの中央値を推定値としてグラフを描いたものである。

 

S3Windowを閉じると、スクリプトの実行終了である。

端末には、以下のようなメッセージが表示されている。

 

 

Results.txt was saved.

 

(stan) *****/mefiles$

 

出力ファイルResults.txtを開くと、以下のような内容である。

 

Input Data File = Data.xlsx

 

           Mean      MCSE    StdDev  ...  ESS_tail  ESS_bulk/s    R_hat

lp__  28.077800  0.043047  1.348590  ...   1266.14     1898.61  1.00371

a      0.093302  0.001147  0.034230  ...   1313.06     1480.39  1.00503

b      0.615792  0.001877  0.055588  ...   1281.25     1445.08  1.00501

sgm    0.324856  0.000961  0.035652  ...   1559.21     2333.04  1.00312

 

[4 rows x 11 columns]

 

 

 

 

リストS1 モデル(10.2.3)のStanスクリプト(ME.stan

 

data {

    int N;

    array[N] real Sts;

    array[N] real MEs;

}

parameters {

    real<lower = 0.0> a;

    real<lower = 0.0> b;

    real<lower = 0.0> sgm;

}

model {

    a ~ exponential(0.001);

    b ~ exponential(0.001);

    sgm ~ exponential(0.001);

    for (i in 1:N) {

        log(MEs[i]) ~ normal(log(a * pow(Sts[i], b)), sgm);   //   Eq. (10.2.3)

    }

}

 

 

 

リストS2 リストS1Stanスクリプトを用いるPythonスクリプト(MEBysn.py

 

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import arviz as az

import matplotlib as mplt

import matplotlib.pyplot as plt

import seaborn as sb

 

finnm = input('Input Data File (*.xlsx) = ')

 

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl,  header = None)

data = data_xlsx.values

 

foutnm = input('Output Text File (*.txt) = ')

fout = open(foutnm, 'w')

fout.write('Input Data File = {}\n'.format(finnm))

 

N = len(data) - 1

Sts = [v for v in data[1:, 1]]

MEs = [v for v in data[1:, 2]]

print('N =', N)

print(len(Sts))

print(len(MEs))

 

for st, me in zip(Sts, MEs):

    print(st, me)

 

Data = {'N':N, 'Sts':Sts, 'MEs':MEs}

 

model = CmdStanModel(stan_file="ME.stan")    #   コンパイル

fit = model.sample(data=Data)                #   サンうリング

 

print(fit.summary())

fout.write('\n{}'.format(fit.summary()))

 

inf_data = az.from_cmdstanpy(fit)           #   Arviz(az)用のデータに変換

az.plot_trace(inf_data, figsize=(10,7))     #   トレースプロット

plt.tight_layout()

plt.show()

 

dframe = fit.draws_pd()                     #   Pandas DataFrameへ変換

 

plt.figure(figsize=(10,7))

 

mean_a = dframe['a'].mean()

P025_a, P50_a, P975_a = np.percentile(dframe['a'], [2.5, 50, 97.5])

 

plt.subplot(2,4,(1,2))

 

plt.title('Posterior Distribution of a', fontsize = 20)

sb.kdeplot(dframe['a'], clip = (0.0, max(dframe['a'])),

           label ='Mean = {0:.2f}\nMed = {1:.2f}\n95%CI = [{2:.2f}, {3:.2f}]'.

           format(mean_a, P50_a, P025_a, P975_a))

plt.xlabel('a', fontsize = 18)

plt.ylabel('Density', fontsize = 14)

plt.legend(fontsize = 12)

plt.tight_layout()

 

mean_b = dframe['b'].mean()

P025_b, P50_b, P975_b = np.percentile(dframe['b'], [2.5, 50, 97.5])

 

plt.subplot(2,4,(3,4))

 

plt.title('Posterior Distribution of b', fontsize = 20)

sb.kdeplot(dframe['b'], label = 'Mean = {0:.2f}\nMed = {1:.2f}\n95%CI = [{2:.2f}, {3:.2f}]'.

           format(mean_b, P50_b, P025_b, P975_b))

plt.xlabel('b', fontsize = 18)

plt.ylabel('Density', fontsize = 14)

plt.legend(fontsize = 12)

plt.tight_layout()

 

plt.subplot(2,4,(6,7))

 

plt.title('Scattergram of (a,b)', fontsize = 20)

plt.plot(dframe['a'], dframe['b'], 'bo')

plt.xlabel('a', fontsize = 18)

plt.ylabel('b', fontsize = 18)

plt.tight_layout()

plt.show()

 

min_St = min(Sts)

max_St = max(Sts)

St_Values = np.arange(min_St, max_St, (max_St - min_St) / 1000.0)

 

def PowerLaw(x, a, b):

    return a * (x ** b)

 

plt.figure(figsize=(10,5))

 

plt.subplot(1,2,1)

 

plt.plot(Sts, MEs, 'bo', alpha = 0.5, label = 'Data')

plt.plot(St_Values, [PowerLaw(v, P50_a, P50_b) for v in St_Values], 'k-',

         label = r'$\psi=a\varphi^b$' + ('\na = {0:.3f}'.format(P50_a)) +

                     ('\nb = {0:.3f}'.format(P50_b)))

plt.xlabel(r'$\varphi$', fontsize = 18)

plt.ylabel(r'$\psi$', fontsize = 18)

plt.legend(fontsize = 12)

plt.tight_layout()

 

plt.subplot(1,2,2)

 

St_grid = np.linspace(min_St, max_St, 100)

tn_samples = len(dframe)

n_samples = 500

idx_samples = np.random.choice(range(tn_samples), n_samples, replace = False)

curves_smpl = []

for s in idx_samples:

    va = dframe['a'][s]

    vb = dframe['b'][s]

    phis = PowerLaw(St_grid, va, vb)

    plt.plot(St_grid, phis, ls = ':', lw = 1, alpha = 0.5)

    curves_smpl.append(phis)

curves_smpl = np.array(curves_smpl)

c_s_q1, c_s_q3 = np.percentile(curves_smpl, [25, 75], axis = 0)

plt.plot(St_grid, c_s_q1, c = 'aqua', ls = '--', label = 'Q1')

plt.plot(St_grid, c_s_q3, c = 'lime', ls = '--', label = 'Q3')

plt.plot(St_Values, [PowerLaw(v, P50_a, P50_b) for v in St_Values],

         'r-', lw = 3, label = r'$\psi=a\varphi^b$' + ('\na = {0:.3f}'.format\

                        (P50_a)) + ('\nb = {0:.3f}'.format(P50_b)))

plt.plot(Sts, MEs, 'o', c = 'b', alpha = 0.5, label = 'Data')

plt.xlabel(r'$\varphi$', fontsize = 18)

plt.ylabel(r'$\psi$', fontsize = 18)

plt.legend(fontsize = 12)

 

plt.tight_layout()

plt.show()

 

fout.close()

print()

print(foutnm + ' was saved.')

print()

 

 

 

参考文献

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University

岡本安晴(2014)心理学データ分析と測定.勁草書房

岡本安晴(2019)いまさら聞けないPythonでデータ分析.丸善出版

岡本安晴(2025)感覚・知覚測定法.和氣典二・重野純・村上郁也(編)感覚・知覚心理学ハンドブック 第三版(第2章)、誠信書房

Stevens, S. S. (1975). Psychophysics: Introduction to its Perceptual, Neural, and Social Prospects. Wiley.

Wichmann, F. A. & Hill, N. J. (2001). The psychometric function: I. Fitting, sampling, and goodness of fit. Perception & Psychophysics, 63, 1293-1313.

 

 

 

Home