母平均の区間推定

今回は正規分布に従う母集団から抽出した標本について、その母平均を区間推定する式を紹介します。
テキストは東京大学出版会の統計学入門です(11.5.1 正規母集団の母平均、母分散の区間推定)。

まず、母集団の分散$\sigma$は既知で、平均$\mu$を推定したい場合から。
(どんな場合に、分散が既知になるんだろうかという気もしますが、簡単なので先に紹介します。)

標本を $X = (X_1, X_2, \dots, X_n)$とし、標本平均を$\bar X$とします。
すると、$\bar X$は、正規分布 $N(\mu, \sigma^2/n)$に従うので、標準化すると、次のようになります。
$$
P(-Z_{a/2} \leq \frac{\sqrt{n}(\bar X – \mu)}{\sigma} \leq Z_{a/2}) = 1-\alpha.
$$

ここで、$Z_\alpha$は、パーセント点と呼ばれるもので、
標準正規分布$N(0, 1)$において、その点より上側の確率が$100\alpha\%$となる点です。
統計学入門ではP197の10.2 で導入されている記号です。

この式を括弧内の不等式についてとくと次のようになります。
$$
P(\bar X-Z_{a/2} \cdot \frac{\sigma}{\sqrt{n}} \leq \mu \leq \bar X+Z_{a/2} \cdot \frac{\sigma}{\sqrt{n}}) = 1-\alpha.
$$

結果、$\mu$の信頼区間は、次のようになります。
$$
[\bar X-Z_{a/2} \cdot \frac{\sigma}{\sqrt{n}}, \bar X+Z_{a/2} \cdot \frac{\sigma}{\sqrt{n}}].
$$

さて、次は$\sigma^2$が未知の場合です。通常はこちらだと思います。
この時は、標本分散を使うことになります。標本分散$s^2$の式はいつものこれ。
$$
s^2 = \sum_{i=1}^n\frac{(X_i-\bar X)^2}{n-1}.
$$

すると今度は$\sqrt{n}(\bar X-\mu)/s$は、自由度$n-1$の$t$分布に従うので、
$t$分布のパーセント点を使って、次のようにかけます。
$t_{\alpha}(n)$は自由度$n$の$t$分布の$\alpha$パーセント点です。

$$
P(-t_{a/2}(n-1) \leq \frac{\sqrt{n}(\bar X – \mu)}{s} \leq t_{a/2}(n-1)) = 1-\alpha.
$$
これを先ほど同じように整形すると、次のようになります。
$$
P(\bar X-t_{a/2}(n-1) \cdot \frac{s}{\sqrt{n}} \leq \mu \leq \bar X+t_{a/2}(n-1) \cdot \frac{s}{\sqrt{n}}) = 1-\alpha.
$$

よって、母平均$\mu$の信頼係数$1-\alpha$の信頼区間はこのようになります。
$$
[\bar X-t_{a/2}(n-1) \cdot \frac{s}{\sqrt{n}}, \bar X+t_{a/2}(n-1) \cdot \frac{s}{\sqrt{n}}].
$$

点推定と区間推定

東京大学出版会の統計学入門の10〜11章を読み直したので、これから数回の更新で区間推定についてまとめておこうと思います。

確率分布のパラメーター(正規分布であれば平均や標準偏差)を、標本から定めることを
パラメーター(母数)の推定と言います。

その時、パラメーター$\theta$をある一つの値$\hat\theta$で指定する方法点推定(point estimation)と呼びます。
例えば標本の平均をとってそれを母平均の推定値とするようなことです。

それに対して、区間推定(interval estimation)は、真のパラメーターの値が入る確率が
ある値$1-\alpha$以上と保証される区間$[L, U]$を求めるものです。
言い換えると、次の式を満たす$L,U$を求めるものです。
$$
P(L\leq \mu\leq U) \geq 1-\alpha.
$$

点推定の場合は別途誤差の評価を行う必要がありますが、
区間推定は最初からある程度の誤りがあることを認めた推定法と言えます。

Googleアナリティクスで参照元ページのURLを調べる

このブログも200件以上の記事を書いてようやくオーガニックサーチからの訪問以外のランディング、
つまり他のサイトのリンクからの訪問が発生するようになりました。

Googleアナリティクスの [参照元 / メディア] で調べると、 teratail.com / referral からやってきている人がいます。
そこで、teratailのどこにこのブログへのリンクが貼られているのかきになるのですが、
このままではそれがわかりません。

具体的な参照元のURLを知るためには、セカンダリディメンションか、カスタムレポートを使います。
選択するディメンションは [完全なリファラー] です。(APIでは ga:fullReferrer)

これを指定すると、”メディアによっては” リンク元のURLがわかります。
結果、こちらのページに引用していただけてるのを見つけました。

cabochaのPythonバインディングがエラーとなる

ちなみに、Qiitaからも流入があるのですが、
完全なリファラーに表示されるURLもQiitaのトップページになってしまっていて、
こちらはどの記事からの訪問者なのかわかりませんでした。

オーガニックサーチについても検索エンジンのトップページのURLが表示されるので、
「完全な」とついている割には不完全な情報しか取れないようです。

変動係数

なんとなく統計検定2級の過去問を眺めていたら、変動係数というキーワードが登場してきました。
語感や説明なく使われている様子から基礎的な用語な感じがしたのですが聞いた覚えがありません。
それで、東大の統計学入門(赤本)にも載ってないんじゃないかと思って調べたのですが、バッチリ載っていました。
自分が覚えていなかっただけのようです。(38ページに登場します。)
良い機会なのでここで定義を紹介しようと思います。

変動係数(Coefficient of Variation)とは、標準偏差を期待値で割ったものです。

本に載っている用途をそのまま紹介しますが、
これは二つの分布の中心(期待値)の大きさが著しく異なり、標準偏差を使って散らばり具合を比較できない場合などに使います。
たとえば、県の1人当たり県民所得の平均と標準偏差が次の値だったとします。
1965年は平均26.6万円、標準偏差は7.5万円で、
1975年は平均117.5万円、標準偏差は23.8万円。

これで、標準偏差が約3倍になっているから地域間の格差が広がっているか?という問いを考えると、
平均も約4.5倍になっているのでそう単純に比較できない、という場合に変動係数で比較するという方法が使えます。

1965年の変動係数は$7.5/26.6 = 0.28$で、1975年のそれは$23.8/117.5 = 0.2$なので、相対的な地域間所得格差は小さくなっていると判断できます。

以上のように紹介してみましたが、結構用途が限られる指標だとも感じました。
確率変数$X$に定数$a$を足した分布を考えると値が変わってしまいますし、
平均が$0$の場合は定義できません。
統計学入門の以降のページにもほぼ登場していないようです。

とはいえ、時系列データなどを扱う時、時代によって平均値が大きく変わってしまった場合の散らばりの比較などには
使えそうなので覚えておこうと思います。

Pythonの数値を2進法、8進法、16進法の表記に変換する

以前の記事で、 pythonで、2進法/8進法/16進法で数値を定義する というのを書きました。
今回はその逆に、10進法の数値を2進法、8進法、16進法での表記に変換します。
(なお、結果的にデータ型は文字列になってしまいます。)

これには、組み込み関数として実装されている bin(), oct(), hex()を使います。


num = 23456
print(bin(num))
# 0b101101110100000
print(oct(num))
# 0o55640
print(hex(num))
# 0x5ba0

先頭の0b等の有無の調整や、16進法でアルファベットを大文字/小文字のどちらで表記するかの制御なども行いたい場合、
format関数が使えます。
ずらっと並べると次のような感じ。


print(format(num, "b"))
# 101101110100000
print(format(num, "#b"))
# 0b101101110100000
print(format(num, "o"))
# 55640
print(format(num, "#o"))
# 0o55640
print(format(num, "x"))
# 5ba0
print(format(num, "#x"))
# 0x5ba0
print(format(num, "X"))
# 5BA0
print(format(num, "#X"))
# 0X5BA0

16進法(hex)の場合に、hではなくxを使うところに注意が必要です。

globで手軽にファイル名の一覧を取得する

特定のディレクトリの配下にあるファイルの一覧が欲しい場面というのはよくあります。
サブディレクトリの探索等少々高度なことをする時はもっと違うライブラリを使ったほうがいいのですが、
特定ディレクトリ直下の特定のパターンのファイル名のファイルの一覧を取得する時などは、
glob を使うと便利です。

ドキュメント: glob — Unix 形式のパス名のパターン展開

次の例はカレントディレクトリ直下のテキストファイルをリストアップしたもの。


import glob
glob.glob("./*.txt")
# ['./text1.txt', './text2.txt']

ご覧の通り、ワイルドカードとして*が使えます。また、?も使えます。
パスの指定は相対パス、絶対パスの両方に対応していて、イメージ通りの挙動をしてくれるのでとても手軽です。

scipyで定積分

タイトルの通り、scipyで定積分を計算する方法の紹介です。

とりあえず今回は $\frac{4}{1+x^2}$ を 区間$[0,1]$で積分しみてみましょう。
なお、この答えは$\pi$になります。

scipyで定積分をする時は integrate モジュールに定義されている、quad という関数を使います。
ドキュメント: scipy.integrate.quad


import scipy.integrate as integrate


def f(x):
    return 4/(1+x**2)


print(integrate.quad(f, 0, 1))
# (3.1415926535897936, 3.4878684980086326e-14)

ご覧通り、結果はタプルで戻ってきます。
一つ目の要素が積分の答えであり、確かに円周率ぽい値になっています。
そして、二つ目の要素は、誤差の推定値です。

これはscipyが代数的に積分を計算しているのではなく、
数値計算で結果を返しているため、どうしても誤差が発生するためです。

NetworkXの辺の書式を設定する

前回の記事がNetworkXの頂点の書式についての記事だったので今回は辺の書式です。

ドキュメントは同じ場所です。
ドキュメント: draw_networkx

有向グラフの場合は矢印の形などもう少し項目が多いのですが、
無向グラフの場合は、太さと色とスタイル(単線、破線、ドットなど)が設定項目になります。

このうち、太さ(width) と 色(edge_color) は単一の値だけではなく、
辺の数と同じ長さの配列で指定することができます。

配列で指定した値は、 G.edge で取得できる辺の一覧の順番と対応するように設定されます。
(初めてNetworkXをいじった時はこの辺りをしらず、思った設定ができずにいました。)

細かいですがドキュメントにデフォルトの色は’r’ (赤)と書いてありますが、これはおそらくドキュメントの誤りです。
(何も指定しないと黒になります。)

edge_color (color string, or array of floats (default=’r’))

styleに指定できるのは次の4種類です。(スタイルの指定は辺別ではなく、全体で一つ指定する必要があります。)
solid|dashed|dotted|dashdot

それでは適当なグラフで実行してみましょう。


import networkx as nx
import matplotlib.pyplot as plt

# グラフを生成
G = nx.Graph()
G.add_cycle(range(5))

# 辺の順番を確認
print(G.edges)
# [(0, 1), (0, 4), (1, 2), (2, 3), (3, 4)]

# 可視化
nx.draw_networkx(
    G,
    edge_color=["r", "g", "m", "c", "y"],
    width=[1, 2, 3, 4, 5],
    style= "dashdot",
)
plt.show()

出力がこちら。

0と1を結んでるのが一番細い赤い線で、次に0と4を結んでるが緑の線で、
と G.edges の値と対応して設定が反映されました。

NetworkXの頂点の書式を設定する

NetworkXのグラフを可視化するときに、もっと見た目を変えたいなと思うことは多々あると思います。
そのようなニーズに応えるために頂点と辺のそれぞれにいくつかのオプションが用意されています。
今回はその中から、頂点の書式に関するものをいくつか紹介します。

ドキュメント: draw_networkx

このページの node_xxx の形式で定義されているパラメーターがそれです。
この中から、node_size、node_color、node_shapeを使って、頂点の色や大きさを変えてプロットしてみましょう。
なお、node_shapeはグラフに対して1つですが、node_sizeとnode_colorは頂点ごとに設定を変えることができます。
ついでに font_family も設定して日本語のフォントを表示してみましょう。


# グラフの生成
G = nx.Graph()
for n in "あいうえお":
    G.add_node(n)

# 頂点の順番の確認。(node_size や node_color はこの順番で指定する)
print(G.nodes)
# ['あ', 'い', 'う', 'え', 'お']

nx.draw_networkx(
    G,
    node_shape="s",
    node_size=[100, 200, 300, 400, 500],
    node_color=["r", "g", "m", "c", "y"],
    font_family="IPAexGothic"
)
plt.show()

出力がこちら。

node_size はデフォルトが 300であることに注意しましよう。
node_shape には ‘so^>v<dph8’ のいずれかの文字が指定できます。

matplotlibでTableau風の色を使う

滅多なことでは使う機会がないと思うのですが、matplotlibで使える色名の中に、
Tableauのカラーパレット(Tableau 10)の色を指定するものがあるのを見つけたので紹介します。

参考:matplotlib.colors
こちらのページにひっそりと次の記載があります。

one of the Tableau Colors from the ‘T10’ categorical palette (the default color cycle): {‘tab:blue’, ‘tab:orange’, ‘tab:green’, ‘tab:red’, ‘tab:purple’, ‘tab:brown’, ‘tab:pink’, ‘tab:gray’, ‘tab:olive’, ‘tab:cyan’} (case-insensitive);

定数も用意されているのでそれを確認してみましょう。


import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

mcolors.TABLEAU_COLORS
'''
OrderedDict([('tab:blue', '#1f77b4'),
             ('tab:orange', '#ff7f0e'),
             ('tab:green', '#2ca02c'),
             ('tab:red', '#d62728'),
             ('tab:purple', '#9467bd'),
             ('tab:brown', '#8c564b'),
             ('tab:pink', '#e377c2'),
             ('tab:gray', '#7f7f7f'),
             ('tab:olive', '#bcbd22'),
             ('tab:cyan', '#17becf')])
'''

matplotlibで色を指定する部分に”tab:blue”と入れてやればいつも見慣れたTableauの青が表示されます。
(“tab:olive”は何か違うような気がするのですが僕の環境のせいでしょうか)

せっかくなので、10本の棒グラフを用意して使ってみましょう。


x = range(1, 11)
y = range(10, 0, -1)
fig = plt.figure()
ax = fig.add_subplot(111, title="Tableau Colors")
ax.bar(x, y, color=mcolors.TABLEAU_COLORS)
plt.show()

出力はこちら。