Easy to type

個人的な勉強の記録です。データ分析、可視化などをメイントピックとしています。

ScipyでBrunner-Munzel検定

概要

Brunner-Munzel検定は、不当分散のときに使えるノンパラメトリック検定です。ランクから、得られたデータの中央値に有意差があるか検出します。 ノンパラメトリックの検定では、Wilcoxon-Mann-WhitneyのU検定が非常に有名ですが、この検定は不当分散の場合には第一種の過誤を得る確率が5%にならない場合があり、適切な手法として使えません。 詳しい調査には、以下のような先行記事があります。

hoxo-m.hatenablog.com

Brunner-Munzel検定

二群の平均値(代表値)の差を検定するとき

しかしBrunner-Munzel検定は、有用性のわりにどマイナーな検定なので、Rではlawstatというパッケージで利用が出来ましたが、Pythonにはパブリックで使えるいい感じのライブラリがありませんでした。 今回アップデートされたScipy Version 1.2.0でBrunner Munzel検定が使えるようになりました。使い方の記録です。

使い方

ドキュメントはこちらです。上述した奥村先生の記事と全く同じやり方にしてみます。

from scipy.stats import brunnermunzel

a = [1,2,1,1,1,1,1,1,1,1,2,4,1,1]
b = [3,3,4,3,1,2,3,1,1,5,4]
brunnermunzel(a, b)

結果

BrunnerMunzelResult(statistic=3.1374674823029505, pvalue=0.005786208666151538)

統計値やP-valueも一致しました。デフォルトでは、統計値をt分布に当てはめて推定しています。 原著論文中ではサンプルサイズが50未満であるときにはこちらを利用し、それ以上だと正規分布に対して当てはめています。正規分布に対して使いたい場合はdistributionのオプションをnormalに変更すれば使えます。 検定を両側、片側でやる場合にはalternativeが使えます。

皆様も、ぜひBM検定を使って「まだ不等分散で消耗してんの??」とドヤ顔してください。

本文

さて、この記事の目的はBM検定の宣伝でもあるのですが、真の目的はScipyにCommitする感想です。このBM検定、及びScipyのAPIドキュメント*1は僕が書きました。当該PRはこちらです。 自身のメイン言語がPythonだったこともあり、BM検定は有用そうでしたから、車輪の再開発で消耗するような人を減らしてみたかったという想いで取り組みました。 僕はOSSにコミットするのは、ほぼ初めてだったのでとてもいい勉強になりました。以下箇条書きです。

  • コード規約について
    • 基本的に、WMW検定の書き方へ似せた。ドキュメントも同じ。似た構造の検定なので、悪い戦略ではなかっただろう。
    • Scipyではgitのコミットメッセージにも規約がある。自身はそれまで非常にテキトーなコミットメッセージを書いていたので、一つの目安として自分のその他のプロジェクトでも使うようになった。
    • 自動シンタックスチェッカーを入れてFlake8で静的テストするようになった。
  • コードレビュー
    • しょーもない変数の使い方などについてもコメントをくれた。
    • 科学計算をする上で桁落ちや型変換によるエラーは大きな問題だが、それについて指摘してくれ、気がつくことが出来た。以降わりと気を使っている。
  • マージまで
    • 非常に時間がかかった。
    • 最初にコミットしたのが2017年4月だったわけだが、放置され気味でmasterにマージされるのに1年放置。
    • 何故か当初予定していたMilestone(1.1.0)から一つ先送りにされて、1.2.0でようやくマージ。結局Stagingされるまでは1年と8ヶ月かかった。
    • 1Milestoneが約半年で締められ300commit進むこと、オーガナイザーの苦労を考えると仕方ないが、もうちょっとなんとかならないのかなと思わなくはない(ガンガン貢献してオーガナイザーとして働け、というデカイまさかりが飛んできそう)。

また機会があれば、Scipyに限らずOSSに取り組んでみたいです。

*1:今読むと文法ミスが多くてしんどい。