1AFC rating task for signal detection theory (SDT)
In One-alternative forced choice (1AFC)
task, one of the two stimulus, e.g., noise only stimulus and noise + signal
stimulus, is presented, and the observer is asked to judge which one was
presented. The sensation of the noise only stimulus is represented by 
, and that of the signal stimulus by 
. They are assumed to have the following distribution
Origin and unit of the sensation are
assumed so that 
. Hence we have
The observer is required to respond with rating about his/her judgment. For example, “Certainly noise only stimulus”, “Probably signal stimulus”, and so on.
In K category rating task, let’s the judgment be represented by integer from 1 to K. “1” represents the most certain judgment of
the noise only, and K the most certain judgment of the stimulus. Corresponding
to the category judgment, the sensation dimension is partitioned by 
category boundaries.
When the sensation is between 
and 
, judgment “k” is
done.
So, we have the following equations:
where 
represents the
cumulative distribution function of the standard normal distribution.
The difference between the means of
sensations and is denoted by 
, that is,
Data of 1AFC rating task can be summarized like Table 1.
|
Tabnle 1 Data of 1AFC rating task |
|||
|
|
Category-1 |
… |
Category-K |
|
Noise only |
N(n1) |
… |
N(nK) |
|
Noise + Signal |
N(s1) |
… |
N(sK) |
N(nk) and N(sk) denote the number of judgment category-k for the noise stimulus and the noise+signal stimulus, respectively.
The likelihood function for the data of Table 1 is given by
(1)
Stan script for eq,(1) is shown in Listing 1.
Python script using the Stan script in Listing 1 is shown in Listing 2.
Files of these scripts are archived in the file sdt1afc_files.zip, which can be freely downloaded and used.
For information about how to install Python and PyStan (Python interface to Stan), check this website.
The file name of the script in Listing 2 is main.py.
To execute the script file main.py, run the following code.
The code was run on Python 3.7 virtual environment of Anaconda on windows 10 using PyStan 2.17.1.0.
For an example using PyStan 3, check the latter part of this website.
python
main.py
When the script is executed, the number of judgment categories K is required to be set as follows.
XXXXX\sdt1afc_files> python main.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
K =
As an example, let’s use the data in Table 2.
|
Table 2. Example of 1AFC rating task data. Number of judgments |
||
|
Judgment |
Stimulus presented |
|
|
Noise stimulus |
Signal stimulus |
|
|
Sure Noise |
6 |
1 |
|
Noise |
18 |
3 |
|
Maybe Noise |
34 |
12 |
|
Uncertain |
31 |
28 |
|
Maybe Signal |
8 |
38 |
|
Signal |
3 |
12 |
|
Sure Signal |
0 |
6 |
In case of Table 2, the number of categories is 7, so input 7.
After setting the value of K, number of judgments of each category to the noise stimulus is set by separating each value by space.
Then, numbers of judgments to the signal stimulus are set.
K = 7
Data for Noise = 6 18 34 31 8 3 0
Data for Signal = 1 3 12 28 38 12 6
After setting the values for the signal stimulus, press the Enter key, then calculation starts.
After MCMC sampling, graph of the posterior
distribution of 
is displayed (Figure 1).
Figure 1
Close the window of the posterior
distribution of , graph of the posterior distribution of 
is displayed (Figure 2).
Figure 2
After closing the window of the posterior
distribution of , it takes some time for calculation.
After the calculation, graph of ROC curve is displayed (Figure 3).
Figure 3
Close the window of ROC curve, graph of posterior distributions of criterion boundaries is displayed (Figure 4).
Figure 4
Close the window, then execution of the script ends.
Execution of the script produces a text file Output .txt, in which some information of the analysis is stored as follows:
Data for Noise:
[6, 18, 34, 31, 8, 3, 0]
Data for Signal:
[1, 3, 12, 28, 38, 12, 6]
mu:
MAP est. = 1.336
med. = 1.317
mean = 1.321
95% CI = [0.954, 1.715]
sigma:
MAP est. = 1.161
med. = 1.175
mean = 1.183
95% CI = [0.916, 1.505]
C1:
MAP est. = -1.579
median = -1.562
mean = -1.566
95% CI = [-1.953, -1.206]
C2:
MAP est. = -0.725
median = -0.720
mean = -0.721
95% CI = [-0.983, -0.465]
C3:
MAP est. = 0.188
median = 0.184
mean = 0.184
95% CI = [-0.052, 0.418]
C4:
MAP est. = 1.184
median = 1.187
mean = 1.188
95% CI = [0.899, 1.484]
C5:
MAP est. = 2.299
median = 2.291
mean = 2.298
95% CI = [1.844, 2.796]
C6:
MAP est. = 3.067
median = 3.147
mean = 3.169
95% CI = [2.491, 3.987]
The script was executed with PyStan 2.17.1.0 installed in Python 3.7 virtual environment of Anaconda in Windows 10.
For information about installation of PyStan, check the following websites:
http://y-okamoto-psy1949.la.coocan.jp/Python/en/PythonAnaconda/
http://y-okamoto-psy1949.la.coocan.jp/Python/en/BeginningPython/
http://y-okamoto-psy1949.la.coocan.jp/Python/en1/PyStanInstalling/
Listing
1. Stan script for 1AFC rating task data (file name: YN_rating.stan)
data {
int<lower = 3> K;
int<lower = 0>
Freq_Noise[K];
int<lower = 0>
Freq_Signal[K];
}
parameters {
real mu;
real<lower = 0.0, upper =
5.0> sigma;
ordered[K-1] c;
}
transformed parameters {
vector[K] theta_n;
vector[K] theta_s;
theta_n[1] =
normal_cdf(c[1], 0.0, 1.0);
theta_n[K] = 1.0 -
normal_cdf(c[K-1], 0.0, 1.0);
for (k in 2:(K-1))
theta_n[k] = normal_cdf(c[k], 0.0,
1.0) - normal_cdf(c[k-1], 0.0, 1.0);
theta_s[1] =
normal_cdf((c[1] - mu) / sigma, 0.0, 1.0);
theta_s[K] = 1.0 -
normal_cdf((c[K-1] - mu) / sigma, 0.0, 1.0);
for (k in 2:(K-1))
theta_s[k]
=normal_cdf((c[k] - mu) / sigma, 0.0, 1.0) -
normal_cdf((c[k-1] - mu) / sigma, 0.0, 1.0);}
model {
mu ~ normal(0.0, 10.0);
sigma ~ uniform(0, 5);
Freq_Noise ~
multinomial(theta_n);
Freq_Signal ~ multinomial(theta_s);
}
Listing
2. Python script using YN_rating.stan in Listing 1 (file name: main.py)
import pystan
import matplotlib.pyplot as plt
import numpy as np
from pystan import StanModel
import seaborn as sb
import scipy.stats as ss
def CalcMAPEst(samples, a = 0.05,
n_points = 2000):
"""
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
"""
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
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx]
fout = open('Output.txt', 'w')
K = int(input('K = '))
s = input('Data for Noise = ')
Data_Noise = []
for v in s.split():
Data_Noise.append(int(v))
s = input('Data for Signal = ')
Data_Signal = []
for v in s.split():
Data_Signal.append(int(v))
fout.write('\nData for Noise:\n')
fout.write(' {}\n'.format(Data_Noise))
fout.write('\nData for Signal:\n')
fout.write(' {}'.format(Data_Signal))
Data = {'K': K, 'Freq_Noise':
Data_Noise, 'Freq_Signal': Data_Signal}
sm = pystan.StanModel(file = 'YN_rating.stan')
fit = sm.sampling(data = Data, iter =
10000, n_jobs = 1)
print(fit)
mu = fit['mu']
sigma = fit['sigma']
mu_mean = mu.mean()
mu_MAP = CalcMAPEst(mu)[0]
mu_q025, mu_med, mu_q975 =
np.percentile(mu, [2.5, 50, 97.5]) # mu_s[int(len(mu_s) / 2)]
sgm_mean = sigma.mean()
sgm_MAP = CalcMAPEst(sigma)[0]
sgm_q025, sgm_med, sgm_q975 =
np.percentile(sigma, [2.5, 50, 97.5])
# sgm_s[int(len(sgm_s) / 2)]
fout.write('\nmu:\n')
fout.write(' MAP est. =
{0:.3f}\n'.format(mu_MAP))
fout.write(' med. =
{0:.3f}\n'.format(mu_med))
fout.write(' mean =
{0:.3f}\n'.format(mu_mean))
fout.write(' 95% CI = [{0:.3f},
{1:.3f}]\n'.format(mu_q025, mu_q975))
fout.write('\nsigma:\n')
fout.write(' MAP est. =
{0:.3f}\n'.format(sgm_MAP))
fout.write(' med. =
{0:.3f}\n'.format(sgm_med))
fout.write(' mean =
{0:.3f}\n'.format(sgm_mean))
fout.write(' 95% CI = [{0:.3f},
{1:.3f}]\n'.format(sgm_q025, sgm_q975))
sb.kdeplot(mu)
plt.title((r'$\mu$') + (r'($\Delta$')
+ 'm) = ' +
'{0:<.3f}(mode)'.format(mu_MAP),
fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize = 14)
plt.xlabel('$\mu$', fontsize = 16)
plt.tight_layout()
plt.savefig('Figure_mu.png')
plt.show()
sb.kdeplot(sigma)
plt.title("$\sigma$ =
{0:<.3f}(mode) ".format(sgm_MAP), fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize = 14)
plt.xlabel('$\sigma$', fontsize = 18)
plt.tight_layout()
plt.savefig('Figure_sgm.png')
plt.show()
c = fit['c']
c_comp = []
for i in range(3):
c_comp.append([])
for v in c:
for j in range(3):
c_comp[j].append(v[j])
z_values = np.arange(-4.0, 4.001, 0.1)
x_values = 1.0 - ss.norm(0.0,
1.0).cdf(z_values)
y_values = 1.0 - ss.norm(mu_MAP,
sgm_MAP).cdf(z_values)
plt.figure(figsize = (5,5))
plt.plot(x_values, y_values, 'b-')
plt.title("ROC-curve",
fontsize = 18)
plt.xlabel("P(False Alarm)",
fontsize = 14)
plt.ylabel("P(Hit)",
fontsize = 14)
cum_freq_n = []
cum_freq_s = []
for k in range(K):
cum_freq_n.append(Data_Noise[K - 1 - k])
cum_freq_s.append(Data_Signal[K
- 1 - k])
for k in range(1, K):
cum_freq_n[k] +=
cum_freq_n[k-1]
cum_freq_s[k] +=
cum_freq_s[k-1]
x_fa = []
y_ht = []
for k in range(K-1):
x_fa.append(cum_freq_n[k] /
cum_freq_n[K-1])
y_ht.append(cum_freq_s[k] /
cum_freq_s[K-1])
plt.plot(x_fa, y_ht, 'ro', label =
'Data')
print('\nI am calculating...')
c_means = np.zeros(K-1)
c_meds = np.zeros(K-1)
c_MAPs = np.zeros(K-1)
c_q025s = np.zeros(K-1)
c_q975s = np.zeros(K-1)
for k in range(K-1):
c_means[k] = c.T[k].mean()
c_MAPs[k] =
CalcMAPEst(c.T[k])[0]
c_q025s[k], c_meds[k],
c_q975s[k] = np.percentile(c.T[k], [2.5, 50, 97.5])
fout.write('\n')
for k in range(K-1):
fout.write('\nC{}:\n'.format(k+1))
fout.write(' MAP est. = {0:.3f}\n'.format(c_MAPs[k]))
fout.write(' median = {0:.3f}\n'.format(c_meds[k]))
fout.write(' mean =
{0:.3f}\n'.format(c_means[k]))
fout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(
c_q025s[k], c_q975s[k]))
x_fa = []
y_ht = []
for k in range(K-1):
x_fa.append(1.0 -
ss.norm().cdf(c_MAPs[k]))
y_ht.append(1.0 -
ss.norm(mu_MAP, sgm_MAP).cdf(c_MAPs[k]))
plt.plot(x_fa, y_ht, 'gs', label =
'Prediction')
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('Figure_ROC.png')
plt.show()
for k in range(K-1):
sb.kdeplot(c.T[k], label =
'$C_{0}$: mode = {1:.3f}'.format(k+1, c_MAPs[k]))
title_str = "Criterions"
plt.title(title_str, fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize = 14)
plt.xlabel('C', fontsize = 14)
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('Figure_Cs.png')
plt.show()
fout.close()
An example of a script using PyStan 3 is shown in Listing 3.
We cannot run a PyStan 3 script directly on Windows. But a PyStan 3 script can be run on WSL (Windows Subsystem for Linux). When we use PyStan 3, the script must be run on WSL, on which we cannot use a graphic library, e.g, matplotlib.
The script in Listing 2 was divided into two parts, one for MCMC sampling using PyStan 3, and the other for analysis of MCMC sample, which uses graphic libraries.
The script for MCMC sampling is shown in Listing 3, which is run on a Python 3.10 virtual environment of Anaconda on WSL on Windows. PyStan 3 is installed on the Python 3.10 virtual environment with the following command. The Stan script used in Listing 3 is the one in Listing 1.
pip install
pystan
.
Analysis of the MCMC sample is done by the script in Listing 4, which can be run on Windows.
Script files are archived in the ZIP file sdt1afc_ps3_files.zip , which can be downloaded freely.
Run the script in Listing 3 as follows
(py310)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/.../sdt1afc_ps3_files$ python mcmc_ps3.py
When the script starts, the number of categories, and frequencies of the categorical responses are required to be set.
When the data in Table 2 is analyzed, input as follows:
K = 7
Data for Noise = 6 18 34 31 8 3 0
Data for Signal = 1 3 12 28 38 12 6
Notice that values are separated by whitespaces.
After the above setting of input data, the building, then the sampling starts.
The MCMC sample is stored in the files infrnc_d.pkl and d_frame.pk by the following codes.
with open('infrnc_d.pkl', 'wb') as f:
pickle.dump(infrnc_d, f)
and
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
The MCMC samples saved in the files are analyzed by the script in Listing 4.
Run the script file anal_ps3.py, then a trace graph is displayed (Figure B1).
Figure B1
Close the form of Figure B1, then some posterior distributions are displayed.
The graph of ROC curve is also displayed after some time of calculation (Figure B2).
Figure B2
After posterior distributions of criterion Cs are displayed, the script ends.
Some statistics of the posterior distributions are saved in Output_Anal.txt.
Contents of the file Output_Anal.txt are as follows:
Data for Noise:
[6, 18, 34, 31, 8, 3, 0]
Data for Signal:
[1, 3, 12, 28, 38, 12, 6]
mu:
MAP est. = 1.293
med. = 1.314
mean = 1.320
95% CI = [0.951, 1.721]
sigma:
MAP est. = 1.159
med. = 1.172
mean = 1.182
95% CI = [0.913, 1.508]
C1:
MAP est. = -1.550
median = -1.563
mean = -1.568
95% CI = [-1.964, -1.204]
C2:
MAP est. = -0.721
median = -0.720
mean = -0.721
95% CI = [-0.984, -0.466]
C3:
MAP est. = 0.187
median = 0.185
mean = 0.185
95% CI = [-0.048, 0.419]
C4:
MAP est. = 1.187
median = 1.186
mean = 1.189
95% CI = [0.899, 1.493]
C5:
MAP est. = 2.275
median = 2.289
mean = 2.298
95% CI = [1.829, 2.801]
C6:
MAP est. = 3.088
median = 3.145
mean = 3.167
95% CI = [2.472, 3.980]
Listing 3. A script using PyStan 3 (mcmc_ps3.py)
import stan
import arviz as az
import pickle
import pandas as pd
pd.set_option('display.max_columns',
None)
fout = open('Output.txt', 'w')
K = int(input('K = '))
s = input('Data for Noise = ')
Data_Noise = []
for v in s.split():
Data_Noise.append(int(v))
s = input('Data for Signal = ')
Data_Signal = []
for v in s.split():
Data_Signal.append(int(v))
with open('noise.pkl', 'wb') as f:
pickle.dump(Data_Noise, f)
with open('signal.pkl', 'wb') as f:
pickle.dump(Data_Signal, f)
fout.write('\nData for Noise:\n')
fout.write(' {}\n'.format(Data_Noise))
fout.write('\nData for Signal:\n')
fout.write(' {}'.format(Data_Signal))
Data = {'K': K, 'Freq_Noise':
Data_Noise, 'Freq_Signal': Data_Signal}
with open('YN_rating.stan', 'r') as
f:
pm = stan.build(f.read(),
data = Data)
fit = pm.sample(num_samples = 10000)
print('fit =\n', fit)
infrnc_d = az.from_pystan(posterior =
fit, posterior_model = pm )
print('\nSummary =\n',
az.summary(infrnc_d))
with open('infrnc_d.pkl', 'wb') as f:
pickle.dump(infrnc_d, f)
d_frame = fit.to_frame()
print(d_frame.keys())
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
fout.close()
print('Output.txt was saved.')
Listing 4. Analysis of MCMC sample (anal_ps3.py)
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sb
import scipy.stats as ss
import pickle
import pandas as pd
pd.set_option('display.max_columns',
None)
def CalcMAPEst(samples, a = 0.05,
n_points = 2000):
"""
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
"""
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
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx]
fout = open('Output_Anal.txt', 'w')
with open('noise.pkl', 'rb') as f:
Data_Noise = pickle.load(f)
with open('signal.pkl', 'rb') as f:
Data_Signal = pickle.load(f)
K = len(Data_Noise)
fout.write('\nData for Noise:\n')
fout.write(' {}\n'.format(Data_Noise))
fout.write('\nData for Signal:\n')
fout.write(' {}'.format(Data_Signal))
fout.write('\n')
with open('infrnc_d.pkl', 'rb') as f:
infrnc_d = pickle.load(f)
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
az.plot_trace(infrnc_d)
plt.tight_layout()
plt.savefig('FigTrace.png')
plt.show()
mu = d_frame['mu']
sigma = d_frame['sigma']
mu_mean = mu.mean()
mu_MAP = CalcMAPEst(mu)[0]
mu_q025, mu_med, mu_q975 =
np.percentile(mu, [2.5, 50, 97.5])
sgm_mean = sigma.mean()
sgm_MAP = CalcMAPEst(sigma)[0]
sgm_q025, sgm_med, sgm_q975 =
np.percentile(sigma, [2.5, 50, 97.5])
fout.write('\nmu:\n')
fout.write(' MAP est. =
{0:.3f}\n'.format(mu_MAP))
fout.write(' med. =
{0:.3f}\n'.format(mu_med))
fout.write(' mean =
{0:.3f}\n'.format(mu_mean))
fout.write(' 95% CI = [{0:.3f},
{1:.3f}]\n'.format(mu_q025, mu_q975))
fout.write('\nsigma:\n')
fout.write(' MAP est. =
{0:.3f}\n'.format(sgm_MAP))
fout.write(' med. =
{0:.3f}\n'.format(sgm_med))
fout.write(' mean = {0:.3f}\n'.format(sgm_mean))
fout.write(' 95% CI = [{0:.3f},
{1:.3f}]\n'.format(sgm_q025, sgm_q975))
sb.kdeplot(mu)
plt.title((r'$\mu$') + (r'($\Delta$')
+ 'm) = ' +
'{0:<.3f}(mode)'.format(mu_MAP),
fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize = 14)
plt.xlabel('$\mu$', fontsize = 16)
plt.tight_layout()
plt.savefig('Figure_mu.png')
plt.show()
sb.kdeplot(sigma)
plt.title("$\sigma$ =
{0:<.3f}(mode) ".format(sgm_MAP), fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize = 14)
plt.xlabel('$\sigma$', fontsize = 18)
plt.tight_layout()
plt.savefig('Figure_sgm.png')
plt.show()
c = [0]*(K-1)
for k in range(K-1):
c[k] = d_frame[f'c.{k+1}']
c = np.array(c).T
print('c =\n', c)
z_values = np.arange(-4.0, 4.001,
0.1)
x_values = 1.0 - ss.norm(0.0,
1.0).cdf(z_values)
y_values = 1.0 - ss.norm(mu_MAP,
sgm_MAP).cdf(z_values)
plt.figure(figsize = (5,5))
plt.plot(x_values, y_values, 'b-')
plt.title("ROC-curve",
fontsize = 18)
plt.xlabel("P(False
Alarm)", fontsize = 14)
plt.ylabel("P(Hit)",
fontsize = 14)
cum_freq_n = []
cum_freq_s = []
for k in range(K):
cum_freq_n.append(Data_Noise[K - 1 - k])
cum_freq_s.append(Data_Signal[K - 1 - k])
for k in range(1, K):
cum_freq_n[k] +=
cum_freq_n[k-1]
cum_freq_s[k] +=
cum_freq_s[k-1]
x_fa = []
y_ht = []
for k in range(K-1):
x_fa.append(cum_freq_n[k] /
cum_freq_n[K-1])
y_ht.append(cum_freq_s[k] /
cum_freq_s[K-1])
plt.plot(x_fa, y_ht, 'ro', label =
'Data')
print('\nI am calculating...')
c_means = np.zeros(K-1)
c_meds = np.zeros(K-1)
c_MAPs = np.zeros(K-1)
c_q025s = np.zeros(K-1)
c_q975s = np.zeros(K-1)
for k in range(K-1):
c_means[k] = c.T[k].mean()
c_MAPs[k] =
CalcMAPEst(c.T[k])[0]
c_q025s[k], c_meds[k],
c_q975s[k] = np.percentile(c.T[k], [2.5, 50, 97.5])
fout.write('\n')
for k in range(K-1):
fout.write('\nC{}:\n'.format(k+1))
fout.write(' MAP est. =
{0:.3f}\n'.format(c_MAPs[k]))
fout.write(' median = {0:.3f}\n'.format(c_meds[k]))
fout.write(' mean =
{0:.3f}\n'.format(c_means[k]))
fout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(
c_q025s[k], c_q975s[k]))
x_fa = []
y_ht = []
for k in range(K-1):
x_fa.append(1.0 -
ss.norm().cdf(c_MAPs[k]))
y_ht.append(1.0 -
ss.norm(mu_MAP, sgm_MAP).cdf(c_MAPs[k]))
plt.plot(x_fa, y_ht, 'gs', label =
'Prediction')
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('Figure_ROC.png')
plt.show()
for k in range(K-1):
sb.kdeplot(c.T[k], label =
'$C_{0}$: mode = {1:.3f}'.format(k+1, c_MAPs[k]))
title_str = "Criterions"
plt.title(title_str, fontsize = 18)
plt.yticks([])
plt.ylabel('Density', fontsize = 14)
plt.xlabel('C', fontsize = 14)
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('Figure_Cs.png')
plt.show()
fout.close()
print('Output_Anal.txt was saved.')