Home

マグニチュード推定法

ベイズ分析(PyStan 3 対応)

サンプルスクリプトは、Ubuntu 22/WSL/Windows 10においてPython 3.12.3/Anaconda仮想環境で実行した(2024.07.04)。

pystan 3.10.0arviz 0.18.0を用いた。Arviz 0.18.0は、次のコマンドでインストールした。

pip install arviz==0.18.0

 

 

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

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

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

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

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

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

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

上記の確率モデルをデータに適用して、ベイズ分析によりパラメータ値を求める。ベイズ分析については、<岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版>の題3部、第9章〜第15章、で説明している。ベイズ分析は、PyStan 3を用いているが、WindowsにおいてPyStan 3を利用する場合は、WindowsUbuntuをインストールしてUbuntu(WSL: Windows Subsysytem for Linux)で使用すればよい。UbuntuWindowsへのインストールなどの手続きは簡単である。説明のウェブサイトを用意しているので、参照されたい。

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

 

 

 

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

 

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

 

 

リスト1 モデル(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> v_a;

    real<lower = 0.0> v_b;

    real<lower = 0.0, upper = 1.0> v_c;

    real<lower = 0.0> v_sgm;

}

transformed parameters {

    real a;

    real b;

    real c;

    real sgm;

    a = 10.0 * v_a + 0.001;

    b = 10.0 * v_b;

    c = v_c * minSt;

    sgm = v_sgm + 0.01;

}

model {

    v_a ~ exponential(0.1);

    v_b ~ exponential(0.1);

    v_c ~ beta(1.0, 1.0);

    v_sgm ~ exponential(0.1);

    for (i in 1:N) {

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

    }

}

 

 

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

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

リスト2のPythonスクリプトをWindowsにインストールされたUbuntu 22上のAnacondaPython 3.12仮想環境で実行した。

次のコマンド

 

(py312) ・・・/me_c_PS3files$ python me_c.py

 

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

 

(py312) ・・・/me_c_PS3files$ python me_c.py

Input Data File (*.xlsx) =

 

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

図1

 

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

入力データファイル名を設定入力すると(次の例では、図1のファイルData.xlsxが設定されている)、データが読み込まれ、次に出力用テキストファイル名の入力が求められる。

 

(py312) ・・・/me_c_PS3files$ python me_c.py

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

Output Text File (*.txt) =

 

適当なテキストファイル名を設定する。

 

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

 

出力ファイル名を設定すると、計算が始まり、MCMCサンプリング後、まずトレース図が表示される(図2)。

図2

 

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

図3

 

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

図4のフォームを閉じると、パラメータの事後分布が表示される(図4)。

図4

 

図4のフォームを閉じると、パラメータの事後分布が表示される(図5)。

図5

 

図5のフォームを閉じると、パラメータの事後分布が表示される(図6)。

図6

 

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

図7

 

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

 

 

を描いたものである。第1四分位数Q!と第3四分位数Q3に対応する値が青色と緑色の破線で描かれている。

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

図7のフォームを閉じると、スクリプトの実行終了となり、端末には以下のようなメッセージが出力されている。

 

Results.txt was saved.

 

(py312) ・・・/me_c_PS3files$

 

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

 

 

Input Data File = Data.xlsx

 

Data...

      1049         5

      1888        10

      1049         7

       180       4.5

     ・

     ・

     ・

       583       3.1

       324       2.9

       180       2.7

 

Summary:

       mean      sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat

a     0.176   0.104   0.043    0.362      0.002    0.001    4988.0    5462.0    1.0

b     0.540   0.075   0.392    0.677      0.001    0.001    5094.0    5407.0    1.0

c    53.070  35.620   0.009  115.363      0.501    0.354    4734.0    5885.0    1.0

sgm   0.329   0.035   0.265    0.394      0.000    0.000    7590.0    7972.0    1.0

 

 

読み込まれたデータと、MCMCサンプリングの統計量が出力されている。

 

パラメータcの事後分布(図6)において、「c=0」の可能性が示されたので、次に「c=0」であるモデル(10.2.1)すなわちモデル(10.2.3)による分析を行う。

 

 

 

モデルの場合

 

確率モデル

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

 

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

 

data {

    int N;

    array[N] real Sts;

    array[N] real MEs;

}

transformed data {

    array[N] real LogMEs;

    LogMEs = log(MEs);

}

parameters {

    real<lower = 0.0> v_a;

    real<lower = 0.0> v_b;

    real<lower = 0.0> v_sgm;

}

transformed parameters {

    real<lower = 0.0> a;

    real<lower = 0.0> b;

    real<lower = 0.0> sgm;

    a = 10.0 * v_a + 0.001;

    b = 10.0 * v_b;

    sgm = 10.0 * v_sgm + 0.001;

}

model {

    v_a ~ exponential(0.1);

    v_b ~ exponential(0.1);

    v_sgm ~ exponential(0.1);

    for (i in 1:N) {

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

    }

}

 

 

リスト3を用いるPythonスクリプトは、リスト4のように用意した。使用法は、リスト2と同じである。スクリプトファイルなどは、まとめて圧縮ファイルme_PS3files.zipとした。ダウンロード展開すれば自由に利用できる。Windowsで実行するときは、Anaconda/Ubuntu/WSL/WindowsPython 3.11仮想環境を利用すればよい。

以下のように、pythonコマンドでスクリプトファイルme.pyを実行すると入力データファイル名を聞いてくる。

 

(py311) ・・・/me_PS3files$ python me.py

Input Data File (*.xlsx) =

 

入力データファイルはリスト2の場合と同様で、図1の形式で用意する。

入力データファイル名を設定すると、次に出力ファイル名の設定が求められる。適当なテキストファイル名を設定する。

 

(py311) ・・・/me_PS3files$ python me.py

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

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

 

出力ファイル名を設定すると、計算が始まる。

MCMCサンプリングが終わると、まずトレースが表示される(図8)。

図8

 

図8のフォームを閉じると、パラメータの散布図が表示される(図9)。

図9

 

図9のフォームを閉じると、パラメータの事後分布が表示される(図10)。

10

 

10のフォームを閉じると、パラメータの事後分布が表示される(図11)。

11

 

11のフォームを閉じると、図12のグラフが表示される。

12

 

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

 

 

を描いたものである。第1四分位数Q!と第3四分位数Q3に対応する値が青色と緑色の破線で描かれている。

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

 

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

端末画面には、下記のようなメッセージが表示される。

 

Results.txt was saved.

 

(py311) ・・・/me_PS3files$

 

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

 

Input Data File = Data.xlsx

 

Data...

1049   5

1888   10

1049   7

  ・

  ・

  ・

583   3.1

324   2.9

180   2.7

 

Summary:

      mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat

a    0.093  0.034   0.038    0.154      0.000    0.000    5237.0    6214.0    1.0

b    0.616  0.055   0.515    0.720      0.001    0.001    5244.0    6029.0    1.0

sgm  0.324  0.035   0.261    0.390      0.000    0.000    7370.0    6850.0    1.0

 

 

 

参考文献

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

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

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.

 

 

 

リスト2 リスト1を利用するPythonスクリプト(me_c.py

 

import numpy as np

import pandas as pd

import stan

import matplotlib as mplt

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

 

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

 

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

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

    print(st, me)

    fout.write(f'{st:>10}{me:>10}\n')

 

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

 

with open('MEc.stan', 'r') as fstan:

    sm = stan.build(fstan.read(), data = Data)

fit = sm.sample(num_samples = 5000)

 

print(fit)

 

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

smry = az.summary(i_data, var_names = ['a', 'b', 'c', 'sgm']) 

print('Summary =\n',smry)

fout.write(f'\nSummary:\n{smry}\n')

 

az.plot_trace(i_data, var_names = ['a', 'b', 'c'])

plt.tight_layout()

plt.show()

 

d_fm = fit.to_frame()

 

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

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

ax.scatter(d_fm['a'], d_fm['b'], d_fm['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()

 

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

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

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

sb.kdeplot(d_fm['a'], clip = (0.0, np.max(d_fm['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.show()

 

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

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

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

sb.kdeplot(d_fm['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.show()

 

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

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

 

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

sb.kdeplot(d_fm['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()

 

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)

 

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

tn_samples = len(d_fm)

n_samples = 500

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

curves_smpl = []

for s in idx_samples:

    va = d_fm['a'][s]

    vb = d_fm['b'][s]

    vc = d_fm['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 = 'b', ls = '--', label = 'Q1')

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

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

         'r-', 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.legend()

plt.show()

   

print()

print(foutnm + ' was saved.')

print()

 

 

 

リスト4 リスト3を用いるPythonスクリプト(me.py

 

import numpy as np

import pandas as pd

import stan

import matplotlib as mplt

import matplotlib.pyplot as plt

import seaborn as sb

import arviz as az

 

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

 

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

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

    print(st, me)

    fout.write(f'{st}   {me}\n')

 

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

 

with open('ME.stan', 'r') as fstan:

    sm = stan.build(fstan.read(), data = Data)

fit = sm.sample(num_samples = 5000)

 

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

smry = az.summary(i_data, var_names = ['a', 'b', 'sgm']) 

print('Summary =\n',smry)

fout.write(f'\nSummary:\n{smry}\n')

 

az.plot_trace(i_data, var_names = ['a', 'b', 'sgm']) 

plt.show()

 

d_fm = fit.to_frame()

 

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

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

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

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

plt.show()

 

mean_a = np.mean(d_fm['a'])

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

 

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

sb.kdeplot(d_fm['a'], clip = (0.0, max(d_fm['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.show()

 

mean_b = np.mean(d_fm['b'])

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

 

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

sb.kdeplot(d_fm['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.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)

 

 

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

tn_samples = len(d_fm)

n_samples = 500

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

curves_smpl = []

for s in idx_samples:

    va = d_fm['a'][s]

    vb = d_fm['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 = 'b', ls = '--', label = 'Q1')

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

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

         'r-', 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.legend()

plt.show()

 

 

fout.close()

print()

print(foutnm + ' was saved.')

print()

 

 

 

Home