Pythonによるノンパラメトリック検定

Pythonで並べ替え検定のコードを描いたので晒す。たぶん、scipy.statsとかにあると思うのだが、マニュアルを読むのが面倒くさかった。その記録。コードが欲しい人だけのは、目次のところからコードの項目までどうぞ（ただし、ちゃんと動く保証はないのであしからず）。暇だったので、御託を山ほど書いてしまったが、読む価値はないと思われるので「今日から統計の勉強はじめました（今、中学生）」みたいな人でもなければ、何のためにもならないと思う。

統計的検定

折角の記事だし、入門なので統計的検定について。さて、統計的検定とは何か？それは、あなたの「かせつ」を正当化するために用いられるマジックである。

たとえば、あなたが『幸せになれる！イマイチパッとしないあなたのための、人生哲学セミナー』を主宰したとする。が、いまいち客の入りが悪い。そこで、「私のセミナーに参加すると幸せになれる」という仮説を立てて、それを「でーた」によって客観的に証明することで、集客アップをしようと考えた。

ところで、普通こういう場合には「二重盲検法」というとっても慎重な方法が用いられるのだが、今回の題材は「幸せになれるセミナー」なので、そこまで慎重になる必要もない¹。そこで、一年以上セミナーに参加しているメンバーに「幸せですか？」と訊いてデータを収集するというアドホックな方法を用いることにした。すると、８０％の参加者が幸せですと答えたので、貴方は自信満々で、私のセミナーに参加すると幸せになれますと主張することにした。

ところが、どこからか、ちょっとだけ統計学をかじったやつ（私）がやってきて「コミュニケーションの時間を作っている人の幸福度の平均は7.46点だってさ²。だから、幸せですか？って聞いたとき、八割が幸せだって答えても、差があるとは言えないと思うの」と言い出した。というわけで、貴方はインターネットで「とうけいてきけんてい」という記事（まさにこの記事です！）にたどり着いた…のかもしれない。というわけで、統計的検定とは「差があるかどうか」を調べる「すうがくてきな」方法のことである（すげー、大雑把）。

帰無仮説

では、「差がある」というということをどのように証明すればいいのだろう？

結論から言うと、統計を用いるかぎり「それは、ぜったい、無理」である。

なにしろ、「幸せかどうか？」なんて問題には何の客観的な構造もないのだから、土台無理なのだ。一方で、貴方が工学者で「より効率的な発電機」を発明したという場合は別である。この場合、設計図があれば、物理学が正しいという前提の元で（そして、それはほとんどまず間違いなく正しいので）効率を計算してどの程度の差があるかくらいは証明できる。でも、それがあなたが想定した通り３０％の効率アップなのか、出来たものを動かしてみると２５％程度なのか…それには誤差が付きまとうので、いずれにしろ統計は必要だ。というわけで、問題の構造が割れていても「効率が３０％アップした」かどうかについては、（統計を用いるので）客観的に絶対きっかり証明することは論理的にできない。

では、証明できるのはどんなことかというと、「差がないと言えるかどうか？」である。統計学の教科書に「とは言えない」などという言い回しが頻発するのは、このことによる。で「差がある」ことを証明したい場合はどうすればいいのだろう？その場合、まず次のような仮説を立てる。

帰無仮説：「私のセミナーに参加者と、一般の人々の間の幸せ度には差がない」

というわけで、証明したいことと反対の仮説を立てる。そして、それを否定すればよい。「え？それって、同じことじゃないの？」と思った人は、論理学の修行が足りないので、論理学の復習をする必要がある。実際、機械学習なんかを勉強しようとしても、根本の論理学を知らないと「混同行列」とかのところで、訳が分からなくなるので、下の記事にある図で、今のうちに勉強しておいた方が良いよ（ただ、認知心理学の研究によると人間は「AならばB」と訊くと、うすぼんやりと「じゃあ、BならばAなんでないの？」と思い込む癖があるようなので、こんがらがるのは致し方がないことなのだそうである）。

というわけで、「とうけいてきに」あなたのセミナーの参加者は幸せであるということを証明しようとすると、まず、逆のことを証明しようとして、わざと失敗するという方法を用いなければならない。そしてそれは、論理的に厳密に言うと、「あなたのセミナーの参加者は（非参加者に比べて）幸せである」ということとは全く別のことである。

母集団

ところで、「差がある」というとき、それは何と何との差なのだろう？それを知るために必要なのが、母集団（ここで、エコーがかかる）。

「ぼしゅうだん」とは、いうなれば、この世の神秘、宇宙のどこかにあると言われる「真の」分布である。そして我々の手元にあるデータは、真の分布である母集団からサンプリングされたものに過ぎぬ。たとえば、あなたのセミナーの参加者が30人だとしよう。そして、30人のうち8割が「幸せだ」と答えたのだから、セミナーの参加者のうち80％は幸せだ、といいたくなるがそうではない。なぜならば、母集団からのサンプル（実際の有限の数のデータ）には必ず偏りがあることが知られているからである。一方で、貴方が比較対象として用いる一般の人の幸せ度に関しても、やはり母集団があり、貴方はトランプからカードを引くように、母集団からデータをサンプリング（収集）しているにすぎないのである！

おお、なんと恐ろしい…その母集団の性質をきっかり知るためには、無限のサンプルが必要なのである。もちろん、セミナーの参加者が無限になれば、貴方の収益も無限になってウハウハなのだが、それは現実には起きえないことである。夢を見るのはやめて現実を直視しよう。

というわけで、サンプルを用いて、母集団の性質を推測する必要がある。そして、母集団がある性質に従っている場合は（一例として、たとえば、母集団が正規分布の場合とか！）、サンプルを用いることで母集団のパラメータを近似できることが「しょうめい」されているのだが、それには常に近似による誤差が付きまとう。よく、無限のサンプルを有限のサンプルで近似することによる誤差、などと言われるが、その誤差とはこのことである。

というわけで、「差があるかどうか」ということを証明するために計算する、二つの母集団の差を有限のサンプルで近似したものには、誤差が付きまとう。そのため、誤差を考慮したうえで「統計的に証明」しなければならないと分かった。が、もし誤差が結構大きかったら、どうすればいいのだろう？今回差が出ましたが、誤差の範囲内でしたｗ　みたいなことになると困る。なんとか、ぐうの音も出ないような数字を出してやりたい。

偶然ではない確率、それがp値

そこで、誤差があっても差があるかどうかを、確率的に検証するという戦略を使う。その確率こそ、p値（ここで、エコーがかかる）。

p値とは、あらゆる分野で乱用される便利な数字である。あなたが高校時代、不良グループのパシリにされていた「あいつ」よりよっぽど便利なので、誰もが見境なく「p値が有意でした」と結論することで、自説が正しいと証明できてしまう。まさに、諸刃の剣。

では、このp値とはなにか？それは、ものすごくざっくり言うと、今回サンプルから計算された「差」が、誤差を考慮すると「たまたま」だった、ということがどの程度ありそうか、という確率である。もしたまたま差が出ただけなのであれば、たまたまなので、あなたのセミナーの参加者が幸せなのは、セミナーのおかげ「とは言えない」ということになってしまう（もちろん、セミナーのおかげ「かも」知れず、どちらかは分からない）。

というわけで、貴方の当面の目標は、とっても低いp値を出して、差が出たのは「たまたま」ではありませんと主張することである。差が出たのが「たまたま」ではないのなら、うちの参加者は統計的に「セミナーの非参加者に比べて同程度にしか幸せでない」とは言えないと主張することができる。論理的には厳密でなくとも、先に述べた人間の性質を考えると、このことから直ちに「あなたのセミナーの参加者は、統計的に有意に幸せである」と他人に誤解させることができるはずだ。素晴らしい。

ちなみに、これは貴方のセミナーの参加者の集団としての幸せ度とは異なるので注意。たとえば、貴方のセミナーにはとても不幸な人々の集団ととても幸せな人々の集団、二つの集団の混ぜ合わせなかもしれない。それでも、幸せな集団の方が多数派であれば、無作為に選んだ一般人と比較した場合のp値は有意になるだろう。分布にピークが二つあるようなケースである。単にp値を計算しただででは、そこまでは分からない。そして、サンプルが少ない場合、残念ながら分布の形ははっきりしないことが多いので、たとえp値が有意になったとしても、新しく参加した人が「誰であれ」平均より幸せになる見込みが等しくあるとは必ずしも言えない。つまり、幸せである確率が高いからと言って、それが常に集団としての幸せの質を現しているとは限らないのである。ちなみに平均値は、少数のとても幸せな人がいると歪むこともある。どういうことかというと、たとえ幸せ度を1~10で表現するとして、10人中、幸せ度10の人が1人、5の人が4人、1の人が5人だとする。この時平均は10+5*4+5で35になるので3.5点だが、5点より幸せ度が高い人は1名しかいない（10点満点で、平均が3.5点だからと言って、3人くらいは幸せな人がいるのかな？という計算は常に成り立つわけではないということ）。が、平均値は計算できるし、それが有意であればp値は低くなる。

サンプルの偏り

一方で、先ほどからちょっとずつ問題をごまかしてきたのだが、あなたのセミナーの参加「すると」幸せになれる、ということと、あなたのセミナーの参加者「は」幸せである、ということとは全く別の問題である。統計的な推測ができたからと言って、直ちに因果関係を主張できるとは限らないということだ。これは、質とはまた別の話、より厄介な因果関係の問題である。

もちろん、貴方が量子力学の研究者かなんかで（もちろん、そんな頭のいい人は私の記事は読まないだろうが）、プロセスの因果関係を記述する、信頼のおけるモデルが手元にあるとする。量子力学はたぶん、信頼のおけるモデルである（そのため、量子コンピュータなどというみょうちくりんなものが、現実世界に出現することになった）。そういう場合は、実験結果の統計から胸を張って「量子力学は確率を使うけど、うすぼんやりではありません！ちゃんと因果関係があります！」ということができる。でも、そういう場合は稀である。あとDNAとかなんとか、そういうのも因果関係が一応はある場合に相当するだろう。けれども、それ以外の場合、因果関係を証明することは困難な場合が多い。

たとえば、貴方が今セミナーに在籍している人に、「幸せかどうか」を訊くことで、あなたのセミナーの参加者の母集団としよう、と考えたとする。でも、幸せになれずに辞めた人だって、結構いるはずである。その場合、生存者バイアスという偏りが発生する。つまり、幸せじゃない人がセミナーを辞めた結果、たまたま何らかの原因で幸になれた人だけが、あなたのセミナーに「生存者」として残されるという可能性である。このとき、その人が幸せになった原因は、あなたのセミナーとは何の関係もないかもしれないが、当人からしてみればセミナーに参加したころから、たまたま幸運に恵まれるようになったのである。彼らが、「セミナーのおかげで、良いことが続くなぁ」などとのんきに信じていても責められはしまい。一方、特に幸せにもなれなかったか、特に不幸になった人はさっさと「幸せになれるセミナー」など辞めてしまうであろう。というわけで、あなたのセミナーに現在在籍している人は、無作為に抽出した人々に比べて、幸せであるかもしれない。

もちろん、それは「たまたま」だが、その誤差はいくらサンプルを収集しても、なくならない。なぜなら、そもそも対象としている母集団が違うからである。そして、「幸せになれるセミナー」を不幸になって辞めた人の追跡調査、などという大変な作業をするだけの余裕が、貴方にあるとも思えない。

じゃあ、事前・事後のデータを取って比べればよいのではないか？という気もする。そうすれば、「いんがかんけい」を証明できるのかも。そこで、貴方はセミナーの参加者に事前にアンケートを取り、セミナーを辞退した人にも辞める時にアンケートを取り、生存者にもアンケートを取ったところ、p値が有意だと分かった！これなら完璧でしょう、といいたくなるがそうではない。もしかすると、貴方が統計を取った時期は、景気回復期だったのかもしれない。その場合、あなたのセミナーの参加者が幸福になったのは、景気が回復したせいであり、セミナーのおかげではないと分かる。

（物理法則は、一年や二年で変わるものではないので）量子力学を研究している人はともかく、対象がしょっちゅう変化する場合は、こういうことも考慮しなければならない。というわけで、通常の場合「対象群」という集団の追跡も行う。対象群とは、今回の場合、あなたのセミナーとは縁もゆかりもない人たちである。この人たちを無作為にサンプルすれば、「一般人」の母集団からのサンプルとみなせる。そのうえで、あなたのセミナーの参加者は、一般人（の母集団からのサンプル）に比べて、ある期間に有意に幸せになれました、などと主張することも可能である。

だが、しつこいようだが、本当にそうなのだろうか？もしかすると、あなたのセミナーの参加者は、セミナーに参加したせいで「幸せ」の定義が、比較対象の一般の人々と変わってしまっただけなのでは？と、疑いだすときりがない。そこで、投薬調査などを行う場合（プラシーボ効果を排除するため）二重盲検などで、慎重に比較対象を行う。それは、薬が効くか効かないかは、我々が知り得るかどうかは別として、おそらく客観的な事実だからである。でも、「幸せになれるセミナー」で「幸せ」になったかどうかに、客観的な基準なんてあるのだろうか？

というわけで、今回はあくまでも、あなたのセミナーの参加者は、何が原因かは知らないが、一般人とくらべて、現状、幸せであるかどうかということに的を絞ることにする。いずれにしろ、「母集団」がどんな分布かなど、有限の、しかも偏りがあるかもしれないサンプルしか手に入らない我々には、うかがい知れないことである。

ノンパラメトリック検定

というわけで、なんだか「セミナーに参加すると幸せになれる」などということなど、いくら統計を使っても証明できなさそうな気がしてきたが、この期に及んでもう一つだけ問題がある。それは、母集団の分布が、どんな分布なのか分からないと、上手いことp値を計算できないという問題である。そもそも、たとえ前提がはっきりしたとして、そんなものを計算するには数学をいっぱい勉強しなくちゃならない。それは嫌である。

通常、p値を計算する場合、母集団が正規分布である場合とかを想定することが多い。なぜなら、身長とか、体重とか、大きな外れ値があまりでなさそうなものは、正規分布することが多いからである（身長３メートルくらいならいるけど、１５メートルとかはさすがにあり得ない）。実際、なんかしらの真の値（母集団の平均値）が存在して、その周りにプラスとマイナスの誤差が均等にちらばり、分散が有限で、誤差が指数的に減少するようなもの（こうやって書きだすと随分条件が多いなぁ）は、正規分布になる。そういう場合、平均と分散を手元のサンプルから求めて母集団の推定値とし、偏りを修正して（といっても、Nで割るところをN-1で割るとかその程度）やることで、母集団の平均の差がどのような誤差の分布になるのかを計算できるので便利である。そして、母集団が本当に正規分布なら、それが、一番の早道なのだが、そのためには色々と勉強する必要があるし、母集団が正規分布じゃないケースで困るので、ここではどんな分布でも使える方法を使うことにする。

一方で、我々が推定しようとしている母集団の平均の差の「誤差」の分布というものが存在するのであった。もし、母集団「そのもの」が手元にあれば誤差もクソもなくなり、母集団Aの平均から母集団Bの平均を引いて、誤差ゼロで差を出せるのだが、そのためには無限のサンプルが必要だったんだね。というわけで、母集団の平均（あるいは、何らかの統計値）の差には常に、サンプルによって母集団を推定しなければならないことによる誤差が付きまとう。一方で、この誤差自体は、母集団が正規分布じゃなくても（ずっと緩いある一定の条件を満たすという前提の下で）サンプルをたくさん集めれば正規分布に従うことが知られている。それを中心極限定理という。こういう想定の下で、検定を行うのが、二番目の早道である。が、これでもまだ、誤差の分布として正規分布を使うので数学のお勉強を必要とする。

そして、ここに第三の道がある！とりあえず分布がなんとかとかは、気にしないで行こうや( ´∀｀ )という道である。その方法をノンパラメトリック検定と呼ぶ。一方、上記のように賢い計算をして「この誤差の分布は○○分布（正規分布以外の分布でも良い）になるから～」と何らかの分布を仮定して検定する方法を、パラメトリック検定と言う。

私たちは、小難しい計算などせずに「p値はこのくらいです」とのたまいたいので、ここではノンパラメトリック検定をやってみることにする。

並べ替え検定

そんな風にコンピューターを使ってごり押しでp値を計算する方法、その一つが、並べ替え検定（ここで、エコーがかかる）。

この並べ替え検定とは、言ってみれば、とりあえず計算してみる方式である。母集団の分布が何か気にする必要がなく、小難しい確率分布の計算もいらないという数学音痴の私たちが、大雑把にp値を計算するのにもってこいの方法である。ただし、母集団の分布が分かっているほうが、正確に早く計算できるので、分かっているならその方がいい。でも、分かっているかどうか、どうやってわかるのだろう？私のようなボンクラは、そのあたりが心配でもある。そこで、次善の策に甘んじることにする。

さて、今回は幸福度の差を検定したいのだった。というわけで、一般人の幸福度の平均と、セミナーの参加者の幸福度の平均に差があるかを調査するという風に目標を定める。もちろん、平均じゃなくて中央値でもなんでもよく、それぞれの値に対してp値を計算することはできる。ただ、平均の差の検定というのが最もポピュラーっぽいので、そうしよう。

さて、ここに母集団が二つある。母集団Aは「貴方のセミナーとは縁もゆかりもな人の幸福度の（真の）分布」である。母集団Bは「貴方のセミナーに参加した人たちの幸福度の（真の）分布である」この二つの分布が、どんな形で、パラメータがどうなっているかは分からない。それは、絶対に分からないのだが、少なくともそこから抽出されたと思しきサンプルがある。

さて、このサンプルを使って、出来る限りのことしよう。まず、「差がある」と言いたいので「差がない」と仮定せねばならぬのだった。そこで、差がないと仮定する。そして、真の母集団については知りえないので、とりあえず、手元のサンプルを母集団の代表値と考えることにしよう。このとき、AとBに差がないのだから、この二つからどんな組み合わせでデータを取り出しても、平均の差は出ないだろうということができる。というわけで、本来の並べ替え検定では、実際にあらゆる組み合わせを計算したのち、組み合わせの総数で割るという処理をするのだが、データ数が500とかになったらどうするのだろう？その場合、順列の数はとてつもないことになりそうである。

そこで、モンテカルロ近似を使うことにする。モンテカルロ近似とは、とてもじゃないが計算できない分量の計算をしなければならないとき（たとえば、めちゃくちゃ複雑な連続関数の積分とか）統計的なサンプリングを行うことで、全体の値を近似的に求めるという方法である。ただし、ちゃんと使えば、サンプリングの回数をいくらでも増やせるならば、近似による誤差をいくらでも小さくできるという性質を持つ。今回は、並べ替えが凄い数になりそうなので、そこをモンテカルロ近似で計算することにしよう（が、これがちゃんとした方法かどうかの確信はちょっとないかも…識者の方教えてくれるとありがたす。たぶん大丈夫だと思うのだが）。

というわけで、まず、Aからのサンプルと、Bからのサンプルを「まぜこぜ」にしたサンプルを作る。もし、AとBの平均に差がないのなら、ここからどんな組み合わせで取り出して平均を計算しても、差はないはずなのであった。というわけで、そういう計算をたくさんしてみよう。すると、色々な組み合わせ（ほんとうなら、「あらゆる」組み合わせについて計算するのだが、そこはモンテカルロ近似にするのでした）に対して、平均の差を計算することができる。

すると、この色々な組み合わせに対する平均の差、の分布が得られる。これは、たぶん正規分布なのだろうが、そうでなくてもよい（いや、はず。自信なし）。というわけで、これは先ほどから考えていた、母集団をサンプルで近似することによって生じる、母集団の統計量（平均とかメディアンとか）の差の誤差の分布ということになる。

これで、分布が手に入ったのだからこっちのものである。あとは、大元の平均の差（Aから抽出されたサンプルから普通に計算した平均と、Bから抽出されたサンプルから普通に計算した平均の差）が、この誤差の分布と比べてどのあたりにあるか？が計算出来ればよい。

正確な「並べ替え検定」とちがって、モンテカルロ近似してしまったので、「実は近似による誤差」がさらに上乗せされているのだが、一万回くらい計算すればそれなりの値になるし、それでもデータが数百の時、全部の順列に対して計算するよりはマシである。というわけで、そのようにコードを書いてみよう。

コード

コードは下記の通り

def permutation_test(func, y, z, B, smooth=False):
    n, m = len(z), len(y)
    x = np.concatenate([y,z])
    xs = np.repeat(np.array([x]), B, axis=0)
    # use permutation instead of bootstrapping
    ssamp = np.apply_along_axis(np.random.permutation, 1, xs)    
    ztheta = np.apply_along_axis(func, 1, ssamp[:,n:])
    ytheta = np.apply_along_axis(func, 1, ssamp[:,:n])
    theta = ztheta - ytheta
    diff = func(z) - func(y)
    return np.sum(theta > diff) / B

ちなみに、ここでnp.random.permutationを使っているのは、「並べ替え」だからである。ここを「復元抽出」にするとブートストラップ検定っぽくなる。復元抽出と、復元しない抽出の間にどのような差があるのか？理論的には知らないのだが、復元抽出したほうがちょっと緩めのp値になることが多い。というわけで、p値とかいっても、計算の方法によって出てくる値が結構ちがーう。まぁ、我々は「かがくしゃ」でもなんでもないので、こんなもんでいいんじゃないかな？

なぜこんなうんこ記事を書くのか？

多分10年くらい前、「データを出して議論しよう」という、そういう文章をネットで読んだ記憶がある。そして、その時はその通りだと思った。ところが、最近、あちこちをみると「とりあえずデータを出せば」よいみたいな感じになりつつある気がする。ぶっちゃけ、とりあえずのデータなど何の根拠にもならないのではないかと思う。注釈のところで文句を言った幸福度指数的なものもそうだ。考えてみて欲しいのだが、北朝鮮に住んでたって、政府の要職で平壌の一等地にすんでいるのなら、私のようにクソ田舎の最底辺の日本人でいるより幸せに決まっている。あるいは、つい先日まで飢えていた人たちが、「なんとか食べれる」ようになったのなら、私以下の生活でも十分幸せに感じることだろう。

ようするに、「なんとか財団」みたいな組織に山のように金をつぎ込んだところで、世の中クソもよくならないのではないか、ということである。実際、SDGsとか何とかいったところで、「とりあえずのデータを出して」たいして益にもならないことに金をつぎ込む方式のビジネスになり果てるのが関の山なのではないだろうか…。

いや、なんか、最近後ろ向きなんです。なんかすべてが疑わしいもんねー。だって、温暖化でやがて人類は滅びるってときに、夜中まで煌々と電気つける？いや、電力需給をある程度維持しないと困るってのは分かるけどさ、節電の時は消すじゃん。全部消せとは言わんけど、減らすことくらいできるのに、全然減らないとこ見ると、いい加減疑わしいと思ってしまうよ（温暖化ではなくて、それを解決しようと努力しています、みたいな話の方ね）。

以上。

1. 本当のところ、こんな「胡散臭い」例を持ち出したのは、私が間違ったことを言っても、「でも、これは胡散臭い例だからｗ」と言ってごまかすことができるためである。 ↩

2. https://prtimes.jp/main/html/rd/p/000000802.000000983.html　
   こんな調査があるんですね、知らんかった。ただ「幸せですか？」と訊いて幸せだと答えたからと言って、幸せかというとどうなんだろう？と思うのは私だけでしょうか。基準だって社会状況に応じてコロコロ変わりそうなので、そもそも「客観的な幸せ度の分布」なんてないんじゃないのかなぁ。一方で、「世界幸福度ランキング」なんてのもあってよくテレビで見かけますがアメリカが16位、メキシコ46位、日本が54位だそうです。いやね、日本は決して高い方じゃないと思うよ、幸福度。でも、もし米国やメキシコの方が幸福なら、なぜメキシコの人たちがアメリカに逃げ出す傍ら、海外に逃げ出す日本人は少ないのでしょう？まあ、27位のシンガポールあたりは、日本から「逃げ出した」富裕層が結構移住しているらしいので、日本よりずっと住みやすいのでしょうね。色々うるさいし、お高そうですが。ただ、台湾は26位ですから、物価の差を考えてると、そんなに幸福度が低いなら、日本から台湾に逃げ出す人がもっといてもいいのではないかと思います。というわけで、「統計を取ったから」といって常に信頼できるか？といわれると、私はそんなことは全然ないと思っています（治験とか、製品検査は別ね。ああいうのは、それなりに信頼できるのだろうけど）。ま、この先は分からないですけどねｗ。20年後には日本から中国に移民してたりしてね。 ↩