恒常法（累積相対度数分布の応用例）

「講義で使える統計素材」シリーズ．今回は，恒常法の説明に必要な図．累積相対度数分布がどのように使用されるのかの例を示すために，弁別閾を調べる恒常法の理論を説明する図を作成しました．

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
import scipy.special as sp
from scipy.optimize import curve_fit

データ¶

In [2]:

data = np.random.normal(170, 10, 100)
print(data)

[178.23506769 170.1704592  169.78269334 170.16925294 181.68554718
 182.19836113 165.3217599  173.93168465 179.08260265 156.89873167
 173.53483302 168.88214726 158.84182373 159.63656197 174.24969275
 173.60929878 183.58539171 172.60132639 187.00964687 175.78404439
 162.6424379  155.84551999 168.2784263  182.79075235 179.44768162
 167.90100493 172.49833271 179.49396604 178.39308857 159.39351474
 163.69230087 167.11013214 161.24449198 176.92952196 168.95127965
 159.6792695  162.16841122 178.06181173 163.1898473  154.93114517
 170.60755762 164.89443918 175.62336281 155.68623881 190.54284347
 186.30781366 159.30964145 178.06951279 177.9710841  146.9442003
 168.38945467 163.27516543 159.6228302  178.695995   170.76631585
 175.86974312 175.93361529 166.43408183 181.7583454  169.21823094
 158.41846824 164.416249   174.541496   183.49841993 171.19171771
 163.40127166 158.29216176 170.47000112 155.35854266 178.00088891
 175.79582669 151.39724417 164.56033404 152.87512497 169.39851671
 152.43474633 159.08325988 166.39205544 174.20996785 159.27190077
 161.35205236 160.80971674 171.1003515  164.75232694 174.21162011
 155.20125793 181.99297242 167.60451108 175.56255822 181.49034339
 169.46381792 164.63820281 164.38397251 145.27589225 163.06470854
 177.42362611 173.87562038 174.01026355 174.19112714 162.88848047]

累積相対度数分布（観測データ）¶

In [3]:

step = 10
bins = range(100, 250, step)
n, bins_ = np.histogram(data, bins)

x = bins[:-1]
y = np.cumsum(n/100)

plt.figure(figsize=(4,4))
plt.ylim(0, 1)
plt.xlim(100, 250)
plt.plot(bins[:-1], y, linestyle='None', marker='.', markersize=10)

plt.savefig('plot_out.svg')

フィッティング結果¶

In [4]:

def err(x, mu,sigma):
    return (1+sp.erf((x-mu)/(np.sqrt(2)*sigma)))/2

popt, pcov = curve_fit(err, x, y, p0=[150, 5])
mu = popt[0]
sigma = popt[1]

print("mu = " + str(mu))
print("sigma = " + str(sigma))

x_ = range(100, 250, 1)
plt.figure(figsize=(4,4))
plt.ylim(0, 1)
plt.xlim(100, 250)
plt.plot(bins[:-1], y, linestyle='None', marker='.', markersize=10)
plt.plot(x_, err(x_, mu, sigma))
plt.plot(x_, st.norm.pdf(x_, mu, sigma))
plt.vlines(mu, ymin=0,ymax=1,linestyle="dashed")
plt.hlines(0.5, xmin=100,xmax=250,linestyle="dashed")

plt.savefig('plot_out.svg')

mu = 158.84691629643382
sigma = 9.941433791491729