2021-12-06

[JavaScript][TypeScript] 配列の直積（すべての順列）を得る関数productを作る

はじめに

Pythonの場合

TypeScriptで書く

JavaScriptの場合

まとめ

はじめに

JavaScript/TypeScriptで複数の配列の直積（すべての順列）が欲しいという場面があった。

Pythonなら組み込み関数productで解決できるのだが、JavaScript/TypeScriptではそのような関数を見つけられなかったので、自作してみることにした。

Pythonの場合

Pythonでは、標準ライブラリitertoolsのproduct関数を使うと、複数の配列の直積を得ることができる。

返り値はイテレータなので、リストに変換して表示させれば、目的の配列が得られていることが分かる。

gist.github.com
出力結果

[(1, 4), (1, 5), (2, 4), (2, 5), (3, 4), (3, 5)]

product関数を使えば、単純に配列の直積（順列）を得ることもできるし、複数のfor文を一つに減らしてコードの見通しを良くすることもできる。

cyanatlas.hatenablog.com

TypeScriptで書く

引数となる配列の数が分からないため、ここでは再帰を使って書くことにした。

大きな配列を代入する場合はメモリに注意されたい。

gist.github.com

上のコードをutils.tsというファイル名で保存し、同じディレクトリで次のコードを実行すると、目的の配列が得られていることが確認できる。

gist.github.com
出力結果

[ [ 1, 4 ], [ 1, 5 ], [ 2, 4 ], [ 2, 5 ], [ 3, 4 ], [ 3, 5 ] ]

まとめ

JavaScript/TypeScriptでもPythonで提供されているような直積を求める関数productを書くことができた。
引数となる配列の数が不明なため、再帰を使って書いた。
ここでは引数を二次元配列としたが、残余引数構文を使った書き方もありそう。
使用する際は、組み合わせ爆発や再帰処理に伴うメモリ不足に注意が必要。

2020-11-10

信号検出理論に関するノート

信号検出理論統計

はじめに

問題の設定と数理モデル

弁別力の指標 d'（Sensitivity index）

受信者操作特性（Receiver operating characteristic: ROC）

参考文献

はじめに

信号検出理論（Signal Detection Theory）という理論に興味を持ったので、調べつつざっくりまとめてみる。

もともとはノイズの中から真の信号を取り出すための工学理論らしいのだが、その後に心理学や臨床試験、品質管理などの分野に応用されていったようだ。

en.wikipedia.org

問題の設定と数理モデル

ある人が二種類の異なる刺激を受けとり、それを識別しなければならないという状況を考える。例えば、真水または食塩水を渡され、それを舐めてみて真水か食塩水か答えるという状況である。

結果は四パターンに分類される。

真水を真水と答える：Correct rejection
真水を食塩水と答える：False alarm
食塩水を真水と答える：Miss
食塩水を食塩水と答える：Hits

受け取る刺激は、本来の刺激にノイズが加わった形で与えられるとしよう。例えば、食塩水を舐めたときの刺激の強さの確率分布は、ある平均と分散を持つ正規分布に従っていると仮定する。

f:id:cyanatlas:20201110153442j:plain

もし二種類の刺激の分散が大きければ、刺激の確率分布は重複することになり、二種類の刺激を完全に正確に識別することはできなくなる。

例えば上の図では、刺激の強さが2.5のとき、真水（青）である確率は50%、食塩水（橙）である確率は50%となり、真水か食塩水かを完全に識別することはできない。

そこで、完全な識別を諦め、ある閾値を定めて、閾値より大きいか小さいかで刺激を区別するという方法をとることにしよう。

f:id:cyanatlas:20201110153804j:plain

上の図では、閾値を3.0（赤）に設定している。すなわち、3.0より小さければ真水（青）、3.0より大きければ食塩水（橙）と判定する、という識別ルールを採用している。

このときの刺激の種類別の正答率を計算してみると以下のようになる。

真水が与えられたとき
- 真水と答える（Correct rejection）：98%
- 食塩水と答える（False alarm）：2%
食塩水が与えられたとき
- 真水と答える（Miss）：9%
- 食塩水と答える（Hits）：91%

真水が与えられたときの正答率は98%、食塩水が与えられたときの正答率は91%となる。刺激が強くても真水と判定するような識別ルールを採用しているため、真水の正答率は高い一方で、食塩水の正答率が低いという結果になっている。

弁別力の指標 d'（Sensitivity index）

上記の問題設定と数理モデルを踏まえて、現実の問題を考えてみよう。

現実には、二種類の刺激がどれほど違うのかは分からない。むしろ、どれほど違うのかを調べたいという場合がある。このようなときに、上記の問題設定と数理モデルを用いて、二種類の刺激の違いの大きさを測定する。

二種類の刺激の違いを、平均値の差で評価することにしよう。これをdプライムと呼ぶ。先程の図では、正規分布の平均値の差は5.0であったので、dプライムは5.0だったということになる。

ではどうやってdプライムを計算するか。先程はdプライムと分散、識別閾値を与えた上で刺激別の正答率を計算した。そこでこの手順を逆に行えば、dプライムを求めることができる。

具体的には、二種類の刺激を被験者に与えて、それぞれの正答率を測定する。その正答率からdプライムを計算する。

具体例で考えてみる。真水の正答率が80%、食塩水が90%だったとしよう。いずれの刺激も標準正規分布（分散が1の正規分布）に従うとすれば、80パーセンタイルと90パーセンタイルを足し合わせれば、dプライムが求まる。この例では、dプライムは2.12、識別水準は0.84と計算される。

en.wikipedia.org

受信者操作特性（Receiver operating characteristic: ROC）

上で説明したdプライムは直観的であるものの、刺激の確率分布が正規分布であり、かつ二種類の刺激は等分散であるという仮定を含んでいた。しかし、このような条件がいつでも満たされるとは考えにくい。

そこで、dプライムの代わりにROCカーブというものを用いて刺激の違いの大きさを評価する。ROCカーブは、false alarmの割合を横軸に、hitの割合を縦軸にとって描いた曲線のことで、この曲線が上に凸であればあるほど刺激の違いが大きいことを意味する。

f:id:cyanatlas:20201110170726j:plain

二種類の刺激に対する識別閾値を変えると、それに従ってFalse alarm率とHits率が変化する。識別閾値を負の無限大から正の無限大にわたって変化させれば、ROCカーブを描くことができる。

ja.wikipedia.org

en.wikipedia.org

参考文献

畑江敬子 (1993)「信号検出理論の官能検査への応用」調理科学 Vol. 26, No. 1

www.jstage.jst.go.jp

石田翼「信号検出理論の指標をめぐって」 http://www5e.biglobe.ne.jp/~tbs-i/psy/tsd/tsd.html

www5e.biglobe.ne.jp

2020-10-21

異なるN人のグループ分けの総数が知りたい

数学 Python itertools

はじめに

例：4人のグループ分けの総数

Pythonで計算してみる

n=4の場合

n=8までの場合

計算方法の難点

ベル数を使って計算する

はじめに

識別可能なN人を任意のグループに分けるとき、その分け方の総数を知りたい。

はじめにPythonを使って力業で数え上げて計算してみた。次に、ベル数を使った計算を行った。

例：4人のグループ分けの総数

a, b, c, dの4人がいたとしよう。分け方は以下の15通り（のはず）。

グループの分け方	メンバーの分け方
oooo	abcd
ooo \| o	abc \| d
	abd \| c
	acd \| b
	bcd \| a
oo \| oo	ab \| cd
	ac \| bd
	ad \| bc
oo\|o\|o	ab\|c\|d
	ac\|b\|d
	ad\|b\|c
	bc\|a\|d
	bd\|a\|c
	cd\|a\|b
o\|o\|o\|o	a\|b\|c\|d

Pythonで計算してみる

n=4の場合

以下の手順でグループ分けの総数を求めた。

人数分のグループ（空リストn個からなるリスト）を用意する。
一人ずついずれかのグループに配属させる（リスト内リストにappendする）。
全員がいずれかのグループに配属したら、グル－プ分けの状態を整頓する（sorted関数でリスト内要素を並び替える）。
得られたグループ分けを記録する（group_listに追加する）。
これ（1~4）をすべての配属パターンについて行う。
group_listから重複を取り除く（一度set関数で集合にする）。
group_listの要素数を数える（len関数）。

簡単に計算するために標準ライブラリitertoolsに含まれる関数を使う。
docs.python.org

import string
from itertools import product


def main(n):
    assert n <= len(string.ascii_letters), ValueError
    menbers = [s for s in string.ascii_letters[0:n]]
    group_list = []
    assignment_list = product(range(n), repeat=n)
    for assignment in assignment_list:
        groups = [[] for _ in range(n)]
        for e, i in enumerate(assignment):
            groups[i].append(menbers[e])
        groups = '|'.join(sorted([''.join(agroup) for agroup in groups]))
        group_list.append(groups)
    group_list = sorted(list(set(group_list)))
    print(len(group_list))
    print(group_list)
    return


if __name__ == '__main__':
    main(n=4)

実行すると次のようになる。

15
['a|b|c|d', '|ab|c|d', '|ac|b|d', '|ad|b|c', '|a|bc|d', '|a|bd|c', '|a|b|cd', '||abc|d', '||abd|c', '||ab|cd', '||acd|b', '||ac|bd', 
'||ad|bc', '||a|bcd', '|||abcd']

全部で15通りなので、手計算の結果と一致する。

n=8までの場合

N人の場合を計算してみる。値が大きくなると計算時間が非常に大きくなるので、n=8まででテスト。

import string
from itertools import product
import pandas as pd


def calc(n):
    assert n <= len(string.ascii_letters), ValueError
    menbers = [s for s in string.ascii_letters[0:n]]
    group_list = []
    assignment_list = product(range(n), repeat=n)
    for assignment in assignment_list:
        groups = [[] for _ in range(n)]
        for e, i in enumerate(assignment):
            groups[i].append(menbers[e])
        groups = '|'.join(sorted([''.join(agroup) for agroup in groups]))
        group_list.append(groups)
    group_list = sorted(list(set(group_list)))
    print(len(group_list))
    print(group_list)
    return len(group_list)


def main(n):
    df = []
    for k in range(n):
        df.append([k, calc(k)])
    df = pd.DataFrame(df, columns=['n_persons', 'n_groups'])
    df.to_csv('result.csv', index=False)
    return


if __name__ == '__main__':
    main(n=9)

力業でやっているのでまあまあな時間がかかる。出力されるCSVファイルを結果として示す。

n_persons	n_groups
0	1
1	1
2	2
3	5
4	15
5	52
6	203
7	877
8	4140

人数が増えるとグループ分けの総数は爆発的に増える。

f:id:cyanatlas:20201021113213j:plain

計算方法の難点

上で示したPythonコードでは、 $n^n$ 通りのグループ分けを計算した後で、重複を取り除くという方法をとっている。

そのため、最初に行うグループ分けの数が、 $n = 9$ のときでさえ3.8億通りであり、 $n$ が大きくなると計算負荷は爆発的に増加してしまう。

手元のノートパソコンでは、 $n \geq 9$ のときは応答が返ってこなくなってしまったので、計算できなかった。使用するメモリ量を減らすような書き方をしないといけないかもしれない。

ベル数を使って計算する

異なるn個のものをグループに分けるやり方の総数は、ベル数と呼ばれている。

実際、wikiに記載されているベル数を見ると

1, 1, 2, 5, 15, 52, 203, 877, 4140, ...

となっていて、先にPythonで数え上げたときの結果と一致している。

ja.wikipedia.org

グループ分けの一覧ではなく、グループ分けの総数だけが知りたいのであれば、数え上げまでせずにベル数を計算してやれば十分だ。

from scipy.special import comb
import pandas as pd
import matplotlib.pyplot as plt


def bell(n):
    if n == 0:
        return 1
    elif n == 1:
        return 1
    else:
        res = 0
        for k in range(n):
            res += comb(n - 1, k, exact=True) * bell(k)
        return res


def main(n):
    df = []
    for i in range(n + 1):
        df.append([i, bell(i)])
    df = pd.DataFrame(df, columns=['n', 'bell'])
    df.to_csv('bell.csv')

    fig, ax = plt.subplots()
    ax.plot(df['n'], df['bell'], marker='o')
    ax.set_xlabel('n')
    ax.set_ylabel('Bell number')
    fig.tight_layout()
    fig.savefig('bell.jpg', dpi=300)
    return


if __name__ == '__main__':
    main(n=10)

f:id:cyanatlas:20201021122511j:plain

数え上げる方法よりも早く計算できたため、より大きなnに対しても計算が可能になった。

2020-10-20

[統計] 統計でおすすめの本（主に回帰）

統計 R 書評

はじめに

諸注意

対象者について

統計の基本～単回帰

統計学入門

一般化線形モデル（GLM）

データ解析のための統計モデリング入門

一般化線形モデル入門

自然科学の統計学

R

アート・オブ・Rプログラミング

はじめに

統計学を勉強するために読んだ本の中から、オススメの本をいくつか紹介しつつまとめてみる。

統計学と一口に言っても、トピックは多岐にわたるので、この記事では回帰に関する記述が豊富な書籍に絞った。

諸注意

本稿の筆者は統計学の専門家ではありません。
基本的には独立なデータに対する回帰の話題に絞っています。例えば時系列解析の書籍などは本稿では扱っていません。

対象者について

本の良し悪しは絶対的なものでなく、往々にして目的依存であると考えている。そこで、本を紹介する前に、私がどのような目的で統計学の本を読んでいるか説明しておく。

統計手法がどのような目的を持つか知りたい（どんなタイプのデータから何を導けるか）
統計手法がどのような数学的理論に基づいているのか知りたい（「このコマンドを打てばP値が出ます」では不満。とは言え、数理統計学の本を読んで定理と証明を厳密に知りたいというほどではない）
得られた解析結果が統計学的にどう解釈されるか知りたい（有意差がある・ない、AIC最小ってつまりどういう"意味"？）

私と似たような目的を持って統計学の本を探している人には、以下で紹介する本が役に立つかもしれない。

統計の基本～単回帰

統計学入門

統計学入門 (基礎統計学Ⅰ)

発売日: 1991/07/09
メディア: 単行本

オススメ度　★★★☆☆

良いところ

検定や区間推定の細かいところまで詳しく書かれている。例えば、t統計量を用いた区間推定をどのように行うのか、その計算の詳細まで書かれており、付録の表を使えば手計算でも区間推定が可能なほど詳しい。
統計学での考え方をはっきり書いている。例えば、区間推定で得られた区間がどういう意味なのか、詳しく説明されている。
個人的には、第1章で紹介されている統計学の歴史がとても興味深かった。統計学が一見するとさまざまな手法の寄せ集めにも見える理由に納得がいった。

良くないところ

書名に「入門」とあるものの、書き方が固く読みづらい。もちろん、固い本に慣れている人ならあまり支障はないかもしれない。
話の行く末が見えづらい。本の前半は用語の紹介や古典的手法の説明が続くので（非常に大事なことが書かれているものの）、退屈な印象を受ける。
基本的かつ古典的な手法しか載っていない。最終章まで読んでも、一変数の直線回帰（単回帰）までしかできるようにならない。もちろん、新しい手法を理解するために古典的な手法を抑えておくべきではある。

総評
統計学の考え方をちゃんと学びたい人には、自信を持ってオススメできる。この本が読める人には1冊目としてオススメしたいが、固すぎて読めないという人も多いと思う。

一般化線形モデル（GLM）

一般化線形モデルは正規線形モデル（直線回帰）の拡張。正規線形モデルを勉強したいという場合でも、その拡張である一般化線形モデルを勉強しておき、正規線形モデルをその特殊な例とみなす方が、見通しが良いように感じる。また、分散分析や共分散分析は、正規線形モデルに基づく分析手法であるため、これらの手法を学ぶ場合にも一般化線形モデルに触れておいて損はないと思う。本によっては分散分析も扱っている。

データ解析のための統計モデリング入門

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者:久保拓弥
発売日: 2012/05/19
メディア: 単行本

オススメ度　★★★★☆

統計を扱ったインターネットの記事ではよく見かける本で、緑本と呼ばれて親しまれているようだ。

良いところ

前半でGLMとGLMM、後半でベイズ階層モデルが扱われており、統計モデルを概観することができる。とくに前半のGLMの話は非常に分かりやすい。
実際に動くRのコードも載っており、自分で手を動かしながら学ぶことができる。
AICについてのまとまった記述がある。AICについては、非専門家向けに日本語で書かれた解説は少ないように思うので、これは貴重だと思う。

良くないところ

数学的な話・理論的な話は省略されている。その辺りは無理にこの本で理解するよりも、他の本を読んだ方が良い。例えば下で紹介するDobson『一般化線形モデル入門』など。
統計モデルを使った回帰の話に特化している。例えば相関に関する話題はない。まあこの点は良くないところというわけではないが･･･

総評
とても読みやすい本なので、GLMを勉強する際の1冊目としてオススメできる。

一般化線形モデル入門

一般化線形モデル入門原著第2版

作者:Annette J.Dobson
発売日: 2008/09/08
メディア: 単行本

オススメ度　★★★★★

数式を使いながら一般化線形モデルを解説している本。

良いところ

一般化線形モデルを数式を使ってちゃんと定式化し、最尤推定なども数式で示している。式展開は丁寧で分かりやすい。本文の説明が丁寧なので、途中の式を読み飛ばしながらでも十分内容を理解できる。
指数型分布族、デザイン行列（計画行列）、連結関数（リンク関数）、スコア統計量、情報量といった用語は三章までで分かりやすく説明されている。
最終章には、相関のあるデータに対する回帰手法が紹介されている。

良くないところ

他の本に比べると数式が多いので（とは言え数理統計の本ほどではないが）、とっつきにくさは感じる。ただ、言葉の説明が曖昧で分かりにくいところは、数式を見に行けば厳密な意味が理解できるという利点はある。

総評
数式が多くてとっつきにくいかもしれないが、とても良い本なので是非読んで欲しい。とは言え、いきなり読むのは難しいので、RでGLMを使えるくらいにはGLMを理解してから読むことをオススメしたい。

自然科学の統計学

自然科学の統計学 (基礎統計学)

発売日: 1992/08/01
メディア: 単行本

オススメ度　★☆☆☆☆

上で紹介した『統計学入門』の続編。本全体としては必ずしも統計モデルに焦点を当てているわけではないが、第2章で正規線形モデル、第3章で分散分析、第4章で一般化線形モデルが扱われている。

良いところ

第4章ではカウントデータに直線回帰をしてはいけない理由が丁寧に説明されていたり、最尤法がなぜ良いのか（最尤推定量は推定量としてどういった点で優秀なのか）が説明されていたりと、他の本にはあまり書かれていないような細かい話まで抑えている。

良くないところ

前の『統計学入門』よりも数学的な話が多くなり、さらに固く読みづらい（実のところ第2章は難しくてまだ理解できてない･･･）。

総評
推定量の話などは一読の価値があると思うが、万人にオススメできる本ではないと感じる。数理統計学の話題を数理統計学の本を勉強することなしに理解したいという人には良いかもしれない。

R

プログラミングが得意な人はRをわざわざ勉強する必要はないかもしれない（ネット上にも情報は豊富にある）。しかし実際には、R特有の落とし穴や頻出バグがあるので、一度体系的に勉強してみても良いと思う。自分の場合は、「Rでの正しいコーディングルールが知りたかった（スペースの入れ方や慣習的な命名規則）」というのと「グローバル変数の慣習的な扱いが知りたかった（ネットで探すとグローバル変数だらけのコードが多くてちょっと気持ち悪かった）」というのがモチベーションとしてあった。

アート・オブ・Rプログラミング

オススメ度　★★★☆☆

アート・オブ・Rプログラミング

作者:Norman Matloff
発売日: 2012/09/26
メディア: 大型本

いくつかRの本を開いてみたが、この本が一番良さそうだと感じた。信頼と実績のオライリー。

ただし、これはRの文法書であって統計の本ではないことに注意。統計処理を行うRのコード集が欲しいという人には向かない。

2020-10-18

[統計] GLMの最尤推定を自分でやってみる

Python 統計 statsmodels

はじめに

ロジスティック回帰

対数尤度

対数尤度の計算

statsmodelsでGLM

実際の最尤推定の手順

メモ

関連記事

はじめに

GLMで対数尤度を自分で計算して自分で最尤推定をやってみる、という回。

ロジスティック回帰

今回はロジスティック回帰を行ってみる。

ガス濃度とカブトムシの死亡数のデータを使う（Dobson, 2008, p.144）。

x	n	y
1.6907	59	6
1.7242	60	13
1.7552	62	18
1.7842	56	28
1.8113	63	52
1.8369	59	53
1.8610	62	61
1.8839	60	60

統計モデルは次のような式で表される。

$y_i \sim Binomial(n_i, \pi_i)$

$\pi_i \sim N(\beta_1 + \beta_2 x_i, \sigma^2)$

対数尤度

2パラメータのロジスティック回帰の場合の対数尤度関数を示す (Dobson, 2008, p.145)。

$l = \sum_{i=1}^N \left[y_i (\beta_1 + \beta_2 x_i) - n_i \log [1 + exp(\beta_1 + \beta_2 x_i) ] + \log\begin{pmatrix} n_i \\ y_i \end{pmatrix} \right]$

対数尤度の計算

パラメータの値を変えて対数尤度を計算してみる。

from itertools import product

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
from scipy.special import comb


def loglikelihood_function(Y, N, X, beta1, beta2):
    res = 0
    for y, n, x in zip(Y, N, X):
        res += y * (beta1 + beta2 * x) - n * np.log(1 + np.exp(beta1 + beta2 * x)) + np.log(comb(n, y))
    return res


def calc_loglikelihood():
    df = pd.read_csv('beetle.csv')
    Y = df['y'].values
    X = np.ones(shape=[len(df), 2]) * np.nan
    X = df['x'].values
    N = df['n'].values
    RESOLUTION = 100
    beta1_arr = np.linspace(-80, -40, RESOLUTION)
    beta2_arr = np.linspace(15, 55, RESOLUTION)
    res = []
    for beta1, beta2 in tqdm(list(product(beta1_arr, beta2_arr))):
        loglikelihood = loglikelihood_function(Y=Y, N=N, X=X, beta1=beta1, beta2=beta2)
        res.append([beta1, beta2, loglikelihood])
    res = pd.DataFrame(res, columns=['beta1', 'beta2', 'loglikelihood'])
    res.to_csv('loglikelihood.csv')
    print(res[res.index == res['loglikelihood'].idxmax()])
    return


def plot_heatmap():
    df = pd.read_csv('loglikelihood.csv')
    df = df.pivot(index='beta2', columns='beta1', values='loglikelihood')
    df.to_csv('heatmap.csv')
    fig, ax = plt.subplots()
    sns.heatmap(data=df, ax=ax, cmap='YlGnBu')
    ax.set_xticklabels(['{:.2f}'.format(float(x.get_text())) for x in ax.get_xticklabels()])
    ax.set_yticklabels(['{:.2f}'.format(float(x.get_text())) for x in ax.get_yticklabels()])
    fig.tight_layout()
    fig.savefig('heatmap_loglikelohood.jpg', dpi=300)
    return


if __name__ == '__main__':
    calc_loglikelihood()
    plot_heatmap()

         beta1      beta2  loglikelihood
4947 -60.20202  33.989899      -18.72801

f:id:cyanatlas:20201018141256j:plain

図では分かりにくいが、対数尤度が最大になるのは $\beta_1=-60, \beta_2=34$ あたり。

statsmodelsでGLM

statsmodelsライブラリでロジスティック回帰を行ってみる。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf


def glm_statsmodels():
    df = pd.read_csv('beetle.csv')
    df['z'] = df['n'] - df['y']
    glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
    fit = glm.fit()
    print(fit.summary())
    return


if __name__ == '__main__':
    glm_statsmodels()

                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:             ['y', 'z']   No. Observations:                    8
Model:                            GLM   Df Residuals:                        6
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -18.715
Date:                Sun, 18 Oct 2020   Deviance:                       11.232
Time:                        13:56:22   Pearson chi2:                     10.0
No. Iterations:                     6
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    -60.7175      5.181    -11.720      0.000     -70.871     -50.563
x             34.2703      2.912     11.768      0.000      28.563      39.978
==============================================================================

となって、対数尤度は $\beta_1=-60, \beta_2=34$ あたりで最大となることが分かる。この結果は、先に計算した結果と一致する。

実際の最尤推定の手順

ここでは、パラメータに対して総当たり的に対数尤度を計算し、対数尤度が最大になる場所を探して最尤推定量を求めた。しかし、この方法は非効率的である。

そこで、山登り的に計算を進めて最尤推定量を求める方法が考えられる。実際には、ニュートン法（ニュートンラプソン法）をGLM向けに改良したスコア法が用いられ、対数尤度の一階微分（スコア統計量）が0になる場所が探索される。

メモ

このデータと最尤推定の例はDobson (2008)で紹介されているが、p.146の表7.3に示されている対数尤度は桁が1つ間違っていると思われる。具体的には、-186.235となっているところは-18.6235だろう。

cyanatlas.hatenablog.com

2020-10-18

[統計] GLMでの区間推定

Python statsmodels 統計

はじめに

GLM

GLMのパラメータの最尤推定量が従う分布

例：パラメータ2つのロジスティック回帰

信頼区間の検証

正規分布近似とWald統計量

参考文献

関連記事

はじめに

前回の記事では、正規分布の平均を区間推定する方法を説明した。

cyanatlas.hatenablog.com

同様に、一般化線形モデル（GLM）でもパラメータを区間推定することがあるが、この区間推定をどうやっているのか気になったので調べてみた。

GLM

パラメータを $p$ 個持つGLMを考える。

$Y_i \sim \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p$

正確に書けば次のようになる

$g(E[Y_i]) = \boldsymbol{x}^{T} \boldsymbol{\beta}$
$\boldsymbol{x}^{T} = (x_{i1}, x_{i2}, \dots, x_{ip})$
$\boldsymbol{\beta}^{T} = (\beta_1, \beta_2, \dots, \beta_p)$

ここで $g$ はリンク関数、 $\boldsymbol{x}$ は説明変数の行列、 $\boldsymbol{\beta}$ はパラメータのベクトルである。 $i$ 番目の応答変数 $Y_i$ の期待値は、リンク関数を介して、説明変数とパラメータの線形結合によって説明される。

GLMのパラメータの最尤推定量が従う分布

GLMのパラメータ $\boldsymbol{\beta}$ の最尤推定量 $\boldsymbol{b}$ は、近似的に正規分布に従うことが分かっている (Dobson, 2008, p.87)。

$\boldsymbol{b} \sim N(\boldsymbol{\beta}, \mathcal{I}^{-1})$

ここで $\mathcal{I}^{-1}$ は、フィッシャー情報行列 $\mathcal{I}$ の逆行列である。フィッシャー情報行列は、対数尤度から計算される量であり、対数尤度関数が最尤推定量の周辺でどれだけ凸型になっているかを表している（詳しくはフィッシャー情報量 - Wikipedia）。

したがって、最尤推定量は、真の値を平均とし、フィッシャー情報行列の逆行列を分散共分散行列として持つ正規分布に従うということである。上の式を推定量と真の値の差で書き直せば、次のようになる。

$\boldsymbol{b} - \boldsymbol{\beta} \sim N(0, \mathcal{I}^{-1})$

この式から、最尤推定量が真の値とどれだけ離れている可能性があるかを評価することができるので、信頼区間を求めることができる。

例：パラメータ2つのロジスティック回帰

例として次のデータのようなデータ（ $N=40$ ）を考える。

x	y
2.614672698355135	1
0.2311712788458009	0
4.130713626870605	1
3.4980589040126513	0
1.2242786762488267	1
0.5373907435482667	0
1.8970854531195753	1
0.3183052275680498	0
3.9466709572816496	0
3.5309785781654632	1
3.3761561081893414	0
0.720001443910317	0
2.064750554734893	0
2.009122195532624	0
2.2860103892990615	0
4.7927043238872375	1
1.1305070181779313	0
4.027287985243035	0
2.3250037833200765	0
4.554878548309989	0
0.6777507705415364	0
1.0177827317855659	0
4.042219354035978	1
1.5374167164460428	1
0.419902823072304	0
3.2341776629722996	0
1.2807428825075295	0
1.1185882637563382	0
0.37086337116461865	1
0.6935115140150583	1
4.775067677568351	1
3.5673200749785567	1
3.4595515158432804	1
3.195628563216748	0
0.3674178782304305	1
0.4163850367959582	0
0.23154185542532957	0
0.11134517325856519	0
4.287761228825091	1
2.2963424298511153	0

f:id:cyanatlas:20201018014744j:plain

パラメータを2つ持つロジスティック回帰を考える。

$logit(q_i) = \beta_0 + \beta_1 x_i$

$i$ 番目のサンプル $Y_i$ が1となる確率 $q_i$ は、切片 $\beta_0$ とパラメータ $x_i$ の係数 $\beta_1$ によって決まるというモデルである。

ちなみにデータを生成するときに使ったパラメータ（真のパラメータ）は $\beta_0 = -2.0, \beta_1 = 0.8$ である。

ではPythonのstatsmodelsライブラリを使って、ロジスティック回帰を行ってみよう。0の数もカウントする必要があるので、0の数をカウントした列を追加して、ロジスティック回帰を行う。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv('test.csv')
df['z'] = 1 - df['y']
glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
fit = glm.fit()
print(fit.summary())

先に示したデータに対してロジスティック回帰を行い、print(fit.summary())で出力される結果を以下に示す。

                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:             ['y', 'z']   No. Observations:                   40
Model:                            GLM   Df Residuals:                       38
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -24.577
Date:                Sun, 18 Oct 2020   Deviance:                       49.153
Time:                        01:57:24   Pearson chi2:                     40.2
No. Iterations:                     4
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -1.4937      0.648     -2.306      0.021      -2.763      -0.224
x              0.4327      0.231      1.871      0.061      -0.021       0.886
==============================================================================

サマリーの右下に2つのパラメータの信頼区間が表示されている。今回はせっかくなので、この値をフィッシャー情報行列から求めてみよう。

statsmodelsでは、GLMオブジェクトのinformationメソッドを使うことで、フィッシャー情報行列を得ることができる（おそらく符号が逆になっていると思われたので、ソースコードの中では符号を反転させている）。

以下は $\beta_1$ の区間推定を行う例。フィッシャー情報行列の逆行列を求め、その $(1, 1)$ 成分（ $\beta_1 - b_1$ の分散）から標準偏差を求める。±標準偏差×1.96の範囲に95%が含まれるので、ここから信頼区間を得ることができる。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv('GLM_test/test.csv')
df['z'] = 1 - df['y']
glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
fit = glm.fit()
x_est = fit.params['x']  # Maximum likelihood estimation of beta1
fim = - glm.information(params=fit.params)  # Fisher Information Matrix (FIM)
inv_fim = np.linalg.inv(fim)  # Inverse of FIM
std1 = np.sqrt(inv_fim[1, 1])  # Standard error of beta1
lower = x_est - 1.96 * std1
upper = x_est + 1.96 * std1
print('[{:.3f}, {:.3f}]'.format(lower, upper))

出力結果は

[-0.021, 0.886]

となって、先にサマリーで確認した信頼区間と一致する。

信頼区間の検証

GLMのパラメータの信頼区間の求め方が分かったところで、95%信頼区間が本当に95%信頼区間になっているか、データの生成と信頼区間の推定を繰り返して検証してみる。

先程と同じく $N=40$ のデータにロジスティック回帰を行って $\beta_1$ の信頼区間を求める。これを $M=1000$ 回繰り返して、信頼区間の的中率を計算する。

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from tqdm import tqdm

# Parameters of a population
beta0 = -2.0
beta1 = 0.8

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 1000

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def logistic(x):
    return 1 / (1 + np.exp(-x))


def make_data(n):
    x_arr = np.random.uniform(low=0, high=5, size=n)
    err_arr = np.random.normal(loc=0, scale=0.3, size=n)
    probs = logistic(beta0 + beta1 * x_arr + err_arr)
    y_arr = np.array([np.random.binomial(n=1, p=p) for p in probs])
    df = pd.DataFrame({'x': x_arr, 'y': y_arr})
    return df


def calc_confint(df):
    df['z'] = 1 - df['y']
    glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial())
    fit = glm.fit()
    x_est = fit.params['x']  # Maximum likelihood estimation of beta1
    fim = - glm.information(params=fit.params)  # Fisher Information Matrix (FIM)
    inv_fim = np.linalg.inv(fim)  # Inverse of FIM
    std1 = np.sqrt(inv_fim[1, 1])  # Standard error of beta1
    lower = x_est - 1.96 * std1
    upper = x_est + 1.96 * std1
    print(fit.summary())
    return (lower, upper)


def estimation():
    df = []
    for i in tqdm(list(range(M))):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower <= beta1) and (beta1 <= upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(beta1, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を示す。

f:id:cyanatlas:20201018022108j:plain

この例では、的中率は96.5%であり、おおよそ95%になっている。

正規分布近似とWald統計量

以上の方法では、GLMのパラメータの最尤推定量は真の値を平均とする正規分布に従うという法則を用いていた。しかし、これは近似であるということに注意が必要である。データの数が少ないときには近似の精度が悪くなるので、正規分布に近似する代わりに、パラメトリックブートストラップを行う方が良いかもしれない。

最尤推定量と真の値の差が正規分布に従うのであれば、その二乗和はカイ二乗分布に従う（パラメータは独立でないので、分散共分散行列の逆行列を挟む必要がある。詳しくはDobson(2008, p.87)など）。この考え方で求まるのが、Wald統計量であり、このWald統計量を用いてWald検定が行われる。したがって、Wald検定は上で説明した信頼区間と類似した考え方だろうと思うのだが、詳細を知らないので確実なことは言えない。信頼区間の場合は分布の両側から2.5%を削って求めたが、Wald検定の場合はパラメータ $\beta$ が0から離れているかがポイントになるので、0に近い方から分布を5%削って有意性を調べるのだろうか。この辺りは久保(2012, p.51)に簡単な解説がある。とは言え、多くの場合はWald検定ではなく尤度比検定を使うだろうと思うので、この辺りの詳細はあまり気にする必要はないかもしれない。

参考文献

Annette J, Dobson (2008)『一般化線形モデル入門』（訳：田中豊・森川敏彦・山中竹春・冨田誠、原著第二版）共立出版
久保拓弥 (2012)『データ解析のための統計モデリング入門　一般化線形モデル・階層ベイズモデル・MCMC』（シリーズ確率と情報の科学）岩波書店

cyanatlas.hatenablog.com

2020-10-14

[統計] 信頼区間とは何か

統計 Python

はじめに

点推定と区間推定

区間推定の定義

区間推定の妥当性

正規分布に従う母集団の平均を区間推定する（t 統計量を用いる方法）

正規分布に従う母集団の平均を区間推定する（ブートストラップ法）

まとめ

参考資料・参考文献

関連記事

はじめに

与えられたデータの性質を知るために、回帰分析を行って係数（パラメータ）を推定するということが行われる。

このとき、得られた推定結果をp値や信頼区間、AICといったさまざまな指標で評価する。

どの評価指標を使うべきかについては、これまでに数多くの議論がされてきた。ここでその詳細は追わないが、その議論の中で、p値よりも解釈の分かりやすい信頼区間を使うべきだとの指摘もあった。

では信頼区間とは何だったか。どうやって求めるのか。これが今回のテーマである。

点推定と区間推定

母集団のパラメータをどう推定するか。まず思いつくのは、最も可能性の高そうなある1つの値を示す方法である。これは点推定と呼ばれる。

点推定 point estimation：母集団のパラメータをある1つの値で指定する方法。

しかし、推定には常に誤差が伴うものである。点推定の結果が正しいとは限らないので、推定誤差の大きさも評価したい。

そこで、推定誤差の大きさも同時に評価できるような、区間推定という方法が考えられる。

区間推定 interval estimation：真のパラメータの値が入る確率がある値 $1-\alpha$ 以上と保証される区間 $[L, U]$ を示す方法。ここで $1-\alpha$ は信頼係数と呼ばれる。

区間推定に関する上記の説明は、直観的には理解しやすい。しかし、立ち止まってよくよく考えてみると、分からないところがある。

区間推定では「真のパラメータがある区間に入る確率」なるものを考えている。真のパラメータの値は、我々が調査する前からただ1つに決まっているはずなのに、である。

例えば、目の前に密閉されて中の見えない箱が1つあったとしよう。ある調査の結果、「その箱の中に猫が入っている確率は95%以上であると保証されます」と報告された。これは一体どう理解するべきだろうか。箱の中には最初から猫がいるかいないか決まっているはずであり、箱を開ける瞬間に猫がいるかどうか確率的に決まるとは考えにくい（少なくとも我々が通常認識している時空間的スケールでは）。では、同じ箱が100個あったら、95個には猫が入っているということだろうか？しかし、目の前には箱がたった1つあるだけである。

区間推定の定義

区間推定の正確な定義を見よう。『統計学入門』のp.225から引用する。

　同一の母集団から抽出した標本でも、標本ごとに信頼区間の推定値は変化する。 $\theta$ は未知ではあるが決まった定数である。したがって、一つの標本から信頼区間を具体的な数値として推定してやれば、これは信頼区間に含まれるか含まれないかのいずれかしかない。すなわち、具体的に数値として計算した現実の信頼区間に対して、” $1-\alpha$ の確率で $\theta$ を含む”ということはない。信頼区間の意味は、繰り返し多くの異なった標本について信頼区間をここで述べた方法によって何回も計算した場合、 $\theta$ を区間内に含むものの割合が $1-\alpha$ となるというころである。

母集団から取り出された $n$ 個の標本に基づいて、母集団のパラメータ $\theta$ を区間推定するということを、 $m$ 回繰り返すという状況を考えてみよう。

1回目の区間推定の結果を $[L_1, U_1]$ 、 $i$ 回目の区間推定の結果を $[L_i,U_i]$ と書くことにしよう。

すると $m$ 回の区間推定のうち、区間内に真のパラメータ $\theta$ が含まれる回数（ $\theta \in [L_i, U_i]$ である回数）は、 $m$ が十分大きいとき、おおよそ $(1 - \alpha)m$ 回になるということである。

区間推定の妥当性

区間推定の正確な定義を参照した上で、区間推定の妥当性について考えてみる。

ここでは信頼係数 $1-\alpha$ を0.95としておこう。

区間推定では「95%の確率で真のパラメータはだいたいこの区間にあるよ」ということが示される。

では「95%の確率で」というのがどういう意味かというと、区間推定の定義に立ち戻れば、「同じ母集団から今回と同じ数だけデータを取ってきて同じように区間推定した場合、95%の区間推定は真の値を含んでいるよ」ということである。

したがって、今回の区間推定が的中しているのか失敗しているのかは分からないけれども、区間推定が失敗する確率（5%）が無視できるとすれば、真のパラメータは今回の区間推定で得られた区間に入っているはずである。

以上が定義に基づく区間推定の解釈だが、正直言ってなかなかややこしい。

以下では具体的に区間推定を行う方法を説明する。

正規分布に従う母集団の平均を区間推定する（t 統計量を用いる方法）

正規分布 $N(\mu, \sigma^2)$ に従う母集団の平均を、母集団から取り出した $n$ 個の標本 $(X_1, X_2, \dots, X_n)$ から推定する場合を考える。

標本平均 $\bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i$ は $t$ 分布に従うことが知られている。この結果を応用すると、母平均 $\mu$ の信頼係数 $1-\alpha$ の信頼区間は

$\displaystyle [\bar{X} - t_{\alpha/2} (n - 1) s / \sqrt{n}, \bar{X} + t_{\alpha/2} (n - 1) s / \sqrt{n}]$

となる。ここで $s$ は標本分散の平方根であり、

$\displaystyle s^2 = \frac{1}{n-1} \sum_{i=1}^{n}(X_i-\bar{X})^2$

から計算される。また、 $t_{\alpha} (n - 1)$ は自由度 $n-1$ の $t$ 分布の上側確率 $100\alpha$ %のパーセント点である。

では実際に計算してみよう。以下のPythonコードでは、平均 $\mu=2.456$ かつ標準偏差 $\sigma=0.327$ の正規分布に従う集団から、 $N=40$ 個の標本を取り出し、信頼係数 $1-\alpha=0.95$ で母平均の区間推定を行うという操作を $M=100$ 回繰り返している。

import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt


# Parameters of a population
mu = 2.456
sigma = 0.327

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 100

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def make_data(n):
    return pd.DataFrame({'x': np.random.normal(loc=mu, scale=sigma, size=n)})


def calc_confint(df):
    n = len(df)
    x_bar = df['x'].mean()
    t_alpha_half = stats.t.ppf(q=1 - alpha / 2, df=n - 1)
    s2 = np.sum((df['x'].values - x_bar * np.ones(n)) ** 2) / (n - 1)
    s = np.sqrt(s2)
    lower = x_bar - t_alpha_half * s / np.sqrt(n)
    upper = x_bar + t_alpha_half * s / np.sqrt(n)
    return (lower, upper)


def estimation():
    df = []
    for i in range(M):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower < mu) and (mu < upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(mu, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を図にして示す。

f:id:cyanatlas:20201014000003j:plain

垂直な青線で各回の区間推定の結果を表示している。水平な赤線は真の母平均を示す。青線が赤線と交点を持つ場合、その回の区間推定は真のパラメータを含んでいた（的中した）ということになる。逆に、青線が赤線と交点を持たない場合は、その回の区間推定は真のパラメータを含んでいなかった（失敗した）ということになる。

この例では、100回の区間推定のうち、95回では区間推定で得られた区間に真の値が含まれていた。この結果は、信頼係数 $1 - \alpha$ と対応している（乱数によってばらつきがあるので、必ず100回中95回というわけではない）。

ではもっと回数を増やしてみよう。 $M=1000$ 回の例を示す。

f:id:cyanatlas:20201014002510j:plain

この場合は、1000回中967回で区間推定が的中していた。

もっと増やしてみよう。最後は10000回。

f:id:cyanatlas:20201014005904j:plain

10000回中9537回で的中していた。おおよそ95%である。

正規分布に従う母集団の平均を区間推定する（ブートストラップ法）

母集団が正規分布に従う場合、標本平均は $t$ 分布に従うことがあらかじめ分かっていた。上の例では、それを利用して、比較的簡単な計算で区間推定を行うことができた。

しかし、標本から得られた統計量の従う分布が分かっていない場合も多い。また、標本数が多い場合は漸近的に正規分布に従うが、標本数が少ないときは正規分布に近似すると精度が悪いという場合もある。

そこで、統計量の従う分布が数学的に分かっていなくても、コンピュータで乱数を生成することで注目する統計量の分布を生成してしまおう、というのがブートストラップ法である。

ここでは、先程と同じように、正規分布 $N(\mu, \sigma^2)$ に従う母集団の平均を、母集団から取り出した $n$ 個の標本 $(X_1, X_2, \dots, X_n)$ から推定する事例を、ブートストラップ法を使って考えたい（もちろん、本来であれば、この事例には上で説明したような $t$ 分布を使う方法の方が良い）。

以下のPythonコードでは、平均 $\mu=2.456$ かつ標準偏差 $\sigma=0.327$ の正規分布に従う集団から、 $N=40$ 個の標本を取り出し、信頼係数 $1-\alpha=0.95$ で母平均の区間推定を行うという操作を $M=100$ 回繰り返している。

ただし、区間推定の方法は先程とは異なる。得られた $N=40$ 個の標本の平均と（標本）分散を持つ正規分布から $N=40$ 個の乱数を生成してその平均を計算する、という操作を $N_{bootstrap} = 2000$ 回繰り返すことによって、 $N=40$ 個の標本の平均が従う確率分布を計算し（この分布は $t$ 分布と同じになる、はず）、この分布に基づいて区間推定を行っている。

実行には時間がかかるので注意。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm


# Parameters of a population
mu = 2.456
sigma = 0.327

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 100

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def make_data(n):
    return pd.DataFrame({'x': np.random.normal(loc=mu, scale=sigma, size=n)})


N_bootstrap = 2000


def calc_confint(df):
    samples = [np.mean(np.random.normal(loc=df['x'].mean(), scale=df['x'].std(), size=len(df))) for _ in range(N_bootstrap)]
    lower = np.quantile(samples, q=alpha / 2)
    upper = np.quantile(samples, q=1 - alpha / 2)
    return (lower, upper)


def estimation():
    df = []
    for i in tqdm(list(range(M))):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower < mu) and (mu < upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(mu, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を示す。おおよそ95%で的中している。

f:id:cyanatlas:20201014005200j:plain

$M=1000$ 回の場合。こちらもおおよそ95%で的中している。

f:id:cyanatlas:20201014013759j:plain

※ $M=10000$ 回で計算してみたところ、93.99%で的中していた。上の二つの例も踏まえると、95%よりも少ない方にバイアスがかかっているかもしれない。原因が分かったら追記したい。

まとめ

区間推定で得られた区間の厳密な解釈には注意が必要（95%の”確率”で区間は真のパラメータを含んでいる）
区間推定には、あらかじめ統計量の従う確率分布が分かっている場合（ $t$ 分布など）には、その分布を利用して区間推定が行われる。
確率分布が分かっていない場合でも、ブートストラップ法によって区間推定を行うことができる。ただし、計算時間は長くなる。

参考資料・参考文献

東京大学教養学部統計学教室編 (1991)『統計学入門』（基礎統計学Ⅰ）東京大学出版会
Annette J, Dobson (2008)『一般化線形モデル入門』（訳：田中豊・森川敏彦・山中竹春・冨田誠、原著第二版）共立出版
久保拓弥 (2012)『データ解析のための統計モデリング入門　一般化線形モデル・階層ベイズモデル・MCMC』（シリーズ確率と情報の科学）岩波書店
フリーソフトによるデータ解析・マイニング　第44回 https://www1.doshisha.ac.jp/~mjin/R/Chap_44/44.html

cyanatlas.hatenablog.com

グループの分け方	メンバーの分け方
oooo	abcd
ooo \| o	abc \| d
	abd \| c
	acd \| b
	bcd \| a
oo \| oo	ab \| cd
	ac \| bd
	ad \| bc
oo\|o\|o	ab\|c\|d
	ac\|b\|d
	ad\|b\|c
	bc\|a\|d
	bd\|a\|c
	cd\|a\|b
o\|o\|o\|o	a\|b\|c\|d