Erlang: 組合せをやってみる

さて、少しやる気も出てきたので、 Erlang の勉強編です。

私は大学生の頃、(実験系の)心理学科だった...というのもありまして、ゼミの内容はマジに「統計学」みたいなものでした。この中で取り上げられていたのが、「パーミュテーション検定（並べ替え検定）」という技法です。世の中にはいろいろな検定手法があるものですが、要するにこの「パーミュテーション検定」というのは、

帰無仮説＝実験結果で実験群と対象群とでは、その平均値の違いは偶然の産物に過ぎない

というケースで、「ということは、実験群と対象群の差は見かけだけなんだから、その区別をなしにして、全部バラバラに組み合わせた時の平均値の違いを見る、と考えてみれば、その時に実験で得られた値を越えるような組み合わせは、いくらでもあるはず...」というかたちで考えるものです。だから、もしそれを超える平均値の出現頻度が５％以下ならば、「95％有意である」と結論できることになります。

実際、こういうケースでは普通、分散分析で分散の値を比較して検定するのですが、この「パーミュテーション検定」だと、分散分析の前提である「正規分布をしていること」を全然考慮しなくていいわけです。また、これだと、たとえば実験結果が「○×」でしか得られないような実験であっても、「○×の個数の偏り」が有意かどうかさえ、検定できてしまう（いわゆる「ノンパラメトリック検定」というものです）...というような、「検定」というものの根本から考えたような手法なのです。

ですから、これ「すべての組合せを求めて、それぞれに対して独立に判定を行い、その総計が全体に対して何％になるか」を計算する、という問題になるわけです....特に Erlang っぽいネタですよね！！！

でしかも、「組合せを生成する」というのはちょっと書いてみる価値のある、アルゴリズムの面倒なプログラムでもあります。実は Erlang の場合、順列を求めるだけなら、こんなコードで書けちゃいます。

....とはいえ、これはちょい反則です。また、たとえばこの順列の結果から「前からk個取って、重複は無視／削除する」というアルゴリズムがないわけではないですが、順列の個数 >> 組合せの個数 なので、あまりに不効率です。まあ、ちょいとここから書いてみましょう。

-module(comb).
-export([comb/2,stripcomb/1]).
-import(lists,[map/2]).

%% 外部公開の呼び出しインターフェイス
%% comb( リスト, リストから「取る」数 )
comb([], _) -> [];
comb(_, 0) -> [];
comb(X,K) when is_integer(K), length(X) < K -> [];
%% 実際の計算がある部分は下請けに回す
comb(X,K) -> comb2( initcomb(X), K).

% 準備：[1,2,3] -> [{1,[2,3]},{2,[3]},{3,[]}] のようなタプルに変換する
initcomb([]) -> [];
initcomb( [H|L] ) -> [{[H],L} | initcomb(L)].

% タプル別に必要なだけ組合せを生成する
comb2([{H,I}|J],K) when length(H) == K ->  %% 長さが求めるだけになった
    stripcomb([{H,I}|J]);
comb2(X,K) -> comb2( nexts(X), K ).  %% まだまだ必要

%% タプル第1要素が求める長さになっているので、第1要素だけのリスト（結果）にする
stripcomb(L) -> lists:map(fun(X)->element(1,X) end, L).

%% タプル要素ごとのループを回す
nexts([]) -> [];
nexts([H|T]) -> nexts2(H) ++ nexts(T).

%% タプル第1要素について、１つ長くなった組合せを生成し、リストで返す
nexts2({_,[]}) -> [];
nexts2({A,[L|M]}) -> [{A++[L],M} | nexts2({A,M})].

少々ややこしいですが、こんなところかな。最初アルゴリズムを確認するために、Scheme で書いたのですが、これ Erlang の方がずっと見やすく、理解しやすいです。特に「タプル」があって、見た目上リストと区別できるので、「書きやすい」印象が強いです。これで、

comb:comb( [1,2,3,4,5,6], 3 ).

みたいな呼び出しで、 ₆C₃ の組合せ（[[1,2,3],[1,2,4],...[4,5,6]]）が全部得られます。

では、今度は実際の検定です。測定自身は、使い捨て関数を書いて、引数で渡す、というかたちにしましょうか。

totalPerm:totalPerm(fun(X)->lists:sum(X)/length(X)>=25 end,[5,10,15,20,24,10,25,25,40,25],5).

これだと、実験群 5, 対象群 5 のデータに対して、実験群の平均値が 25 だったので、それを超える平均値の組合せの数と、その割合を返す、といった格好がいいでしょう。

{25,252,9.92063} %% {超える個数, 全体の数, パーセンテージ}

という感じですかね。では、実装は、

-module(totalPerm2).
-import(comb).
-export([totalPerm/3,calc/3]).

%% 外部公開。組合せを生成して下請けに回す
%% totalPerm( 判定述語, リスト, リストから「取る」数 )
totalPerm(Pred,L,K) -> X = comb:comb(L,K), tpLoop( Pred, X, length(X) ).

%% すべての組合せを試す
tpLoop(_, [], All) -> receiveResult(All,0,0); %% プロセス起動終りなので「受け取り」を起動
tpLoop(Pred, [H|T], All) -> invokeProcess(Pred,H), tpLoop(Pred, T, All ).

%% Erlang プロセスを起動する [自身のPid, 判定述語, 試す組合せ]
invokeProcess(Pred, L) -> spawn(totalPerm2,calc,[self(), Pred, L]).

%% Erlang プロセスで実行される関数。組合せを判定述語で評価し、
%% 結果を親に送り返す
calc(P, Pred, L) -> P ! evalList(Pred(L)).

%% 子プロセスで計算された結果を受け取る
receiveResult(All,All,Result) -> {Result, All, Result / All * 100}; %% 全部受け取った
receiveResult(All,Got,Result) -> receive Res -> receiveResult(All,Got + 1, Result + Res) end. %% まだあれば、更に受け取る

%% 述語での true, false を加算に変換
evalList(true) -> 1;
evalList(_) -> 0.

といったところでしょうか。意外にキレイに書けるものです....ちょっとびっくりするくらいです（まあ、他のマシンで依頼計算....まではしてないです。それでも大差はないですね)。このプログラムだと、Erlang のプロセスモデルを使い倒してます（それぞれの具体的な計算がプロセスで行われる）が、全然実行時にはストレスは感じません。もし、具体的な計算が重くなるようならば、分散させるのが良い手になってくるわけです....そういうあたりのバランス感覚が面白いことになるのでは？などとも感じました。