potisanのプログラミングメモ

趣味のプログラマーがプログラミング関係で気になったことや調べたことをいつでも忘れられるようにメモするブログです。はてなブログ無料版なので記事の上の方はたぶん広告です。記事中にも広告挿入されるみたいです。

Python&SciPy&matplotlib 正規分布の曲線と下側・上側2.5%領域の塗りつぶし

SciPyとmatplotlibの学習で書いたコードです。実行すると標準正規分布と下側・上側2.5%領域を描画します。

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

# 正規分布で使う平均値と標準偏差
mean_sd = (0, 1)

# グラフのプロット
plt.rcParams["font.family"] = "MS Gothic" #日本語表示対応
fig, ax = plt.subplots()
x = np.linspace(-5, 5, 1000)
#ax.axhline(0, color="gray") #x軸
#ax.axvline(0, color="gray") #y軸
ax.plot(x, norm.pdf(x, *mean_sd), color="black") #正規分布(平均0、標準偏差1)

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm
# 下側2.5%の塗りつぶし
x = np.linspace(*norm.ppf([0.0000001, 0.025], *mean_sd), 100)
ax.fill_between(x, norm.pdf(x, *mean_sd), alpha=0.8)
# 上側2.5%の塗りつぶし
x = np.linspace(*norm.ppf([1 - 0.025, 0.9999999], *mean_sd), 100)
ax.fill_between(x, norm.pdf(x, *mean_sd), alpha=0.8)

ax.set_title("標準正規分布のPDF&下側・上側2.5%領域")
ax.set_xlabel("確率変数")
ax.set_ylabel("確率")
ax.set_ylim(bottom=0) #グラフの描画後に呼び出します。
plt.show()

mean_sd1/2の代わりにPDFやその逆関数ラムダ式で定義することもできます。

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

# 正規分布のPDFとPDFの逆関数
f     = lambda x: norm.pdf(x, 0, 1)
f_inv = lambda x: norm.ppf(x, 0, 1)

# グラフのプロット
plt.rcParams["font.family"] = "MS Gothic" #日本語表示対応
fig, ax = plt.subplots()
x = np.linspace(-5, 5, 1000)
#ax.axhline(0, color="gray") #x軸
#ax.axvline(0, color="gray") #y軸
ax.plot(x, f(x), color="black") #正規分布(平均0、標準偏差1)

# 下側2.5%の塗りつぶし
x = np.linspace(*f_inv([0.0000001, 0.025]), 100)
ax.fill_between(x, f(x), alpha=0.8)
# 上側2.5%の塗りつぶし
x = np.linspace(*f_inv([1 - 0.025, 0.9999999]), 100)
ax.fill_between(x, f(x), alpha=0.8)

ax.set_title("標準正規分布のPDF&下側・上側2.5%領域")
ax.set_xlabel("確率変数")
ax.set_ylabel("確率")
ax.set_ylim(bottom=0) #グラフの描画後に呼び出します。
plt.show()

説明

plt.rcParams["font.family"] = "MS Gothic"

pltmatplotlib.pyplotの慣用的な別名(import matplotlib.pyplot as plt)です。plt.rcParamsにはグラフの描画で使う様々なパラメーターが保管されており、plt.rcParams["font.family"] = "MS Gothic"ではフォントファミリーにMSゴシックを指定します。日本語を表示する場合、このような日本語対応フォントの指定が必要です。

norm.pdf(x, *mean_sd)

タプルと展開演算子を組み合わせることで引数をまとめて管理できます。次のコードはほぼ等価です。

norm.pdf(x, *mean_sd)
norm.pdf(x, mean_sd[0], mean_sd[1])

norm.ppf([0.0000001, 0.025], *mean_sd)

norm.ppfは与えられた面積に対応するパーセント点(確率変数)を返します。理論的にはnorm.ppf([0, 0.025], *mean_sd)としたいですが、norm.ppf(0, *mean_sd) = -Infとなりnp.linspaceに渡せません。一方でnorm.ppf(0.0000001, *mean_sd)は具体的な値なので0.0000001を使います。0.0000001はグラフが崩れない範囲で0に近い適当な値でも構いません。

ライブラリのインストール

SciPyとmatplotlib、それにNumPyをインストールするにはコマンドプロンプトPowerShellで次のコマンドを実行します。

pip install numpy
pip install scipy
pip install matplotlib

公式ドキュメント