Up、en

Psychometric関数（PF）のベイズ推定

２件法、３件法、４件法データの分析

本ウェブサイトのプログラムは、利用者の責任において、自由に御利用下さい。

PyStan 3 対応のウェブサイトを用意した。

２件法より３件法、４件法の方が統計学的に望ましいことを示すシミュレーションが岡本（2017）で報告されている。

なお、psychometric関数の説明を＜岡本安晴「心理学データ分析と測定」＞で行っている。

Psychometric関数（精神測定関数／心理測定関数）のベイズ推定により、閾値（丁度可知差異；Just Noticeable Difference: JND）および主観的等価点（Point of Subjective Equality: PSE）を求めるPythonプログラムを作成してみた。プログラムの使用法などは、本ウェブサイトの後半で説明することとし、まず、背景となる考え方について説明する。

Psychometric関数（ＰＦ）を求めるpsychophysics（精神物理学／心理物理学）実験としては、反応が「大きい／小さい」など（本ウェブサイトでは、物理刺激の比較判断を「大きさ」の次元を例として用いている。他の感覚、例えば「重さ」であれば、「重い／軽い」判断となる）の２カテゴリ強制選択法（２件法）が標準であるが、「わからない／同じ」という判断を認める３件法もある。３件法については、データ分析を単純な平均値に基づいて行う場合には、感覚の弁別力（閾値）の推定値が判断の基準の影響を受けるということが知られているが（岡本、1995）、この問題は確率モデルを設定して分析を行えば回避できる。設定された確率モデルに対して最尤法を用いた場合（岡本、1994、1995）、ベイズ的方法を用いた場合（岡本、2005a）において、感覚の弁別力（閾値）と判断の基準が確率モデルによって分離されて適切に分析できることが報告されている。主観的等価点PSEを求めることを考慮して、２件法と３件法の混合型ともいえる４件法（「大きい」、「どちらかといえば大きい」、「どちらかといえば小さい」、「小さい」の４カテゴリを用いる方法）もある（Okamoto, 2012）。２件法、３件法、４件法、いずれの場合も分析の基本となるのは物理刺激と反応の関係を表すＰＦである。

２件法のモデル

まず、２件法の場合について考える。比較刺激の方が標準刺激より大きいと判断される（と表記）確率をの関数として表したものがＰＦであるが、心理学ではＰＦとして累積正規分布を用いることが多い。累積標準正規分布をで表すと

と書ける。およびは、標準刺激によって決まる値であり、１つの固定したの値について考えるときは、ＰＦはの関数となる。ＰＦをの関数とみてをで表すと、式（１）は

となる。２件法なので

である。

は主観的等価点を表し

である。

　弁別閾（difference threshold; Gescheider, 1999）を表す丁度可知差異（JND; 岡本（2014）、p. 169参照）としては、の倍を用いる。ここで、は、確率0.75に対応する標準正規分布における値である。すなわち、

である。プログラムでは

として計算している。すなわち、

である。

　第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、

とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式

で与えられる。ここで、

とおいた。

　ベイズ的分析では、式（２）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2009）、pp. 153-157」で簡単に説明している。事後分布を岡本（2011）ではグリッド法で求めたが、本ウェブサイトのプログラムでは、３件法、４件法の場合に合わせてＭＣＭＣのMetropolis-within-Gibbsあるいはcomponent-wise version of MH algorithmと呼ばれている方法を用いている。この方法については、「岡本（2014）、pp. 227-232」で説明している。なお、上下法データにベイズ的分析を適用するときには注意が必要である。上下法データは、データ数が少なく偏っている可能性があるが、このため、無情報事前分布を用いると事後分布が拡散することがある。これを避けるためには、弱い事前分布を用いるか、データの収集法に注意する。岡本（2011）では、上下法を２系列用意する二重上下法を用いることにより、適切な事後分布を得ている。

３件法のモデル

　３件法データの分析モデルは、２件法の拡張として考える。３件法では、判断は「標準刺激より比較刺激の方が小さい；」、「とは同じ／どちらが大きいかわからない；」、「よりの方が大きい；」の３つのカテゴリから選ばれる。これらを、「」と「あるいは（と表記）」、および「あるいは（と表記）」と「」の２通りに分けて考える。

「」と「」の組み合わせに対するＰＦを

とおく。は、判断「」と「」を分ける判断の基準値である。

「」と「」の組み合わせに対するＰＦを

とおく。は、判断「」と「」を分ける判断の基準値である。

判断「」の確率は、次式

で与えられる。また、

である。

　JNDは、２件法の場合と同様に

で与えられる。

　PSEについては、２つの判断の基準値との中央に主観的等価点PSEがあると考えられるときは、次式

でPSEが推定されるが、Okamoto (2012)は４件法のデータにおいて式（３）が成り立たない例を報告している。

　いま、第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、

とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式

で与えられる。ここで、

とおいた。

　ベイズ的分析では、式（4）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2009）、pp. 153-157」で簡単に説明している。本ウェブサイトのプログラムでは、事後分布をMetropolis-within-Gibbsあるいはcomponent-wise version of MH algorithmと呼ばれている方法を用いて求めているが、この方法については、「岡本（2014）、pp. 227-232」で説明している。

４件法のモデル

　４件法では、判断は「標準刺激より比較刺激の方が小さい；」、「どちらかといえばよりの方が小さい；」、「どちらかといえばよりの方が大きい；」、「標準刺激より比較刺激の方が大きい；」の４カテゴリから１つを選ぶという形で行われる。これらの判断を以下の３つの２件法データの形式に分けて考える。

（A）「」と「またはまたは（と表記）」

（B）「または（と表記）」と「または（と表記）」

（C）「またはまたは（と表記）」と「」

形式（A）によるデータに対するＰＦを

とおく。形式（B）と（C）に対しても同様に

とおく。

このとき、

である。

　上記モデルにおいて、PSEはによって与えられる。すなわち、

である。ここで、

は仮定されていないことに注意。式（５）が成り立つときは、２つの判断の基準値はPSEから等距離に設定されていることになる。Okamoto (2012)では、４件法データの分析において、式（５）の成り立たない例が報告されている。

　JNDは、２件法および３件法の場合と同様に

で与えられる。

　いま、第試行において提示された比較刺激を、それに対する判断（反応）をで表す。ただし、

とおく。このとき、Ｎ回の試行に対する確率は、独立性の仮定のもとで次式

で与えられる。ここで、

とおいた。

　ベイズ的分析では、式（６）を尤度関数として事後分布が与えられるが、ベイズ的分析については「岡本（2009）、pp. 153-157」で簡単に説明している。本ウェブサイトのプログラムでは、事後分布をMetropolis-within-Gibbsあるいはcomponent-wise version of MH algorithmと呼ばれている方法を用いて求めているが、この方法については、「岡本（2014）、pp. 227-232」で説明している。

比較刺激の範囲と間隔について

恒常法の場合

　比較刺激の範囲と刺激間間隔について少し書いておく。刺激値の範囲は、比較刺激の方が「大きい」および「小さい」の判断が出る範囲であればよい。確率モデルに基づいて分析が行われるので、全系列法のように「大きい」あるいは「小さい」の判断の確率が１になる刺激値を含む必要はない。刺激値の間隔については、岡本（2006a）は、JNDのと3倍とで分析に差がないことを報告しているが、刺激値間隔の選定の参考になると思う。

適応的方法の場合

　試行の始めにおいては、明らかに標準刺激との差が分かる比較刺激値からスタートし、比較刺激値の増減の幅も最初は大きいものから始めて、徐々に小さくする方法がとられるが、このとき、最終的な比較刺激値の増減の幅は過度に小さくする必要はない。これは、確率モデルを用いた分析を行うからで、上の恒常法の場合と同様である。また、先行試行での反応に提示刺激が依存していることを隠すために、２系列の試行を用意してランダムに混合して実施するなどの工夫が必要である。適応的方法としては、Ψ法と呼ばれているエントロピー最小化基準によるものがあり、２件法（岡本、2006b）、３件法（岡本、2004）、４件法（Okamoto, 2012）の場合について説明している。単純な適応的方法では、ターゲットとなる１つのパラメータに焦点を当てた提示刺激の選択が行われるが、Ψ法では、パラメータ全体の事後分布の最適化がエントロピー最適化基準によって行われる。

Fechner問題について

　閾値と感覚尺度の関係を扱う問題はFechner問題と呼ばれているが、これについてはOkamoto (2001, 2005b)で論じている。

引用文献

岡本安晴（1994）「全系列データの最尤法による分析」金沢大学文学部論集 行動科学科篇、pp. 63-71。

岡本安晴（1995）「全系列法データの最尤法による分析について」金沢大学文学部論集 行動科学科篇、pp. 25-30。

Okamoto, Y. (2001). On Fechner’s problem. Studies and Essays, Behavioral Sciences and Philosophy: Kanazawa University, 1-25.

岡本安晴（2004）「３件法に対するBayesian Up-down法」日本行動計量学会第３２回大会発表論文抄録集、pp. 4-5。

岡本安晴（2005a）「恒常法データ（３件法）のベイズ的分析」日本心理学会第６９回大会発表論文集、p. 565。

Okamoto, Y. (2005b). Generalization of Fechner’s approach. Japan Women’s University Journal: Faculty of Integrated Arts and Social Sciences, 107-121.

岡本安晴（2006a）「適応的恒常法の可能性について」日本心理学会第７０回大会発表論文集、p. 559。

岡本安晴（2006b）「計量心理学」培風館

岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館

岡本安晴（2011）「上下法のベイズ的分析」基礎心理学研究、44−55。

Okamoto, Y. (2012). An experimental analysis of psychometric functions in a threshold discrimination task with four response categories. Japanese Psychological Research, 368-377.

岡本安晴（2014）「心理学データ分析と測定」勁草書房

岡本安晴（2017）「恒常法における反応カテゴリ数の影響」基礎心理学研究、85−91．

プログラム

２件法データのプログラム、３件法データのプログラム、４件法データのプログラムと順番に説明する。

２件法データの分析プログラム

　２件法データの分析プログラムをサンプル入力データファイルとともに圧縮ファイルPM2Cat.zipにまとめた。ファイル名PM2Cat.zipをクリックしてダウンロード・解凍すれば利用できる。プログラムを実行すると、入力データファイル名を聞いてくるので、設定する（図１）。

図１

入力データファイルは、図２の形式で用意する。

・

・

・

図２

先頭がスラッシュ / で始まる行の次の行からデータを

データID　　　比較刺激値　　　　反応（判断）

の順序で１組ずつ書いていく。データIDは各行のデータを区別するものであるが、文字列でよく、また同じものがあってもよい。反応は、「比較刺激値の方が小さい」というときは１、「比較刺激の方が大きい」というときは２を設定する。

最後のデータの次の行は、スラッシュ / で始まる行を置いて、データの終わりとする。図２のデータファイルはプログラムと一緒に圧縮ファイルPM2Cat.zipに含まれている。

入力データファイル名を設定してEnterキーを押すと、次に出力ファイル名を聞いてくる（図１）。適当なテキストファイル名を設定する。出力ファイル名の設定後、計算が始まる。

図３

図３に示されているように、「The main MCMC ended.」という出力が表示されたら、MCMCによる事後分布からのサンプリングが終了であり、このサンプルに基づいて分析結果が表示される。まず、図４のようなグラフが表示される。このウィンドウは、他のウィンドウの後ろに表示されると気が付かないので注意が必要である。

図４

青の曲線がＰＦである。このＰＦはの事後分布の中央値をPSEの点推定値とし、の事後分布の中央値をJNDの点推定値として描かれている。緑色の小円はデータである。ただし、小円が重ならないように、ランダムな小さい値を加えて表示してある。すなわち

としたとき、

はランダムな小さい正数である。

図４において、フロッピーディスクのアイコンをクリックすると画像を保存することができる。右上のX印のクリックでウィンドウが閉じられ、次のウィンドウが表示される（図５）。

図５

図５は、の事後分布を示すものである。図５のウィンドウを閉じると、図６のウィンドウが表示される。

図６

図６には、JND（）の事後分布が示されている。図６のウィンドウを閉じると、出力ファイルが保存されたというメッセージが出力されて（図７）、プログラムの終了となる。

図７

出力ファイルをエディタで開くと、図８のような内容になっている。

・

・

・

図８

先頭に入力データファイル名が書き出され、続いて読み込んだデータが書き出されている。最後に、PSEの事後分布の中央値、四分位数Q1とQ3,、そして９５％CIが出力されている。JNDの事後分布も同様の統計量が出力されている。

プログラムのメインソースファイルは以下のようである。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ mainprg.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

from multiprocessing.pool import Pool

import math

from mcmcFunc import *

from RNGen import *

from subprg import *

#

# Yasuharu Okamoto, 2016.11

#

if __name__ == '__main__':

X, n, meanSt, sgmSt, fout, nm_fout = InputData()

NSimu = 1000

params = {'jump' : 0,

'n' : n,

'X' : X,

'NSimu' : NSimu,

'initMu' : meanSt,

'initSgm' : sgmSt,

'genSgmMu' : sgmSt/10.0,

'genSgmSgm' : sgmSt/10.0}

#

# Adjust the parameters in MCMC

#

istep = 0

while True:

print("\nStep-{} started.".format(istep))

Mu, Sgm, cnt_mu, cnt_sgm = mcmcFunc( params )

sum_mu = 0.0

sum_sgm = 0.0

for t in range(1, NSimu + 1):

sum_mu += Mu[t]

sum_sgm += Sgm[t]

params['initMu'] = sum_mu / NSimu

params['initSgm'] = sum_sgm / NSimu

print("meanMu = {}".format(params['initMu']))

print("meanSgm = {}".format(params['initSgm']))

ck_cnt = 0

r_mu = cnt_mu / NSimu

print("r_mu = {}".format(r_mu))

v, ck_cnt = scale_sgm( r_mu, ck_cnt )

params['genSgmMu'] = params['genSgmMu'] * v

print("genSgmMu = {}".format(params['genSgmMu']))

r_sgm = cnt_sgm / NSimu;

print("r_sgm = {}".format(r_sgm))

v, ck_cnt = scale_sgm( r_sgm, ck_cnt )

params['genSgmSgm'] = params['genSgmSgm'] * v

print("genSgmSgm = {}".format(params['genSgmSgm']))

print("ck_cnt = {}".format(ck_cnt))

if ck_cnt == 0:

break

istep += 1

#

# The main MCMCs

#

NSimu = 3000

params['NSimu'] = NSimu

n_mcmc = 4

ary_params = []

for i in range(0, n_mcmc):

temp = params.copy()

temp['jump'] = i + 1

ary_params.append(temp)

print("\nThe main MCMCs started.")

pool = Pool()

results = pool.map(mcmcFunc, ary_params) # Multiprocessing

print("\nThe main MCMCs ended.\n")

mrg_mu = []

mrg_sgm = []

for temp_Mu, temp_Sgm, cnt_mu, cnt_sgm in results:

mrg_mu += temp_Mu[1:NSimu + 1]

mrg_sgm += temp_Sgm[1:NSimu + 1]

print("Length(mrg_mu) = {}".format(len(mrg_mu)))

print("Length(mrg_sgm) = {}".format(len(mrg_sgm)))

mrg_mu.sort()

mrg_sgm.sort()

DispResults( NSimu, n_mcmc, X, n, mrg_mu, mrg_sgm, fout)

fout.close()

print("\nOutput file {} was saved.".format(nm_fout))

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

メインのMCMCにおいてマルチプロセッシングが利用されているが、子プロセスのprintによる出力は、端末（コンソールウィンドウ）でプログラムを実行すると、見ることができる（図９）。

図９

３件法データの分析プログラム

　２件法データの分析プログラムと同様であるので、簡単に説明する。プログラムファイルとサンプル入力データファイルは圧縮ファイルPM3Cat.zipにまとめた。ファイル名PM3Cat.zipをクリックしてダウンロード・解凍すれば利用できる。プログラムを実行すると入データファイル名の入力が求められる（図１０）。

図１０

入力データファイルは図１１に示す形式で用意する。

・

・

・

図１１

先頭がスラッシュ / で始まる行の次の行から１行に1組ずつのデータを

　　　データID　　　　提示刺激　　　　判断（反応）

の形式で書いていく。データIDは文字列として扱われ1つ1つのデータに付けるが、同じものがあってもかまわない。反応（判断）は、比較刺激の方が小さいなら１、同じなら２、大きいなら３を設定する。データの最後は、スラッシュ / で始まる行を置く。

入力データファイル名を設定すると、次に出力用テキストファイル名の設定が求められる。適当な出力テキストファイル名を設定すると計算が始まる。メインのMCMCが終了すると、まず２本のPM関数を描いたウィンドウが表示される（図１２）。

図１２

データ点は

とランダムな小さな値を加えて、データの小円の重なりを避けている。

次に、２つの基準値C1とC2の事後分布のグラフが表示される。

図１３

最後に、JNDの事後分布が表示される（図１４）。

図１４

プログラムの実行終了後、出力ファイルを開くと図１５のようになっている。

・

・

・

図１５

基準値C1とC2、およびJNDの事後分布の中央値、四分位数Q1とQ3、９５％CI区間が出力されている。

メインプログラムのソースコードを以下に示す。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝ mainprg3Cat.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

from multiprocessing.pool import Pool

import matplotlib.pyplot as plt

import numpy as np

import math

from mcmcFunc3Cat import *

from RNGen import *

from subprg import *

#

# Yasuharu Okamoto, 2016.11

#

if __name__ == '__main__':

X, n, meanSt, sgmSt, fout, nm_fout = InputData()

NSimu = 1000

params = {'jump' : 0,

'n' : n,

'X' : X,

'NSimu' : NSimu,

'initC1' : meanSt - sgmSt,

'initC2' : meanSt + sgmSt,

'initSgm' : sgmSt,

'genSgmC1' : sgmSt/10.0,

'genSgmC2' : sgmSt/10.0,

'genSgmSgm' : sgmSt/10.0}

#

# Adjust the parameters in MCMC

#

istep = 0

while True:

print("Step-{} started.".format(istep))

C1, C2, Sgm, cnt_C1, cnt_C2, cnt_sgm = mcmcFunc( params )

sum_C1 = 0.0

sum_C2 = 0.0

sum_sgm = 0.0

for t in range(1, NSimu + 1):

sum_C1 += C1[t]

sum_C2 += C2[t]

sum_sgm += Sgm[t]

params['initC1'] = sum_C1 / NSimu

params['initC2'] = sum_C2 / NSimu

params['initSgm'] = sum_sgm / NSimu

print("meanC1 = {}".format(params['initC1']))

print('meanC2 = {}'.format(params['initC2']))

print("meanSgm = {}".format(params['initSgm']))

ck_cnt = 0

r_C1 = cnt_C1 / NSimu

print("r_C1 = {}".format(r_C1))

v, ck_cnt = scale_sgm( r_C1, ck_cnt )

params['genSgmC1'] = params['genSgmC1'] * v

print("genSgmC1 = {}".format(params['genSgmC1']))

r_C2 = cnt_C2 / NSimu

print("r_C2 = {}".format(r_C2))

v, ck_cnt = scale_sgm( r_C2, ck_cnt )

params['genSgmC2'] = params['genSgmC2'] * v

print("genSgmC2 = {}".format(params['genSgmC2']))

r_sgm = cnt_sgm / NSimu;

print("r_sgm = {}".format(r_sgm))

v, ck_cnt = scale_sgm( r_sgm, ck_cnt )

params['genSgmSgm'] = params['genSgmSgm'] * v

print("genSgmSgm = {}".format(params['genSgmSgm']))

print("ck_cnt = {}".format(ck_cnt))

if ck_cnt == 0:

break

istep += 1

#

# The main MCMCs

#

NSimu = 3000

params['NSimu'] = NSimu

n_mcmc = 4

ary_params = []

for i in range(0, n_mcmc):

temp = params.copy()

temp['jump'] = i + 1

ary_params.append(temp)

print("\nThe main MCMCs started.")

pool = Pool()

results = pool.map(mcmcFunc, ary_params)

print("\nThe main MCMCs ended.")

mrg_C1 = []

mrg_C2 = []

mrg_sgm = []

for temp_C1, temp_C2, temp_Sgm, cnt_C1, cnt_C2, cnt_sgm in results:

mrg_C1 += temp_C1[1:NSimu + 1]

mrg_C2 += temp_C2[1:NSimu + 1]

mrg_sgm += temp_Sgm[1:NSimu + 1]

print("Length(mrg_C1) = {}".format(len(mrg_C1)))

print('Length(mrg_C2) = {}'.format(len(mrg_C2)))

print("Length(mrg_sgm) = {}".format(len(mrg_sgm)))

mrg_C1.sort()

mrg_C2.sort()

mrg_sgm.sort()

DispResults( NSimu, n_mcmc, X, n, mrg_C1, mrg_C2, mrg_sgm, fout )

fout.close()

print("\nOutput file {} was saved.".format(nm_fout))

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

４件法データの分析プログラム

　２件法の場合と同様であるので、簡単に説明する。プログラムとサンプル入力データファイルは、圧縮ファイルPM4Cat.zipにまとめた。ファイル名PM4Cat.zipをクリックしてダウンロード・解凍すれば利用できる。プログラムを実行すると、入力データファイル名の設定が求められる（図１６）。

図１６

入力データファイルは、図１７の形式で用意する。

・

・

・

図１７

先頭がスラッシュ / で始まる行の次の行から１行に1組ずつのデータを

　　　データID　　　　提示刺激　　　　判断（反応）

の形式で書いていく。データIDは1つ1つのデータに付けるが、文字列として扱われ、同じものがあってもかまわない。反応（判断）は、比較刺激の方が小さいなら１、たぶん小さいなら２、たぶん大きいなら３、大きいなら４を設定する。データの最後は、スラッシュ / で始まる行を置く。

入力データファイル名に続いて出力用テキストファイル名を設定すると計算が始まる。メインのMCMCが終了すると、PFとデータ点を描いたウィンドウが表示される（図１８）。

図１８

３本のＰＭ関数が描かれている。データ点は

とランダムな小さい値を加えて、データの小円の重なりを避けている。

次に、２つの基準値C1とC2,およびＰＳＥの事後分布のグラフが表示される。

図１９

最後に、JNDの事後分布が表示される。

図２０

プログラムの終了後、出力ファイルを開くと、図２１のように各事後分布の中央値、四分位数Q1とQ3、および９５％CIが出力されている。

図２１

プログラムのメインソースコードファイルを以下に示す。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝ mainprg4Cat.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

from multiprocessing.pool import Pool

from mcmcFunc4Cat import *

from subprgs4Cat import *

#

# Yasuharu Okamoto, 2016.11

#

if __name__ == '__main__':

X, n, meanSt, sgmSt, fout, nm_fout = InputData()

NSimu = 1000

params = {'jump' : 0,

'n' : n,

'X' : X,

'NSimu' : NSimu,

'initC1' : meanSt - sgmSt,

'initMu' : meanSt,

'initC2' : meanSt + sgmSt,

'initSgm' : sgmSt,

'genSgmC1' : sgmSt/10.0,

'genSgmMu' : sgmSt / 10.0,

'genSgmC2' : sgmSt/10.0,

'genSgmSgm' : sgmSt/10.0}

istep = 0

while True:

print("Step-{} started.".format(istep))

C1, Mu, C2, Sgm, cnt_C1, cnt_Mu, cnt_C2, cnt_sgm = mcmcFunc( params )

sum_C1 = 0.0

sum_Mu = 0.0

sum_C2 = 0.0

sum_sgm = 0.0

for t in range(1, NSimu + 1):

sum_C1 += C1[t]

sum_Mu += Mu[t]

sum_C2 += C2[t]

sum_sgm += Sgm[t]

params['initC1'] = sum_C1 / NSimu

params['initMu'] = sum_Mu / NSimu

params['initC2'] = sum_C2 / NSimu

params['initSgm'] = sum_sgm / NSimu

print("meanC1 = {}".format(params['initC1']))

print("meanMu = {}".format(params['initMu']))

print('meanC2 = {}'.format(params['initC2']))

print("meanSgm = {}".format(params['initSgm']))

ck_cnt = 0

r_C1 = cnt_C1 / NSimu

print("r_C1 = {}".format(r_C1))

v, ck_cnt = scale_sgm( r_C1, ck_cnt )

params['genSgmC1'] = params['genSgmC1'] * v

print("genSgmC1 = {}".format(params['genSgmC1']))

r_Mu = cnt_Mu / NSimu

print('r_Mu = {}'.format(r_Mu))

v, ck_cnt = scale_sgm( r_Mu, ck_cnt )

params['genSgmMu'] = params['genSgmMu'] * v

print("genSgmMu = {}".format(params['genSgmMu']))

r_C2 = cnt_C2 / NSimu

print("r_C2 = {}".format(r_C2))

v, ck_cnt = scale_sgm( r_C2, ck_cnt )

params['genSgmC2'] = params['genSgmC2'] * v

print("genSgmC2 = {}".format(params['genSgmC2']))

r_sgm = cnt_sgm / NSimu;

print("r_sgm = {}".format(r_sgm))

v, ck_cnt = scale_sgm( r_sgm, ck_cnt )

params['genSgmSgm'] = params['genSgmSgm'] * v

print("genSgmSgm = {}".format(params['genSgmSgm']))

print("ck_cnt = {}".format(ck_cnt))

if ck_cnt == 0:

break

istep += 1

NSimu = 3000

params['NSimu'] = NSimu

n_mcmc = 4

ary_params = []

for i in range(0, n_mcmc):

temp = params.copy()

temp['jump'] = i + 1

ary_params.append(temp)

print("\nThe main MCMCs started.")

pool = Pool()

results = pool.map(mcmcFunc, ary_params)

print("\nThe main MCMCs ended.")

mrg_C1 = []

mrg_Mu = []

mrg_C2 = []

mrg_sgm = []

for temp_C1, temp_Mu, temp_C2, temp_Sgm, cnt_C1, cnt_Mu, cnt_C2, cnt_sgm in results:

mrg_C1 += temp_C1[1:NSimu + 1]

mrg_Mu += temp_Mu[1:NSimu + 1]

mrg_C2 += temp_C2[1:NSimu + 1]

mrg_sgm += temp_Sgm[1:NSimu + 1]

print("Length(mrg_C1) = {}".format(len(mrg_C1)))

print("Length(mrg_Mu) = {}".format(len(mrg_Mu)))

print('Length(mrg_C2) = {}'.format(len(mrg_C2)))

print("Length(mrg_sgm) = {}".format(len(mrg_sgm)))

mrg_C1.sort()

mrg_Mu.sort()

mrg_C2.sort()

mrg_sgm.sort()

DispResults( NSimu, n_mcmc, X, n, mrg_C1, mrg_Mu, mrg_C2, mrg_sgm, fout )

fout.close()

print("\nOutput file {} was saved.".format(nm_fout))

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Up