物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

Python (numpy+scipy)でプロファイルヒストグラムを描画する

プロファイルヒストグラム(Profile histograms ROOT: TProfile Class Reference を参照) は、2次元データをあるX軸の区画ごとに、Y軸のデータの平均値と標準偏差などを描画したものである。

matplotlib - Plotting profile hitstograms in python - Stack Overflowソースコードを参考に、実際に描画してみる

import scipy.stats
import numpy as np
import matplotlib.pyplot as plt

x = np.random.rand(10000)
y = x + scipy.stats.norm(0, 0.2).rvs(10000)

plt.scatter(x, y, marker=".", s=1)
plt.ylim(-0.5, 1.5)
plt.show()

means_result = scipy.stats.binned_statistic(x, [y, y**2], bins=50, range=[0,1], statistic='mean') #binごとの統計を計算

means, means2 = means_result.statistic # yとy**2の平均値
standard_deviations = np.sqrt(means2 - means**2) # 標準偏差を計算

bin_edges = means_result.bin_edges # binの両端
bin_centers = (bin_edges[:-1] + bin_edges[1:])/2.0 # binの中心を計算

plt.errorbar(x=bin_centers, y=means, yerr=standard_deviations, linestyle='none', marker='.')
plt.ylim(-0.5, 1.5)
plt.show()

f:id:onsanai:20210812131405p:plain

f:id:onsanai:20210812131418p:plain

絵を見るとプロファイルヒストグラムが何かは一目瞭然だろう。

検出効率のような2値を扱う場合も書いてみる。0以上は検出できた=1、0未満は検出できなかった=0とする。

import scipy.stats
import numpy as np
import matplotlib.pyplot as plt

x = np.random.rand(1000)
y = x + scipy.stats.norm(0, 0.2).rvs(1000)
y = np.array([(1 if yi>=0 else 0) for yi in y]) #何かを検出したかどうかを2値で表す
one = np.ones(1000)

plt.scatter(x, y, marker=".", s=1)
plt.ylim(-0.5, 1.5)
plt.show()

means_result = scipy.stats.binned_statistic(x, y, bins=50, range=[0,1], statistic='mean')
means = means_result.statistic #検出効率
ones = scipy.stats.binned_statistic(x, one, bins=50, range=[0,1], statistic='sum').statistic #binに入った数

standard_deviations = np.sqrt(means * (1-means) / ones) # 二項分布の誤差

bin_edges = means_result.bin_edges # binの両端
bin_centers = (bin_edges[:-1] + bin_edges[1:])/2.0 # binの中心を計算

plt.errorbar(x=bin_centers, y=means, yerr=standard_deviations, linestyle='none', marker='.')
plt.ylim(-0.5, 1.5)
plt.show()

f:id:onsanai:20210812154916p:plain

f:id:onsanai:20210812154933p:plain

できた。