続・自由度について（実験編）

この記事は不正確な内容を含みます。下記の新しい方の記事が正確です。（2017/9/18）
nekomosyakushimo.hatenablog.com

ーーーーー

前の記事で、標本から母分散を推定するときに、nではなく自由度であるn-1を用いることについて書いた。今回はそのことをRを用いた簡単な実験で示す。

母平均が100、母分散が225（母標準偏差は15）の正規分布からn=10の標本を抽出して、その分散を求める。この手続きを10000回繰り返して、経験的に標準偏差の標本分布を作ることにする。これを、 n で割って普通に分散を求めた場合と、n-1で割って不偏分散を求めた場合とを比較するのが今回の目的である。

スクリプトは難しいところは何もないが、Rのvar( )関数は不偏分散を求める仕様であるので、普通の分散を求める方では、平均と標本サイズから定義通りの計算をしている。

#標本分散と不偏分散の実験
sample_var <- numeric(10000)
unbias_var <- numeric(10000)

#普通の分散
for (i in 1:10000){
  dat <- rnorm(10,100,15)
  dat_mean <- mean(dat)
  sample_var[i] <- sum((dat - dat_mean)^2)/10
}

#不偏分散
for (j in 1:10000){
  dat <- rnorm(10,100,15)
  unbias_var[j] <- var(dat)
}

得られた結果をヒストグラムにして重ねたのが次のもの*1。

ピンクの方が nで割った普通の分散で、青い方が n-1 で割った不偏分散である。当然のことながら不偏分散の方が平均的に大きくな値を取っている。

それぞれの平均（つまりこの経験的な標本分布における期待値）を求めると以下の通りである。

> mean(sample_var)
[1] 200.6958
> mean(unbias_var)
[1] 223.8554

不偏分散の方は、母分散である225に近い値だが、普通の分散の方はやはり小さく推定してしまっている。普通の分散の方に、n / n-1 をかけると不偏分散による推定量に近づく。

> mean(sample_var) * 10/9
[1] 222.9954

以上、簡易的な実験により経験的に求めた分散の標本分布からも、不偏分散が不偏推定量になっていることが確認できる。

*1:色付きヒストグラムを重ねるときには ヒストグラム | を参照した。

自由度について

統計関係の概念の中でも理解しづらいものの一つに「自由度」がある。得られたデータから、母分散の不偏推定量を求める時に、nではなくて、n -1　で割るだとか説明されて分かったような分かっていないような気になるアレである。それに関連する学習のメモ。

まず、標本平均の標本分布の話から。ある母集団から標本サイズがnである標本を抽出した際に、その標本平均は様々な値をとる。例えば、母平均が100、母分散が225（母集団標準偏差は15）のIQをとる集団があったとして、その集団から適当に10人ぐらいを抽出して標本平均を求めると、当然のことながら母平均と同じ値になることはあまりない（例えば、102かもしれないし90かもしれない）。ただ、母平均から離れた値はとりづらい（例えば標本平均が60なんかになることは稀である）。このようにして、「標本平均がどのような値を取りやすいか」についての分布を標本平均の標本分布と言う。

この標本平均の標本分布は、平均が母集団と同じμ、分散がσ/nを持つ正規分布に従うことが知られている。

IQを例で言えば、この標本平均のとる値は、平均100、分散22.5の正規分布に従い、確率的に変動するということである。

ここまでを前提に自由度の話に入る。標本平均の期待値は、母集団の平均と一致している。こういうものを不偏推定量と言う。では、標本の分散の期待値はどうだろうか。結論から言えば、母分散よりも小さい値になることが分かっている。つまり、標本の分散を母分散の推定量として扱うと、母分散を小さく見積もることとなってしまう。

なぜそのようなことが起きるかというと、標本の分散を算出する際に標本平均を使っていることに由来している。標本平均は先に確認したようにそれ自体が標本抽出のたびに変わる確率的なものであり、μを中心に、σ/nの分散で分布している。その関係もあって、母分散は、標本平均を中心とする分布の分散（標本の分散）と標本平均の標本分布における分散の和と等しくなる。このことを表したのが次式である。

右辺の第１項が標本の分散の部分、第２項が標本平均の標本分布の分散である。この式を次のように変形していく。まず、第２項を左辺に移行して通分する。

左辺のσをくくりだす。

両辺をnで割って（n -1 ）を右辺へ移行して完成。

というわけで、件の n -1 たる数字が導出された訳である。最後の式の右辺で得られた推定量は不偏推定量であり、不偏分散と呼ばれることもある。n - 1というのは、標本の分散から母分散を推定するにあたって、標本平均それ自体のバラツキを勘案し、過小推定を調整したものだと言える。

【参考】
「心理学のためのデータ解析テクニカルブック」（北大路書房）

連続のための修正とは何か

連続のための修正（連続性の補正）について調べたことのメモ。（まだ理解は不完全なのであまり信用できる記事ではありません。）

やや遠回りではあるが、まずは離散変量での検定について考える。

コインを10回投げて表の出る回数（確率変数）というのは、n=10, p=0.5 の二項分布に従う。それぞれの変数のとる確率はdbinom関数で簡単に求められる。以下はその表および確率分布の図である。（z値はあとで使うので今は無視してよい）

表の回数	確率	z値
0	0.001	-3.162
1	0.010	-2.530
2	0.044	-1.897
3	0.117	-1.265
4	0.205	-0.632
5	0.246	0.000
6	0.205	0.632
7	0.117	1.265
8	0.044	1.897
9	0.010	2.530
10	0.001	3.162

さて、ここでこのコインが本当に0.5の確率で表がでるものなのか（細工がされていないか）調べたいとする。帰無仮説をp=0.5として、例えば、表の回数が（0、1）回あるいは（9、10）回だったらこの帰無仮説を棄却するものとしよう（極端に裏ばかりでる、あるいは極端に表ばかりでるようだったらp=0.5ではないとする両側検定）。通常だと有意水準を定めて、その水準を超えるようだったら帰無仮説を棄却とするだろうが、とりあえずここでは、上記の回数だったら棄却すると設定する。このとき、第１種の誤りαの確率は、それぞれの確率変数のとる確率の和である。すなわち、

coin_prob <- dbinom(0:10,size=10,prob=0.5)
sum(coin_prob[1:2]) +sum(coin_prob[10:11])
[1] 0.02148438

であり、0.021%程度の確率で帰無仮説が真なのに、それを誤って棄却してしまうことになる。

ところで、二項分布はnが大きくなれば正規分布に近づくことが知られている。この性質を利用して、得られた値について検定を行おうという考えがでてくる。実現値を x 、平均を μ 、標準偏差を σ とすると、次のように z へと標準化ができる。

この変換をして得られたzの値は、最初に載せた表に記載してある。この z は近似的に、平均0、分散1の標準正規分布に従う。

さて、このzの値を用いて有意確率を計算してみる。表の回数が（0、1）回あるいは（9、10）回というのは, z値が |z| > 2.53 になる確率なので、正規分布が平均に対称なのを利用し、

pnorm(-2.53) * 2
[1] 0.01140625

と求められる。αの確率は、0.011%となり、最初に二項分布を用いて求めたときよりも、第１種の過誤を犯す確率がかなり小さく評価されてしまうこととなった。

これは、もともと二項分布が離散変量に基づく分布にもかかわらず、xを連続変量として扱っていることに原因があるとされている。離散的な整数を連続変量として当てはめているわけだが、それらの離散的な整数と整数の間にも確率があるわけで、この確率を考慮しないと、連続的変量としての分布への当てはまりが悪くなる。

※下の方にその解説がある。
統計学入門−第3章

離散変量をもとに連続変量の分布に近似して検定を行う際には、次に示す連続のための修正（イェーツの修正や半整数補正とも呼ばれる）を施した方が正規分布への近似がよくなるとされる。

今の例で考えると、1.5、8.5のところからz値を求めると正規分布により正確に近似されるはずである。上の式をもとにして z を修正し、その値を用いてαの確率を求めると以下のようになる。

(abs(1 - 10*0.5) - 0.5) / sqrt(10*0.5*0.5)
[1] 2.213594
(1 - pnorm(2.213594)) *2
[1] 0.02685672

αの確率は0.027%で、何も考えずにz値に変換した場合に比べて、二項分布で求めていたときの有意水準αに近づいており、第１種の過誤の過小推定の危険性を避けている。（今度は大きく評価するようになったが。）

とりあえずここまで。

【参考】
『心理学のためのデータ解析テクニカルブック』（北大路書房）

心理学のためのデータ解析テクニカルブック

柴田淳『みんなのPython』

みんなのPython 第4版

『心理学実験プログラミング』（朝倉書店）からPythonに入り、一応簡単な実験プログラムなら書けるようになったものの、仕組みがよく分からぬまま動かしているところもあったので、Pythonという言語そのものについて知ろうと思い読んだ。

全部で13章からなり、前半はPythonをインストールするところから始まり、基本的な関数、組み込み型の使い方へと続く。中盤は、オブジェクト志向の概念を中心に、クラスやモジュール、標準ライブラリの扱い方を説明している。最後には、データサイエンスでの活用例としてNumPyやmatplotlibなどのライブラリでの実践例を紹介している。

自分は、Rしか使ったことがなかったので、メソッドやらクラスやらアトリビュートやら、オブジェクト志向の言語の仕様が分かっていなかったのだが、この本の真ん中のあたりの章（5~9ぐらいまで）を読んでだいぶ見通しがよくなった。

よくよく読むと『心理学実験プログラミング』にも同様のことが書いてある（った）のだが、『みんなのPython』で１章割かれているような内容が、1段落で済まされていたりするため、Python初心者の私のような人間はさらっと未理解のまま流してしまっていたりする。

Pythonそのものへの入門としてはとてもよく出来ているように思ったので、これからPythonをやってみようという人にはオススメできると思う。ただ、これ1冊読んでも、作りたいものまでの距離感というのは分からないような気もするので、作りたいものが明確な人は、それに特化した入門書なりに手をつけて（写経し）、仕組みはそこまで分からないけど動かせるようになってから、その仕組みを詳しく学ぶために本書のような入門書に手をつけても良いのかもしれない。モチベーション的に。

次は、画像処理系あたりに手をつけていこうと思う。

pep8とflake8の導入について

pep8というものがある。自閉症療育に関心がある方なら、TEACCHプログラムの開発したPEP-3（自閉症・発達障害児教育診断検査第３版）というのが思い浮かぶかもしれないが全くもって関係ない。

pep8とは、Pythonでコードを書くときのスタイルガイドである。これに準拠しなかったらプログラムが動かなくなる訳ではないが、一貫したスタイルのコードを書くことで可読性をあげようとかそんな話だと思う。

pep8の日本語訳された文章は以下から読める。

はじめに — pep8-ja 1.0 ドキュメント

自分が書いたコードがこのスタイルに準拠しているか調べるための文法チェッカーというものもある。便利な世の中である。

いくつかあるうちでflake8というものが有名らしいので試しに使ってみる。

導入はpipが入っていれば簡単で、ターミナル上で以下のコマンドでインストールされる。

$ pip install flake8

使い方もまた簡単で、以下のように引数にファイル名を指定してやれば良い。

$ flake8 test.py

例えば、次のようなコードをチェックするとする。

a = [0,1,3  ]

print("Hello World!") #ハロワ

するとこんな感じで、スタイルに沿っていないところを指摘してくれる。

test.py:1:7: E231 missing whitespace after ','
test.py:1:9: E231 missing whitespace after ','
test.py:1:12: E202 whitespace before ']'
test.py:3:22: E261 at least two spaces before inline comment
test.py:3:23: E262 inline comment should start with '# '

何行目・何文字目に、どんな違反なのかが書いてある。エラーコードの頭にあるアルファベットは何についての違反かを示しており、Eはpep8のスタイルに準拠していないことを示している。

例えば、1行目の7文字目のところでは、「カンマの後に空白が入ってないよ」と教えてくれている訳だ。

Error / Violation Codes — flake8 3.4.1 documentation
Introduction — pycodestyle 2.3.1 documentation

このエラーメッセージを元にさっきのコードを修正すると次のようになる。

a = [0, 1, 3]

print("Hello World!")  # ハロワ

これで、pep8に準拠したリーダブルなコードになった訳だ。

このflake8はAtom上で動くパッケージも出ており、linter-flake8という名前である。

GitHub - AtomLinter/linter-flake8: Linting Python files on the fly using flake8 with Atom

特に難しいこともなく、上記のREADMEに従ってinstallするとAtomでコードを書きながらリアルタイムで指摘してくれるようになった。

どうせ書くなら読みやすく美しいコードを、ということでflake8の導入の話でした。

信号検出理論事始め

信号検出理論というものがある。もともとは、レーダーの性能など通信工学的な分野で発展した理論でそれが感覚や知覚の心理学研究に応用されるようになったらしい。

以下の本を参考に、自分の学習のメモを残す。（記事の内容の正しさは保証できないため、しっかりと学びたい人は原典をあたることをおすすめする。）

Detection Theory: A User's Guide

はじめに ~考え方編~

まず信号検出理論の特色を一言で言ってしまえば、それは、観察者の反応を（1）提示される刺激による要因と（2）反応のバイアスによる観察者側の要因の2つに分離して記述できることだと思う。古典的な精神物理学の研究では、刺激側の要因と心理反応の関数関係の記述を主としていたが、反応における認知的なプロセスの部分も考慮に入れてモデル化しているというのは、この理論の強みであるように思う。

これだけだと「何のことだ？」と思う方もいるかもしれないので例を出して考えよう。

例えば、次のような事態を考えてみる。ここに嘘発見器があるとする。発見器をつけた状態で話をすると、発見器は「嘘」か「正直」のどちらかの回答を返すものとする。さて、この発見器の性能を検証しようと考え、次のような実験を行う。発見器をつけた状態で、それぞれ25回の「嘘の話」と「本当の話」をする（こういうのを刺激クラスと呼ぶ）。「嘘の話」に正しく「嘘」と言うことができ、本当の話に対して「正直」ということができる割合が高ければ、この発見器は優秀だと言うことができる。この架空の実験の結果を表で示す。

刺激クラス	「嘘」と反応	「正直」と反応	合計
嘘の話	20	5	25
本当の話	15	10	25

さて、この嘘発見器は優秀なのであろうか。嘘の話に注目すると、25回中20回も嘘を見破っている。確率で言えば80%の確率で嘘を発見しているので優秀だと言えるかもしれない。しかし、本当の話に注目すると、25回中15回も間違って嘘だという反応を返している。確率でいえば60%の確率で間違った判断を下している。これでは、この発見器が優秀だと言いきることはできなかろう。

このように、片方の刺激クラスにのみに注目していたのでは間違った推論をしてしまう恐れがある。したがって、嘘発見器の正確性を評価しようと思ったのであれば、どちらの刺激クラスも考えなければならない。

ところで、嘘発見器の性能を正確性という一つのパラメータで記述してよいものであろうか。ここに嘘発見器Bがあるとして（さっきの嘘発見器をAとする）、先ほどと同様の実験を行った結果が下記の通りだとする。

刺激クラス	「嘘」と反応	「正直」と反応	合計
嘘の話	16	9	25
本当の話	11	14	25

嘘の話に注目すると、64%の正確性で嘘を見破っており、嘘発見器Aよりも嘘を見破る力は弱い。本当の話に注目すると、44%の確立で間違った判断をしてしまっていて、嘘発見器Aより間違って「嘘」だと反応してしまう率は低い。Aの発見器よりも「保守的な」発見器だと言えるかもしれない。

どちらの発見器も、正しく判断できた数（左上と右下のセルの合計）は30であり、同程度の正確性で嘘を検出している。しかし、この２つは嘘検出の正確性では同程度だとしても、その反応のパターンは結構ちがっている。こうした、反応パターンの違いを考慮してモデルを構築するのが信号検出理論の大きな強みだろう。

信号検出理論の用語たち

さて、ここで信号検出理論の用語を導入する。一番最初の嘘発見を例にして話を進める。2つの刺激クラスについて2つの反応パターンがあった訳なので、この組み合わせは、2×2で４つのパターンがあることになる。その用語を以下の表に示す。いかにも、もともとレーダーの研究で使われていたっぽい用語である。

刺激クラス	嘘	正直	合計
嘘の話	ヒット（Hits） （20）	ミス（Misses） （5）	25
本当の話	誤警報（False Alarms） （15）	正棄却（Correct Rejection） （10）	25

この表には、4つの単語があるが実質独立しているのは2つだけなので（例えば、ヒットの数が決まればミスの数は必然的に（合計 ー ヒット）になる）、この表の特徴は2つの値で表すことができる。すなわち、ヒット率をH、誤警報率をFと表すと次のように書くことができる。（F , H）=（.60, .80）

まず、求めるのは検出力を表す指標である。検出力はHが高い時に高くなり、Fが高い時に低くなる（なにせ間違った判断を下している訳だから）ものが望ましい。この条件を満たすシンプルなものは、H - Fである。今の発見器の例だと、 .80 - .60 = .20 になる。

あるいは正しい判断を行なった確率というのでも良い指標であろう。ヒット（H）と正棄却（1 - F）を2で割ると、試行数における正しい判断の割合が算出できる。

今の例だと、( .80 + .40) / 2　= .60　になる。

ここまでで2つの指標を導入した訳だが、信号検出理論でもっとも広く使われる指標はもう少し複雑である。 HとFの確率を、正規分布における確率点に変換した値を計算して得られるd'（読みはディープライム）というものが使われる。

先の例で言うならば、z(H) = 0.842　、 z(F) = 0.253であるから、d' = 0.589 となる。Rで確率点を求めたいのであれば,qnorm()という関数を使うと良い。

> H <- qnorm(0.8)
> F <- qnorm(0.6)
> H - F
[1] 0.5882741

とりあえず、これで検出における正確性を表す指標が手に入ったこととなる。今日はここまでで、次回以降で正確性を表す指標とは違う、反応のバイアスに関わる指標について書いていこうと思う。

テキストエディタとモチベーション

モチベーションっていうのはなかなか制御するのが難しいものであり、かなりのところ環境要因の影響を受けるものである。スキナー大先生のお言葉を借りれば、私的出来事（private event）であり、私的出来事も他者から観察不能なだけで、三項随伴性の支配下にあるのだろう（テキトー）。

何が言いたいかというとテキストエディタを変えるだけで、モチベーションはぐぐぐっと上がるということなのだ。

【Before】

【After】

Atom、かっこいいよ、Atom。