Up

3点識別試験ベイズ分析

d’の回帰分析例

岡本安晴、2021.122022.01

 

3つの刺激の内、2つが同じ刺激Aで1つが異なる刺激Bであるとき、3つの中で他の2つと異なるものとして正しくBが選ばれる(あるいは、Bが他の2つと最も異なるものとして選ばれる)確率を調べる実験(3点識別試験:triangular test, 3AFC Oddity)において、実験協力者の特性(年齢など)によりABの感覚の差がどのように変化しているかを調べる分析例のプログラムを作成してみた。

信号検出理論に基づいて、刺激Aおよび刺激Bの感覚が、それぞれ独立に平均および平均の分散の正規分布(ガウス分布)に従うものとする(サーストンモデル)。2つの刺激Aの感覚をおよびで表し、刺激Bの感覚をで表すと

である。

3点識別試験法において正しくBが選ばれるのは

が成り立つときである(Frijters et al., 1980, p. 177)。

いま、信号検出理論の慣例に倣って、感覚の平均値の差を標準偏差を単位として表し、

とおく。このとき、正答確率(Bが選ばれる確率)と感覚の差の関係は次式で表される(Frijters et al., 1980, eq.(1))。

式(1)において、の正負を入れ換え(で置き換え)ても値は変わらない。したがって、感覚量の座標の向きをBの平均の方がAの平均より大きくなる方向にとってとおいても一般性は失われない。式(1)より、感覚の差に対する正答確率が与えられる。

3点識別試験データが、実験参加者の判断の正誤と参加者の特性のペアで与えられているとする。

例えば、図1のデータ(仮想データ)では、判断の正誤が変数Rの値として1(正答)・0(誤答)で与えられ、判断者の年齢が変数Aの値として記入されている。変数CaseIDはデータの識別値で図1では通し番号が記入されている。CaseIDが1の判断者の反応は正答(1)であり、年齢は33歳である。

図1

 

感覚の差d’と年齢の関係を調べるが、d’には値が負でないという制約があるので、リンク関数(link function; Gelman et al., 2021)を用いてd’の範囲をからに変換して年齢(変数Xで表す)の1次式として表す(単回帰モデル)。すなわち、

とおく。このとき、

である。

式(1)と式(2)より、データに対する尤度はパラメータabの関数として与えられる。

である。上式に基づき、リスト1のMCMCサンプリングにおける事後確率の計算は以下のようになっている。Rsが判断Rを表し、Agsが独立変数Xである。

 

v = my_log(phi(a, 0.0, 10.0)) + my_log(phi(b, 0.0, 10.0))

        for i in range(len(self.Rs)):

            dprime = np.exp(a + b * self.Ags[i])

            if dprime > 10.0:

                dprime = 10.0

            p = self.TblPDprime[f'{dprime:.2f}']

            v += self.Rs[i] * my_log(p) +\

                 (1 - self.Rs[i]) * my_log(1- p)

 

独立変数Agsは標準化された値(プログラム内で自動的に標準化される)であるという前提の下で、事前分布は弱情報事前分布として平均0,標準偏差10の正規分布が設定されている。また、式(1)の計算は時間が掛かるので、予めまとめて計算しておき辞書TblPDprimeに格納している。

リスト1のプログラムでは、MCMCサンプリングはMetropolis-within-Gibbs (Robert & Casella, 2010)で行っている(岡本、2014pp.229-232)。

MCMCサンプリングは、以下のスクリプト(リスト内)によって4本のチェインが同時に並列実行されている。

 

    #

    #         Parameter values for multiproccing, each of which runs one chain in MCMC

    #

    params_t = {'Rs': deepcopy(Rs), 'Ags': deepcopy(Ags), 'TblPDprime': deepcopy(TblPDprime),

                'init_a': init_a, 'init_b': init_b, 'sgm_a': sgm_a, 'sgm_b': sgm_b,

                   'n_samples': n_samples, 'chain_id': None}

    #

    #    List of parameter values for multiprocessing

    #

    params = []

    for i in range(4):

        p_t = deepcopy(params_t)

        p_t['chain_id'] = i

        params.append(p_t)

 

    print('MCMC sampling by multiprocessing started.') 

    with Pool() as pool:

        Results = pool.map(mm.MCMCs, params)    #  Run the multiprocessing

 

4本のチェインが独立に実行されるために必要な独立な乱数生成オブジェクトは、次のスクリプト(リスト内)により用意している。

 

        rn = RNp2to191()

        #

        #  Get independent chains

        #  cf. http://y-okamoto-psy1949.la.coocan.jp/Python/sampleprgs/rngeneration/en/

        #

        for i in range(self.chain_id):

            rn.jump(100)

 

この乱数生成オブジェクトの詳しい説明は、このウェブサイトを参照されたい。

 

リスト1のモジュール(ファイル名mymodule.py)を用いてベイズ分析を行うプログラム(main_tritestdrgrssn.py)はリスト2に示した。

プログラムファイルなどは圧縮ファイルTriTestDRgrssnFiles.zipにまとめた。ダウンロード展開すれば、自由に利用できる。

リスト2のスクリプトを実行すると、入力データファイル名の入力が求められる。

 

(py39) PS D:\xxxxx\TriTestDRgrssnFiles> python .\main_tritestdrgrssn.py

Input file (*.csv) = Data.csv

 

入力データは、上の図1の形式で用意する。1行目に変数名を書く。1列目はデータの識別変数を置く。変数値は、数値でも文字列でもよい。2列目は、判断の正誤を記入する。正答なら1,誤答なら0である。3列目には、回帰分析の独立変数を置く。

作成したデータは、ファイル拡張子として「.csv」を選ぶとCSV形式で保存される。入力データファイルは、CSV形式のファイルである。

入力ファイル名(上の例の場合、Data.csv)を設定してEnterキーを押すと、計算が始まる。

計算に必要な辞書の作成が終わると、MCMCサンプリングが始まるが、まずproposal分布の調整が行われ、その後、メインのサンプリングが4本のマルコフ鎖で並行して行われる。このサンプリングが終わると、まずパラメータaの4本のマルコフ鎖の分布とトレースのグラフが表示される(図2)。

図2

 

トレースのグラフの上にはRhat)の値(Gelman et al., 2014, pp.284-285)が表示されている。

図2のグラフのウィンドウを閉じる(ウィンドウの右上角のX印をクリックする)と、次にパラメータbの分布とトレースが表示される(図3)。

図3

 

図3のグラフのウィンドウを閉じると、4本のサンプリングをまとめたパラメータaの分布とモードが表示される(図4)。

図4

 

図4のウィンドウを閉じると、パラメータbの4本のマルコフ鎖をまとめた分布とモードが表示される(図5)。

図5

 

図5のグラフのウィンドウを閉じると、正答確率の回帰曲線とデータ値のグラフが表示される(図6)。

図6

 

正答確率の回帰曲線は、パラメータの事後分布のモードをパラメータ値として描かれている。データ値は、独立変数を4つの区間に分割したときの各区間内の平均正答率である。データ値を回帰曲線がよく表していると言える。

図6のウィンドウを閉じると、事後分布のモードをパラメータ値とするd’の回帰曲線のグラフが表示される。

図7

 

図7のグラフのウィンドウを閉じると、プログラムの実行終了である。

 

 regrslts.txt was saved.

 

(py39) PS D:\xxxxx\TriTestDRgrssnFiles>

 

スクリプトファイルのフォルダに出力されたファイルregrslts.txtを開くと、回帰モデルの式が出力されている。

 

Regression model of Pc on the standardized independent variable Z:

P(R = 1) = Pc(exp(0.70 + -0.36 * Z))

 

Regression model of Pc on the raw independent variable X:

P(R = 1) = Pc(exp(1.96 + -0.03 * X))

 

Regression model of d' on the standardized independent variable Z:

d' = exp(0.70 + -0.36 * Z)

 

Regression model of d' on the raw independent variable X:

d' = exp(1.96 + -0.03 * X)

 

グラフも画像ファイルとして保存されている。

 

 

参考文献

Frijters, J. E. R., Kooistra, A., & Vereijken, P. F. G. (1980). Tables od d' for the triangular method and the 3-AFC signal detection procedure. Perception & Psychophysics, 1980, 27, 176-178.

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

Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., & Rubin, D.B. (2014). Bayesian data analysis, third edition. CRC Press.

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

Robert, C. P., & Casella, G. (2010). Monte Carlo statistical methods. Springer

 

 

リスト1 MCMCサンプリングなどのモジュール(ファイル名mymodule.py

"""

          Yasuharu Okamoto, 2021.12

"""

import numpy as np

import scipy.stats as ss

import scipy.integrate as si

from RNGen import *           #   Random number generator

 

CumPhi = ss.norm.cdf

phi = ss.norm.pdf

 

def f_phi(z):

    return phi(z)

 

 

class k_Pc:

    def __init__(self, d):

        self.d = d

 

    def f(self,u):

        """

               Frijters, Kooistra & Vereijken.(1980).

               Perception & Psychophysics, 1980, 27, 176-178.

        """

        v = (CumPhi(-u*(3**0.5) + self.d * ((2/3)**0.5))

             + CumPhi(-u*(3**0.5) - self.d * ((2/3)**0.5))) \

            * phi(u)

        return v

 

def f_Pc(d):

    int_k_Pc = k_Pc(d)

    v = si.quad(int_k_Pc.f, 0.0, +np.inf)

    return v[0] * 2.0

 

 

def calcPDprimes():

    """

            Dictionary:  key = d'

                         value = probability of the correct response

    """

    pdprimes = {}

    dprimes = np.linspace(0.0, 10.0, 1001)

    for i, d in enumerate(dprimes):

        if i % 100 == 0:

            print(i, end = '\r')

        p = f_Pc(d)

        pdprimes[f'{d:.2f}'] = p

    return pdprimes

 

 

def scale_sgm(acpt_r, ck):

    """

            Check the size of acceptance proportion

    """

    v = 1.0

    if (acpt_r > 0.5):

        ck += 1

        v = 1.0 + 9.0 * (acpt_r - 0.5) / 0.5

    elif (acpt_r < 0.3):

        ck += 1

        v = 1.0 / (1.0 + 9.0 * (0.3 - acpt_r) / 0.3)

    return v, ck

 

def my_log(v):

    if v <= 0.0:

        return -750.0

    else:

        return math.log(v)

 

class MCMC:

    """

                  Metropolis-within-Gibbs algorithm

                  cf. Robert,C.P. & Casella, G. (2010)

                  Monte Carlo Statistical Methods. p.393

                  岡本安晴 (2014).心理学データ分析と測定.pp.229-232.

    """

    def __init__(self, Rs, Ags, TblPDprime, init_a, init_b,

                 sgm_a, sgm_b, n_samples, chain_id):

        self.Rs = Rs

        self.Ags = Ags

        self.TblPDprime = TblPDprime

        self.n_samples = n_samples

        self.chain_id = chain_id

        self.a_coeff = np.empty(self.n_samples + 1)

        self.b_coeff = np.empty(self.n_samples + 1)

        self.a_coeff[0] = init_a

        self.b_coeff[0] = init_b

        self.sgm_a = sgm_a

        self.sgm_b = sgm_b

        self.naccpt_a = 0

        self.naccpt_b = 0

 

    def LL(self, a, b):

        """

                Log likelihood of the parameters a and b

        """

        v = my_log(phi(a, 0.0, 10.0)) + my_log(phi(b, 0.0, 10.0))

        for i in range(len(self.Rs)):

            dprime = np.exp(a + b * self.Ags[i])

            if dprime > 10.0:

                dprime = 10.0

            p = self.TblPDprime[f'{dprime:.2f}']

            v += self.Rs[i] * my_log(p) +\

                 (1 - self.Rs[i]) * my_log(1- p)

              

        return v

   

 

    def do_mcmc(self):

        """

                Doing MCMC sampling

        """

        import tkinter as tk    #   For message window

       

        if self.chain_id < 0:

            msg0 = 'Adjusting sgm...step-{}'.format(-self.chain_id)

        else:

            msg0 = 'Proc-{}'.format(self.chain_id)

       

        root = tk.Tk()

        root.title(msg0)

        cnvs = tk.Canvas(width = 700, height = 300, background = '#00ffff')

        cnvs.pack()

        cnvs.create_text(350, 150, text = '', font = ('', 90), fill = 'blue',

                         tags = 'MyText')

        

        rn = RNp2to191()

        #

        #  Get independent chains

        #  cf. http://y-okamoto-psy1949.la.coocan.jp/Python/sampleprgs/rngeneration/en/

        #

        for i in range(self.chain_id):

            rn.jump(100)

       

        prev_LL = self.LL(self.a_coeff[0], self.b_coeff[0])   #   The initial values

        den = prev_LL                                         #   The initial value

        for t in range(self.n_samples):

            if (t % 100) == 0:

                cnvs.itemconfig('MyText', text = f'{t}/{self.n_samples}')

                cnvs.update()

            #

            #   The next value for the parameter a

            #

            y = rn.normalMS(self.a_coeff[t], self.sgm_a)

            num = self.LL(y, self.b_coeff[t])

            den = prev_LL

            accpt = math.exp(num - den)

            if rn.uni() < accpt:

                self.a_coeff[t + 1] = y

                self.naccpt_a += 1

                prev_LL = num

            else:

                self.a_coeff[t + 1] = self.a_coeff[t]

                prev_LL = den

            #

            #   The next value for the parameter b

            #

            y = rn.normalMS(self.b_coeff[t], self.sgm_b)

            num = self.LL(self.a_coeff[t+1], y)

            den = prev_LL

            accpt = math.exp(num - den)

            if rn.uni() < accpt:

                self.b_coeff[t + 1] = y

                self.naccpt_b += 1

                prev_LL = num

            else:

                self.b_coeff[t + 1] = self.b_coeff[t]

                prev_LL = den

 

 

        root.destroy()

 

        return self.naccpt_a, self.naccpt_b, self.n_samples, \

               self.a_coeff[1:], self.b_coeff[1:]

 

 

def MCMCs(params):

    """

            The function for muliprocessing

    """

    my_mcmc = MCMC(params['Rs'], params['Ags'], params['TblPDprime'],

                   params['init_a'],

                   params['init_b'], params['sgm_a'], params['sgm_b'],

                   params['n_samples'], params['chain_id'])

    r = my_mcmc.do_mcmc()

    return r

 

 

 

def calcRhat(psi):

    """

             Gelman,A., Carlin,J.B.,Stern,H.S., Dunson,D.B.,

             Vehtari,A, & Rubin,D.B. (2014)

             Bayesian Data Analysis, third edition. pp.284-285.

    """

    psi = np.array(psi)

    psidj = psi.mean(axis = 1)

    psidd = psi.mean()

    n = len(psi[0])

    m = len(psi)

    B = (((psidj - psidd) ** 2).sum()) * n / (m - 1)

    sj2 = []

    for j in range(m):

        sj2.append(((psi[j] - psidj[j]) ** 2).sum() / (n - 1))

    sj2 = np.array(sj2)

    W = (sj2.sum()) / m

    var_hat = W * (n - 1) / n + B / n

    Rhat = (var_hat / W) ** 0.5

    return Rhat       

 

   

def EstMode(samples, a = 0.05, n_points = 10000):

    """

        Calculatte a MAP estimate from a KDE graph on [Lp, Up]

        Lp and Up are 100*a/2 and 100(1-a/2) percentile points of samples

    """

    import numpy as np

    import scipy.stats as ss

    Lp, Up = np.percentile(samples, [100 * a/2, 100 * (1 - a/2)])  # import numpy as np

    coord = np.linspace(Lp, Up, n_points) 

    est_pdf = ss.gaussian_kde(samples).pdf(coord)  # import scipy.stats as ss

    mode_idx = np.argmax(est_pdf) 

    Mode_Est = coord[mode_idx]       

    return Mode_Est

 

 

リスト2 リスト1のモジュールを用いて3点識別試験データを分析するプログラム(ファイル名main_tritestdrgrssn.py

"""

         Yasuharu Okamoto, 2021.12, 2022.01

"""

import csv

from copy import deepcopy

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

from multiprocessing.pool import Pool

import mymodule as mm     #  a personal module, file name: mymodule.py

 

if __name__ == '__main__':

    fin_nm = input('Input file (*.csv) = ')

    with open(fin_nm, 'r') as f:

        Data_in = [v for v in csv.reader(f)]

       

    Data_in = np.array(Data_in)

    Rs = [int(v) for v in Data_in.T[1][1:]]

    Ags = [float(v) for v in Data_in.T[2][1:]]

 

    print(len(Rs), len(Ags))

    Rs = np.array(Rs)

    Ags = np.array(Ags)

    Ags_raw = deepcopy(Ags)

    Ags_mean = Ags.mean()

    Ags_std = Ags.std()

    Ags = (Ags - Ags_mean) / Ags_std

    #   Construction of a dictionary for d's and probabilities of correct response

    print("Preparing the dictionary of d's and probabilties")

    TblPDprime = mm.calcPDprimes()    

 

    #

    #      Adjusting the proposal distributions

    #

    i_step = -1

    sgm_a = 0.5

    sgm_b = 0.5

    init_a = 0.0

    init_b = 0.0

    print('Adjusting the proposals')

    while True:

        print('step-{0:} started.'.format(-i_step))

        my_mcmc = mm.MCMC(deepcopy(Rs), deepcopy(Ags), deepcopy(TblPDprime), init_a, init_b, sgm_a, sgm_b, 2000, i_step)

        n_accpt_a, n_accpt_b, L_chain, chain_a, chain_b = my_mcmc.do_mcmc()

        acpt_r_a = n_accpt_a / L_chain

        acpt_r_b = n_accpt_b / L_chain

                                

        ck_a = 0

        v, ck_a = mm.scale_sgm(acpt_r_a, ck_a)  #   Check the acceptance proportions

        init_a = chain_a.mean()

        if ck_a > 0:

            sgm_a *= v    #   Adjust the sigma of the normal distribution for proposal

       

        ck_b = 0

        v, ck_b = mm.scale_sgm(acpt_r_b, ck_b) 

        init_b = chain_b.mean()

        if ck_b > 0:

            sgm_b *= v

           

        if ck_a + ck_b == 0:   #  The sizes of the sigmas OK?

            break              #  Adjustment ends,

        i_step -= 1

   

    n_samples = 5000      #    Sampling size for each chain in MCMC sampling

           

    #

    #         Parameter values for multiproccing, each of which runs one chain in MCMC

    #

    params_t = {'Rs': deepcopy(Rs), 'Ags': deepcopy(Ags), 'TblPDprime': deepcopy(TblPDprime),

                'init_a': init_a, 'init_b': init_b, 'sgm_a': sgm_a, 'sgm_b': sgm_b,

                   'n_samples': n_samples, 'chain_id': None}

    #

    #    List of parameter values for multiprocessing

    #

    params = []

    for i in range(4):

        p_t = deepcopy(params_t)

        p_t['chain_id'] = i

        params.append(p_t)

 

    print('MCMC sampling by multiprocessing started.') 

    with Pool() as pool:

        Results = pool.map(mm.MCMCs, params)    #  Run the multiprocessing

   

    print('MCMC sampling ended.')

   

    fout_nm = 'regrslts.txt'

    fout = open(fout_nm, 'w')   #   Output text file

 

    #

    #       Posterior distributon and trace of the parameter a

    #

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

    plt.subplot(1,2,1)

    param_a = [[]]*4

    for i in range(4):

        param_a[i] = Results[i][3]

    for i in range(4):

        sb.kdeplot(param_a[i])     #  KDE posteriors for the parameter a

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

    plt.yticks([])

    plt.title('Posterior distributions of a', fontsize = 14)

           

    Rhat_a = mm.calcRhat(param_a)

    plt.subplot(1,2,2)

    for i in range(4):

        plt.plot(range(len(param_a[0])), param_a[i])     #   Traces of sampling of the parameter a

    plt.xlabel('trial', fontsize = 14)

    plt.title(f'Traces of a,  Rhat = {Rhat_a:.2f}', fontsize = 14)

    plt.savefig('FigDistTrace_a.png')

    plt.show()

 

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

    plt.subplot(1,2,1)

    param_b = [[]] * 4

    for i in range(4):

        param_b[i] = Results[i][4]

    for i in range(4):

        sb.kdeplot(param_b[i])

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

    plt.yticks([])

    plt.title('Posterior distributions of b', fontsize = 14)

           

    Rhat_b = mm.calcRhat(param_b)

    plt.subplot(1,2,2)

    for i in range(4):

        plt.plot(range(len(param_b[i])), param_b[i])

    plt.xlabel('trial', fontsize = 14)

    plt.title(f'Traces of b,  Rhat = {Rhat_b:.2f}', fontsize = 14)

    plt.savefig('FigDistTrace_b.png')

    plt.show()

 

    a_whole = []

    b_whole = []

    for i in range(4):

        a_whole += [v for v in param_a[i]]

        b_whole += [v for v in param_b[i]]

       

    a_mode = mm.EstMode(a_whole)

    b_mode = mm.EstMode(b_whole)

    print(a_mode, b_mode)

 

    #       Posterior distribution of a

    sb.kdeplot(a_whole)

    plt.title(f'Posterior distribution of a,  mode = {a_mode:.2f}', fontsize = 16)

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

    plt.yticks([])

    plt.savefig('Fig_a.png')

    plt.show()

 

    #       Posterior distribution of b

    sb.kdeplot(b_whole)

    plt.title(f'Posterior distribution of b,  mode = {b_mode:.2f}', fontsize = 16)

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

    plt.yticks([])

    plt.savefig('Fig_b.png')

    plt.show()

 

    ags_min = np.min(Ags)

    ags_max = np.max(Ags)

 

    def catAgs(v):

        """

                Categorization into 4 categories

        """

        tv = int((v - ags_min) * 4.0 / (ags_max - ags_min))

        if tv > 3:

            tv = 3

        return tv

 

    x_coord_z = np.linspace(ags_min, ags_max, 101)   #  Grid for x-axis

    y_coord_z = np.empty(len(x_coord_z))

    for i, v in enumerate(x_coord_z):

        d = np.exp(a_mode + b_mode * v)              #  d' predicted by the regression model

        y_coord_z[i] = mm.f_Pc(d)                    #  Predicted probability

    plt.plot(x_coord_z, y_coord_z, label = 'Model')

    ags_count = np.zeros(4)

    r_counts = np.zeros(4)       #   Number of responses in each of 4 categories

    rc_counts = np.zeros(4)      #   Number of crrect responses in each of 4 categories

    for rv, av in zip(Rs, Ags):

        ic = catAgs(av)

        r_counts[ic] += 1

        if rv == 1:

            rc_counts[ic] += 1

    ave_rc = rc_counts / r_counts         #  proportions of crrect responses in each category

    ave_ags = np.empty(4)                 #  means of the independent variable in each category

    step_ags = (ags_max - ags_min) / 8.0  #  half of the category width

    for i in range(4):

        ave_ags[i] = ags_min + step_ags + i * step_ags * 2   #  mid points of the categories

    plt.plot(ave_ags, ave_rc, 'o', label = 'data')

    #

    #     Conversion of standardized scores back to raw scores

    #

    x_ags = [(ags_min + i*step_ags*2) for i in range(5)]   #   positions of the raw score labels

    x_ags = np.array(x_ags)

    x_ags_raw = x_ags * Ags_std + Ags_mean

    x_ags_lbl = [f"{int(v)}" for v in x_ags_raw]                #   labels for raw scores

    plt.xticks(x_ags, x_ags_lbl)

    plt.xlabel('Independent variable', fontsize = 14)

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

    plt.title('Regression Curve of Probability', fontsize = 16)

    plt.legend()

    plt.savefig('Fig_regression_p.png')

    plt.show()

   

    fout.write('\nRegression model of Pc on the standardized independent variable Z:\n')

    fout.write('P(R = 1) = Pc(exp({0:.2f} + {1:.2f} * Z))\n'.format(a_mode, b_mode))

   

    a_raw = a_mode - b_mode * Ags_mean / Ags_std

    b_raw = b_mode / Ags_std

    fout.write('\nRegression model of Pc on the raw independent variable X:\n')

    fout.write('P(R = 1) = Pc(exp({0:.2f} + {1:.2f} * X))\n'.format(a_raw, b_raw))

   

    x_coord_z = np.linspace(ags_min, ags_max, 101)   #  Grid for x-axis

    y_coord_d = np.empty(len(x_coord_z))

    for i, v in enumerate(x_coord_z):

        y_coord_d[i] = np.exp(a_mode + b_mode * v)              #  d' predicted by the regression model

    plt.plot(x_coord_z, y_coord_d)

    plt.xticks(x_ags, x_ags_lbl)

    plt.xlabel('Independent variable', fontsize = 14)

    plt.ylabel("d'", fontsize = 16)

    plt.title("Regression curve of d'", fontsize = 16)

    plt.savefig('Fig_regression_d.png')

    plt.show()

 

    fout.write("\nRegression model of d' on the standardized independent variable Z:\n")

    fout.write("d' = exp({0:.2f} + {1:.2f} * Z)\n".format(a_mode, b_mode))

 

    fout.write("\nRegression model of d' on the raw independent variable X:\n")

    fout.write("d' = exp({0:.2f} + {1:.2f} * X)\n".format(a_raw, b_raw))

   

    fout.close()

    print('\n', fout_nm, 'was saved.\n')

 

 

 

Up