【1頁】
混合分布問題
――その基礎からカーネル降下法まで―― Part 1
金田 尚久、新居 玄武
１．はじめに
統計学で、通常用いられる確率分布は、一様分布を除けば、全て、一つの峰を持った分布である。ところが、諸科学の応用においては、二つ以上の峰を持った分布を考えねばならない場合がある。このようなとき、一番単純なモデルは、一つの峰を持った分布を必要な数だけそろえ、そのそれぞれにウェイトをかけて、たし合わせたものである。
)
ここで、一つの峰を持った分布 
の個数nがわかっているなら、それぞれの 
に含まれるパラメターの推定が問題になる。もしn がわかっていないのなら、n をどうやって推定するかが、問題として付け加わる。このような形で与えられた問題を混合分布問題と呼ぶ。19世紀の終わりにカール・ピアソンが、生物学の問題に関連して考察して以来、混合分布問題は、長い研究史を持つ。しかし、その歩みは、決して平坦ではなかった。まず第一に、ピアソンの最初の解法が余りに難解で、普遍性があるかどうかも疑わしかった。この方法では、問題は、積率法を通して、９次の代数方程式を解くことに帰着される。これは、コンピュータのなかった時代には、余りにも、大きな労力である。そこで、もっと計算の負担が少く、見通しのよい解法に改良する試みが現れた。しかし、これらの努力も、推定値の格別な向上には、つながらなかった。この後、研究は長い停滞の時期を迎える。新しい動きが見えたのは、コンピュータが普及し始める1960年代である。これ以後、この問題のために提案されたアルゴリズムは、コンピュータの使用を前提とし、くり返し計算に全面的に依存している。その中で代表的なのは、EMアルゴリズムとマルコフ連鎖モンテカルロ法 (MCMC) である。(いずれも、後章で詳しく説明する。）さて、初めに述べたように、混合分布問題には、n が与えられている場合と、推定しなければならない場合と、二種類ある。これまでの様々な応用例から見る限り、EM とMCMCは、n が固定されている場合には、まずまずの成果が挙がっている。しかし、nを推定する問題には、充分でない。ScottとSzewczyk
(SS) は、2001年に、この２つのアルゴリズムと全く異なる、新しいアルゴリズムを開発した［11］（文献はPart 2の最後にまとめて掲げる）。これは、１次元正規混合分布における、nの推定を目的としたアルゴリズムである。本論文の共著者の一人（金田）は、この方法に関心を持ち、そこで使われる新しいアイディアの解明と、同じアルゴリズムの他の混合分布への拡張を課題として、Ph.D.論文を執筆した［6］。本論文のトピックの一つは、そこで得られた成果を、初めて日本語で紹介することである。しかし、184ページと、かなり大部なPh.D.論文の要約版は、英文でpublish すべく、現在準備中である。本論文は、その要約版の単なる日本語訳ではない。むしろ、本論文においては、金田の得た成果の紹介は、最も主要なポイントに留め、混合分布問題の初歩から、最新の成果までのバランスの取れた解説を目標としている。次章でも述べるように、混合分布問題は、現在、理工学の幅広い分野で、熱烈な関心を持たれている研究テーマである。その簡明で一般的な問題設定から、我々は、社会科学方面にも大きな潜在的応用可能性を期待し得る。しかし、この問題は、現在までのところ、最大の応用分野が画像処理であることも手伝ってか、「先端工学」的なイメージが強い。社会科学への応用の試みが乏しいのも、その辺のところに理由があろう。そこで我々は、経済学部で通常教えられる、数学の知識を前提として、社会科学者のための、この問題のチュートリアルを提供することは、極めて望ましいと考えた。以下、このプランに沿って、混合分布問題の意味・EMアルゴリズム・MCMC・SSのアルゴリズム・金田による研究の順に解説する。
２．混合分布で，何がわかるのか？
1980年代半ばに出版された，Titterington,Smith,and Makov [12]
は，それまでの混合分布研究の総まとめとして，意義深く，2000年にMcLachlan and Peel[8] が出るまで，この分野の標準的モノグラフであった。そのp.16〜p.21には，Pearson以来その当時までの，応用例が，分野別の表としてまとめられている。これを見ると，医学・心理学・動物学・植物学・地質学・工学・ORなど，幅広い分野で，すでに混合分布が応用されていたことがわかる。応用分野は，その後も広がっているが，人文・社会科学の中で，比較的早くから混合分布の応用に熱心だったのは，マーケティングと心理学（精神医学）である。この二つの分野から，それぞれ一つずつ例を取って，混合分布を応用する動機を見てみよう。
例１．キャンディーの購買数に関する研究
DillonとKumar[2]
は，456人の消費者が，一週間に買ったキャンディーの数を調べた。このデータ全体を確率分布として把握するために，キャンディーの箱の数xを横軸に，
（x箱のキャンディーを買った人の数）/ 456 を縦軸にとる。この平面上にデータをプロットしたのが，下の実線と丸で表わされたグラフである。
【3頁】
このデータの形状は示唆深い。一定時間の顧客の数は、Poisson分布によって、よく近似される。
平均 
のとき、この分布は単調減少である。データは 
で最大値を取り、全体に右下がりの傾向を示している。このことは、この期間中に、1箱もキャンディーを買わなかった消費者が、多かったことを意味する。ところが、標準的な最尤法でパラメターを推定すると、
となり、平均は4に近くなり、一つの峰を持ったPoisson分布が現れてしまう。このことから、 のように見えるのは、見かけ上に過ぎず、この経験分布の右側のtailには、相当に確率質量が集積していることがわかる。さらに、小さな上がり下がりを追っていくと、このデータにはx=3,6,8,12及び15から17にかけてのところに、小さな峰がある。このような特徴をつかむには、いくつかのPoisson分布の重ね合わせと考えるのが適当であろう。DillonとKumarは、3章で説明する、EMアルゴリズムによって、混合分布をあてはめた。彼らは6つのPoisson分布の混合が最適と結論している。行列の第1列に平均、第2列にウェイトを書いて、6つのクラスターは、以下のように表せる。
図1の点線とアステリスクのグラフは、この6つのクラスターを、ウェイトを付けて足し上げたものである。最尤法であてはめたPoisson分布が、データの形状からかけ離れているのと比べて、この6クラスター・モデルは、はるかに良くなっている。しかし、x＝3の大きな峰がとらえられず、x=2のところに峰ができてしまっていること、6つのクラスターを使っているにもかかわらず、その他の小さな峰は、ほとんどとらえられていないこと、などを考えれば、まだ改善の余地はありそうである。
例2.
統合失調症の発症時期に関する研究
統合失調症（精神分裂病）は、多くの病気が征服された今日でも、解明の困難な病気である。ガンについては、発生機構に不明な点が多いにせよ、有効な対症療法は、数多く開発された。これに対して、統合失調症の対症療法には、一時的に患者の気分を軽快にする手段しかなく、その発生機構にいたっては、ほとんど何もわかっていない。しかし、19世紀からのデータの蓄積によって、この病気の多くの側面に、他の病気には見られない規則性が存在することが、わかってきた。それらの規則性には、データを取りまとめる際の判断基準にあいまいな点があり、そこを明らかにしないと、夾雑物が混じったままでしか提示できないといわれている。しかし、その存在に関しては、多くの精神科医が認めているので、将来の病気の本質解明に役立つ知見として、期待されている。そのような規則性の一つは、発症時期に関して、男女間に、共通点と相違点があることである。Lewine［7］によれば、統合失調症には早い時期（青年期）に発症する群と、遅い時期（壮年期）に発症する群と、二つのタイプがある。前者は主に男性、後者は主に女性の病気である。Everitt [3]
は、この仮説をデータによって検証した。図2aと2bでは、男女それぞれの発症時期のデータが棒グラフで表わされている。しかし、グラフを見やすくするために、横軸に年齢ではなく、log(年齢)を取り、これを等分の階級に分けて、それぞれの階級に落ちた観測値の数を、男性または女性の総数で割った値が、棒グラフの高さである。（ただし、図2a,bの説明参照）即ち、男女それぞれの経験分布が読み取れるように、調整されたグラフである。
【4頁】
図1の点線とアステリスクのグラフは，この6つのクラスターを，ウェイトを付けて足し上げたものである。最尤法であてはめたPoisson分布が，データの形状からかけ離れているのと比べて，この6クラスター・モデルは，はるかに良くなっている。しかし，x=3の大きな峰がとらえられず，x=2のところに峰ができてしまっていること，6つのクラスターを使っているにもかかわらず，その他の小さな峰は，ほとんどとらえられていないこと，などを考えれば，まだ改善の余地はありそうである。
例２．統合失調症の発症時期に関する研究
統合失調症（精神分裂病）は，多くの病気が征服された今日でも，解明の困難な病気である。ガンについては，発生機構に不明な点が多いにせよ，有効な対症療法は，数多く開発された。これに対して，統合失調症の対症療法には，一時的に患者の気分を軽快にする手段しかなく，その発生機構にいたっては，ほとんど何もわかっていない。しかし，19世紀からのデータの蓄積によって，この病気の多くの側面に，他の病気には見られない規則性が存在することが，わかってきた。それらの規則性には，データを取りまとめる際の判断基準にあいまいな点があり，そこを明らかにしないと，夾雑物が混じったままでしか提示できないといわれている。しかし，その存在に関しては，多くの精神科医が認めているので，将来の病気の本質解明に役立つ知見として，期待されている。そのような規則性の一つは，発症時期に関して，男女間に，共通点と相違点があることである。Lewine[7]によれば，統合失調症には早い時期（青年期）に発症する群と，遅い時期（壮年期）に発症する群と，二つのタイプがある。前者は主に男性，後者は主に女性の病気である。Everitt [3] は，この仮説をデータによって検証した。図2aと2bでは，男女それぞれの発症時期のデータが棒グラフで表わされている。しかし，グラフを見やすくするために，横軸に年齢ではなく，log（年齢）を取り，これを等分の階級に分けて，それぞれの階級に落ちた観測値の数を，男性または女性の総数で割った値が，棒グラフの高さである。（ただし，図2a,bの説明参照）即ち，男女それぞれの経験分布が読み取れるように，調整されたグラフである。
【5頁】
これを見ると，確かに青年期と壮年期に高い山があり，その間はくぼんでいる。しかし，Lewineが言うような，早発型は主に男性の病気であり，後発型は主に女性の病気であるという違いは見出せるだろうか？ むしろ，男性では高齢で発症する者がいるが，女性ではそのようなケースが少いということが，大きな違いではないかと思われる。そもそも，Lewineは
【6頁】“the early onset, typical schizophrenia is
largely a disorder in men, and late onset, atypical schizophrenia is largely a
disorder in women“
と言うのだが，“largely a disorder in men”とか“largely a disorder in women”とは，何を意味するのだろう？ 男性の分布の中で7:3女性の分布の中で4:6のように，早発性の男性と後発性の女性は，両方とも50％を越えているという意味だろうか？ それとも，男性の分布の中で4:6女性の分布の中で3:7のように，男女とも，早発・後発のいずれかに50％以上偏っている場合を許すのだろうか？ このようにLewineの主張にはアイマイな点があるが，Everitt は以下のように議論を進める。まず，男女それぞれのデータに2クラスター正規分布を当てはめ，次の結果を得た。
当てはめの方法はEMアルゴリズムである。後述するように、EMは、全パラメターに適当な初期値を与えた後、繰り返し計算によって、これらを最適な値に近づけていく、アルゴリズムである。Everittは種々の理由から、Pearsonの積率法が、この初期値の決定に利用できるとしている。しかし、これは今日、パッケージなどで一般に用いられている方法ではない。コンピュータの導入によって、Pearsonの時代よりも、計算の負担は楽に乗りこえられるようになったが、3クラスター以上の場合、多次元の場合に、この方法の拡張は、尚、適当でないからである。とにかく、このようにして得られた2クラスター・モデルの比較の対象として、Everittは、1クラスター・モデル（1つの正規分布）を、男女それぞれの分布にあてはめる（推定法と推定値は示されていない）。次に、2クラスター・モデルが本当に必要かどうか知るために、2クラスター・モデル対1クラスター・モデルの検定を行う。混合分布問題は、推定だけでも充分難しいので、検定はさらに難しい。この種の検定も昔から提案されているが、これには、大きな問題が指摘されている。しかし、Everittは、この検定の結果として、女性については2クラスター・モデルが適当であるが、男性については、1クラスター・モデルで充分であるとしている。Everittは、ロンドンの精神病研究所の統計学者であり、医療統計・応用統計の大家である。混合分布問題のモノグラフも出版している。しかし、このデータの分析のしかたは、疑問無しとしない。まず第一に、男性の標本平均は1.5236, 標本分散は0.1734 である。つまり、男性の早発型クラスターの平均・標準偏差と、男性全体の標本平均・標本分散とは、小数点以下4ケタまで一致している！さらに、ウェイトを見ると、早発型が90%を越えている。このような特徴は、早発型クラスターが、ほぼデータ全体を説明するように、推定されてしまっていることを意味する。一方、後発型クラスターの平均は、1.8267であり、その周辺に棒グラフの突出した形状は見当たらない。これで果たして、最適の2クラスター・モデルが得られたと言えるだろうか？このように無意味な付け足しとして、後発型クラスターが、推定されているならば、その必要性を検定によって、否定するまでもないのである。EverittはEM以外の方法で2クラスター・モデルの推定を、やり直してみるべきではなかったか？男女のデータには、確かに共通性がある。その一方によく当てはまる2クラスター・モデルが得られたなら、より良い推定法によって、男のデータにも良く当てはまる2クラスター・モデルが得られるかもしれない。11章で検討するように、EMは万能ではないのである。
というわけで、この例は、混合分布アプローチの有効性を示すには、不満足な点が残る。しかし、敢えてそのような例を引いたのは、同じデータを別の角度から見てみたいからである。より良い推定法によって、男の方からも、女と良く似た2クラスター・モデルが得られたと仮定しよう。次に考えるべきは、これらの2クラスター・モデルが、男女共通の2クラスター・モデルからの偶然の変異に過ぎないかどうか、ということである。我々がフォーマルな取り扱い方の確立しているデータを考察するのなら、まずそのような共通モデルの説明力について考察し、それでダメだとわかってから、男女それぞれの2クラスター・モデルを推定するだろう。しかし、今は、もっと探索的に考えていることとする。さて、そのような共通モデルを考えるためには、男女のデータを一つにプールし、そこにもう一度、2クラスター・モデルを当てはめる。そして、この男女共通のモデルが男女それぞれのデータをうまく説明するかどうかの検定を行えば良い。ところが、この検定は、通常の適合度のカイ2乗検定とは、異なる。通常のカイ2乗検定を行うなら、2つのクラスターをウェイトをつけて垂直方向に足し上げ、これを推定されたモデルとして、検定を行えばよい。ところが、このやり方では2つのクラスターを、せっかく識別した意味がなくなってしまうのである。我々は2クラスター・モデルとして適合しているかどうかを調べたいのだから、１つのクラスターが非常に良く適合していても、もう１つのクラスターが大きく外れていれば、２クラスター・モデルとして適合していないと判断すべきである。ところが、非常に良く適合しているクラスターのウェイトが高い場合には、通常の検定は、もう１つのクラスターが大きく外れていても、全体として適合しているという結論を出す可能性がある。この違いは、クラスターの数が多くなれば、より深刻になるはずである。では、そのように、１つのクラスターも大きく外れていないかどうか、より精密にチェックする検定は、どう構成したら良いのか？実は、これは、未解決の問題である。問題の設定は自然であり、男女差は人文・社会科学の到るところに現れる。男女それぞれの分布が、単純な確率分布で表わせない場合も多いだろう。このような検定が開発されれば、応用の範囲は、相当に広いに違いない。しかし、混合分布では、数理統計の通常の問題に見られるように、推定・検定・信頼区間が、密接に関連しながら発達するというスタイルで、議論が進んでいない。だから、これくらい自然でやさしく見える問題でも、解決を見ていないのである。実際のところ、Everittのアプローチの疑問点のように、推定を行う段階で、既に、克服しなければならない問題がある。そこで、我々も、本論文では推定の問題に集中する。
３．EMアルゴリズム
1章でも述べたように、今日、混合分布の推定に最もよく使われるのは、EMアルゴリズムとMCMCである。本章と次章では、クラスターの数を固定した枠組で、これらの推定法を述べる。EMについてはTitterington,
Smith and Makov[12]がMCMCについてはGilks, Richardson and Spiegelhalter [4]が、簡潔な展望を与えている。そこで、我々も、これらの本の記述に沿いつつ、より平易で、二つの手法の対比が明瞭な解説を試みる。尚、混合分布の中の一つ一つの小さな確率分布は、クラスターともコンポーネントとも呼ばれる。「クラスター」と言えば「分類的」な、「コンポーネント」と言えば「構成要素的」なニュアンスが伴うが、このような違いには余りこだわらず、適宜、同じ意味で使っていく。また、単純で紙数を費やす計算が現れるときは、細かい部分に立ち入らない。
観測値がn個あるデータに、kコンポーネントの混合分布を当てはめる問題を考えよう。この分布全体は
と表せる。それぞれのコンポーネントは、同じ型で異なるパラメターの確率密度関数
を持っている。j番目のコンポーネントのウェイトは
、パラメター・ベクトルは
である。全てのコンポーネントのウェイトとパラメターを、まとめて表現したいときは、文字
を用いる。確率変数
については、もし読者が、多次元分布の期待値の計算になじみがなければ、1次元と思って構わない。しかし、本章では、が何次元でも通用するように、記述を調整してある。そこで、対数尤度関数は、
ここで、現実のデータとしては得られない、仮説的な変数 
を設定する。各
はk次元の縦ベクトルである。の各要素は、
で、添え時は通常の場合と逆順である。はいわゆるindicator vector であり、
が第jコンポーネントから生み出された場合に、
。それ以外のの要素は全て0である。このようながすべての
に対応してわかっているのであれば、
と書けるだろう。ペア
の単なる代理として、
と書く。
に基ずく尤度関数は、
対数を取って、
ここに、
のj番目の要素は
、
のj番目の要素は
である。
Eステップ（期待値の計算）：
なる期待値の関数形を求めよ。
我々は、現在、イタレーションの第m+1段階を処理中であるとしよう。すると、
は前段階で得られた全てのパラメターということになる。
に含まれる
や
を未知のパラメターにしたまま、
の条件付き期待値を求めよ、というのが、ここでの問題の意味である。一見すると、そんなことは絶望的に見える。ところが、以下のように考えれば、この関数形は求まるのである。すでに求めた
から、
ここに、
この左辺はに対応するk次元ベクトルで、そのj番目の要素は、
最後のステップは、
第m段階における第jコンポーネントのウェイト
第m段階のパラメターを前提とした、第jコンポーネントのに
おける値
= 第m段階の全パラメターを前提とした、kコンポーネント・モデ
ルのにおける値
という式の意味による。このように計算を進めて行くためには、もちろん、第1段階のために、
を適当に与えておかなくてはならない。しかし、それさえ与えられれば、に関して何も確実な情報が無いにもかかわらず、その条件付き期待値は求まるというのが、トリックのポイントである。
Mステップ（最大化）： 
を最大化するを求めよ。
このように、各段階ごとに、EステップとMステップを繰り返しながら、EMアルゴリズムは進んで行く。即ちExpectation
Maximization Algorithmである。このアルゴリズムの段階を進むにつれて、尤度は単調増大することが、証明されている。しかし、大域的最大値は保証されない。むしろ、この尤度関数は、たくさんの局所的最大値を持っている。なぜならば、モデルに取りこまれた、たくさんのパラメターが、全て変数として扱われる、多変数関数だからである。初期値の選択に関しては、幾つかの提案がある。しかし、アルゴリズムによって、最終的に選択されたが、初期値
によって異なっていたという、報告は多い。このような問題点から、EMは完璧には程遠いと考えられている。
。
４．MCMC
MCMC（マルコフ連鎖モンテカルロ法、Markov
Chain Monte Carlo）は、混合分布へのベイズ的アプローチの代表的なものである。しかしながら、MCMCは、基本的な考えを、EMから借りているようなところがあり、その共通点とMCMC独自の点とを、明らかにしなければならない。そこで、前章で使用した文字も別の意味で再定義するので、注意して読んで欲しい。混合分布全体は
と書く。ウェイトは
、i とj は前章と逆の意味である。各コンポーネントは非常に広く、exponential family としよう。
と書けば、共役な事前分布が考えられ、
ここに
は定数のベクトルで、
と同じ長さを持ち、
はスカラーの定数である。ウェイト
の共役な事前分布はディリクレ分布
で与えられるとしよう。すると、その密度関数は
これらによって、事後分布は、
ここで、z変数を導入するが、前章とは異なる。ベクトル
において、
は、
が生成されたコンポーネントの番号を表しているとする。もし、このような
の値が全てわかっているのなら、
の式から非常に多くの項を省くことができる。
ここに 
は( )内の条件が満たされているとき1、そうでないとき0を取るindicator
function である。以上、準備したことから、次のようにイタレーションを構成する。
ステップ0：
に初期値を与え、これに基づいて
と 
を計算する。
ステップ1：パラメター・ベクトルとウェイト
を以下の分布形に基づいて生成する。
~ 
~ 
ステップ2：
を以下の分布形に基づいて生成する。
ここに 
ステップ3：と を更新する。
パラメターと
の生成に関しては、ベイジアンの枠組にふさわしい方式が開発されており、
Gibbs sampling という（Gilks, Richardson、Spiegelhalter[4]参照）。このイタレーションを続けていけば、以下の近似が充分よくあてはまる、イタレーション回数に達する。
ここで
は、イタレーション回数が、今述べた境界をM回越えていることを意味する。
, 
は、この意味での第i回イタレーションで生成された
と
である。
はとを変数とする任意の関数である。
ベイジアンの視点が一貫してはいるが、MCMCは内在的な問題から免れているわけではない。特に深刻なのは「ラベル転換問題」である。これまで仮定した通り、各コンポーネントは同じ型の確率分布に属している。この場合、事後分布は、コンポーネント番号（ラベル）を入れ替えても不変である。イタレーションの過程で、コンポーネント番号がふり替わってしまうことは、実際にあり得るのであり、しかも、このことはprocedure の収束を妨げるほど重大な問題である。我々はこの問題に深く立ち入らないが、まだ満足な解決を見ていないということは、言っておこう。
MCMCはEMよりも良いアルゴリズムなのだろうか？文献を調べてみると、この重要な問題に関して、わずかなことしかわかっていないのに驚く。現状を語るには、「競争そのものが、うまくformulate されていない。」としか言いようがない。行きづまりの原因として、少なくとも三つのことが、言えるだろう。第一に、混合分布の問題では、パラメターの数が非常に多く、どんなアプローチを取っても、計算の過程は長い。このような問題では、通常の理論統計学的基準は、直観的なアピールが乏しい。かと言って、それに代わる、混合分布にふさわしい、理論的基準が提案されたわけでもない。第二に、MCMCはベイジアン・アプローチの理論統計学者に好まれるだろうが、利用するユーザーにとっては、不便な点がある。EMでおかしな答えが出てきたときに、ユーザーは、別な初期値を試してみるだろう。しかし、ラベル転換問題の場合には、ユーザーが、その場に応じて処理できるような、解決策がない。第三に、シミュレーションに基づく比較研究も、積極的に行われたわけではなかった。統計の問題で、このような研究が手薄になるのは奇妙なことである。しかし、これには次のような事情があろう。混合分布は非常に大きなバラエティーを許容するので、どのくらいの難易度を、この分野のゴールとするのか、定め難い。そこで、異なった方法をテストするための、代表的な実例を選ぶのも困難なのである。
５．ScottとSzewczyk による新しいアプローチ
前二章で説明した，EMアルゴリズムとMCMCは，クラスター数の固定されたモデルの推定を，主要な目標としている。クラスター数を推定するには，これらのアプローチは，何等かの情報量基準に頼らざるを得ない。まず，あり得ると思われる，クラスター数の範囲を定める。【13頁】次に，その一つ一つについて，クラスター数を固定した推定を行う。この推定値に基づいてICを計算し，これが最適（普通，最小）となるクラスター数を選択する。従って，モデル・パラメターの推定と，クラスター数の選択は，基本的に別々の作業ということになる。Titterington, Smith, and Makov（1985）[12]とMcLachlan and Peel（2000）[8]に載せられている文献を比較してみれば，クラスター数を固定したモデルの推定に関しては，非常にたくさんの研究があるが，クラスター数の決定に関しては，それほど大きな進歩はなかったことがわかる。ICを用いたアプローチの改良が，その主なものである。さらに，前二章で述べた，EMとMCMCの難点に加えて，もう一つの問題を心に留めておかなければならない。それは，計算の実行速度である。最近では，高次元（例えば100次元）で大規模（例えば観測値数100,000）なデータ処理の需要が急増しており，EMの処理速度の遅さが，新たな懸念を生んでいる。クラスター数選択の実際は，単純に，一般的に述べることができる。procedureが，あらゆるクラスター数を通して，ある程度良いフィットを維持し，正しいクラスター数において，取り分け良いフィットを達成するなら，ICは，この正しいクラスター数を選ぶだろう。あらゆるクラスター数において，できる限りフィットが高くなるように設計されたアルゴリズムでは，クラスター数選択のときに，過剰な候補をかかえることになる。（このことの意味は，本論文の11章で明らかとなる。）これらの理由，即ち，既存の方法に残存する問題，大規模多次元データへの対応，クラスター数選択の基本的に簡明な構成は，我々をして，既存アプローチの彫琢よりもむしろ，単純化に向かわしめる。このような方向に最初に乗り出したのが，ScottとSzewczykの論文[11]である。
ScottとSzewczyk（SSと略記）は，彼らのprocedureを2001年に発表したが，その真価は，まだ充分に認められていない。これには，二つの理由が考えられる。第一に，既存の方法と比べて，彼らのアプローチが単純なように見えることである。そこで，混合分布問題（特にクラスター数の決定）を特別な難問と思っている人達に，見過ごされる傾向があった。しかし，実際は，本論文を読み進めばわかる通り，彼らのアルゴリズムと我々によるその拡張のパフォーマンスは，良好である。「これまでの研究者の期待を，はるかに超える」と言っても良いくらいである。SSの論文の第二の問題としては，多くの新しいアイディアが試みられているにもかかわらず，それらは必ずしも良く説明されていないということがある。彼らの新しい概念装置が，混合分布の本質を深くとらえているということを，納得するには，SSの原文をよりていねいに説明し，たくさんのシミュレーションを行わなければならない。当然のことながら，これら，克服すべき点を考慮に入れて，本論文の議論は組み立てられることになる。しかし，それをどうするかを述べる前に，我々は，SSのprocedureを簡単に紹介しておこう。ただし，Phase2は彼らのオリジナルそのままではなく，本論文の議論の構成にふさわしいように，変更が加えてある。こうしたのは，彼らはPhase2において，計算の能率を優先して，省略計算を取りこみ過ぎているからである。結果として，彼らのprocedureでは，Phase2の新概念がPhase3及び4の新概念と比較できないようになっている。我々のPhase2はより原始的（つまり，SSの元の意図に忠実）ではあるが，それ以後のphaseとの比較が容易になっている。本章ではまた，procedureの解説を要点だけに留める。新概念の意味を充分に納得するためには，例が必要である。しかし，2次元の実例は極めて咀嚼しやすいことがわかった。そこで，1次元procedureの詳細は，6-8章で2次元への拡張を考察する際に，もう一度掘り下げる。SSのprocedureは4部から成っている。
【14頁】
Phase1はカーネル推定である。我々は，どのカーネルも同じ大きさの包摂バンド（bandwidth, 「バンド幅」と訳されることがあるが，日本語として座わりが悪いので，「包摂バンド」と訳す）を持った，標準的なカーネル推定を行う。包摂バンドの選択法はクロス・バリデーションである。
ここにhは包摂バンド，はi番目の観測値，Nは観測値数である。このようにカーネルによって密度関数を推定した結果は，Nコンポーネントの混合分布とみなし得る。我々は以下のようにコンポーネント数を削減して，Nコンポーネント・モデルを最適なコンポーネント数のモデルに近づけて行く。削減プロセスは3つのphaseに分かれる。その最初のPhase2では，最も計算効率の高い方法が用いられる。まず，コンポーネント同士の相似性を相似測度によって評価する。
この測度は，2つのコンポーネント（どちらもpdfである）の位置と形を両方取り入れた「相似性」を測っていることを，注意しておこう。次に，相似測度の値が最高のペアを選び，これを1つのコンポーネントに融合する。融合法としては，積率法（Method of Moments, MMと略）を用いる。第一のコンポーネントはでウェイトは
，第二のコンポーネントはでウェイトはとしよう。すると，新しいコンポーネントのパラメターとウェイトは
【15頁】
SSは，この式が何に基づいているのか，明らかにしていない。しかし，我々は，7章において詳しく論ずる。ここでN段階は終了し，N-1段階に入る。そこでは，N-1コンポーネントをN-2コンポーネントに削減する，同じ作業がくり返される。このようにして，クラスター数nを下って行き，Phase2と3の境界のnに達する。この境界とPhase3から4への，もう一つの境界は，アド・ホック的に選ばれている。しかし，SSの例と本論文のシミュレーションから見る限り，これらは適切に選ばれている。Phase3では，融合法として，MMを残すが，相似測度を変更する。今，nクラスター・モデルを処理中で，n-1クラスター・モデルに減らすところだ，としよう。簡単のために，隣り合ったコンポーネント（隣り合っているかどうかは，例えば，各コンポーネントの平均の位置で決める）のみを考え，以下のやり方で，融合するのに最も良いペアを選ぶ。まず我々は，コンポーネント番号1と2を1つに融合する。nクラスター・モデルの内，番号1,2以外の全てのコンポーネントを残し，これに融合したコンポーネントを加えて，n-1クラスター・モデルの候補が一つできあがる。次に，nクラスター・モデルと，この候補モデルの間の相似測度を計算する。ここで用いる相似測度とは，nクラスター vs. n-1クラスター相似測度である。個々のクラスター同士の相似測度に似せて，新しい相似測度はと定義される。ここにλはnクラスター・モデル全体の，γはn-1クラスター・モデル全体のpdfである。次に，もう一度nクラスター・モデルに戻り，そのコンポーネント番号2と3を1つに融合する。コンポーネント番号2,3以外の全てのコンポーネントを残し，もう一つの候補モデルを作る。それから，相似測度を計算する。このようにして進んで行き，全てのペアについてやり終えた後，nクラスター・モデルとの間の相似測度が最高となる候補モデルを，n-1のクラスの最適のモデルとして選択する。
詳しい説明抜きで，ザッとPhase3を眺めただけでも，Phase2と3の違いは何なのか，1 vs. 1相似測度と，n vs. n-1相似測度の違いは何なのか，ということに戸惑いを覚える。ところが，本研究によって，この違いは，極めて重要であることが，わかった。そこで，実際のprocedureの過程の中で，この違いをどうやって確かめるか，それをどう説明し，理論的背景をどう確立するかという問題は，本論文の焦点の一つとなる。
【16頁】適当な境界を越した後，Phase4が始まる。ここでも我々は，クラスター数nを降下していくプロセスを続けるが，もう最適のクラスターに近いと考えられるので，計算時間を費やす方法を用いてもさしつかえない。n vs. n-1相似測度は，そのままに残すが，融合法はL2E最適化に置き換える。L2EはDavid W. Scottによって新しく提案された，パラメター推定の原理である。今，考えている混合分布の枠組では，outlierに対する頑健性を，この方法のメリットと，彼は見なしている。積分において，は推定すべき密度関数，fは真の密度関数，θは推定の目標たるパラメター・ベクトルである。L2Eはこの積分を最小化するを求めることでfを推定しようとする。
この式は真の密度関数を含んでいるので，直接には操作可能でない。しかし，によって近似できることが，わかっている。さらに，正規混合分布の場合，一項目は，積分記号の入らない閉じた形に変形できる。今，紹介している論文とは別に，Scottは，L2Eのためのもう一編の論文[10]を著わしているので，その理論的な性質については，そちらを参照して欲しい。Phase4では，各クラスター数において選ばれたモデルから，情報量基準を計算する。これに基づいて，最適なクラスター数を選び，そのクラスター数で選ばれていたモデルが，最適なモデルである。
[11]において述べられている，その他の結果を，簡単に要約しておこう。彼らは，このprocedureを，コンピュータで実行する上での，細かい問題点について論じている。ただし，1次元のみである。経済学および天文学のデータへの応用を試みている。仮説例から発生させたデータ・セットについては，三つの例を挙げている。ただし，これらは，うまくいった例だけである。一つの例から，たくさんのデータ・セットを発生させて行う，繰り返し実験は，やっていない。高次元データについても，試みたようであるが，方法の詳細は示されていない。彼らのprocedureが，どれくらいうまくいったかについての，定量的な評価は，充分に行われていない。繰り返しシミュレーションによるデータ・セット中，いくつに対して，正しいクラスター数が得られたか，成功率を示す必要があろう。このことは重要である。混合分布のような分野では，全ての方法を評価する基準が確立していない。個々の方法の得失を知るために，多様なシミュレーションを試みなければならない。
SSのprocedureをカーネル降下法（Kernel Reduction, KR）と呼ぼう。彼らの達成したことを踏まえて，本論文では，二つの新しい課題を設定する。第一は，彼らが充分に説明できなかった，新概念の意味を明確にすることである。第二は，1次元正規混合分布のためのprocedureを，2次元正規混合分布に拡張することである。すでに述べたように，2次元の例はわかりやすいので，これからの叙述は，第二の課題を答える順通りとし，第一の課題は，その中の適当な所で答えることとする。尚，Phase3は，データ・サイズもクラスター数も少ない場合には，省略【17頁】できることがわかった。そこで，本論文ではPhase3を飛ばし，Phase4をPhase3と呼ぶことにする。この3 phase procedureはKRの基本的なアイディアを全て含んでいる。
６．2次元のためのPhase1（カーネル推定）
SSのオリジナル・アプローチでは，全てのカーネルに共通の包摂バンドを仮定して，カーネル推定が行われる。通常，これを2次元に拡張するときは，3つのパラメターを推定する。（各カーネルのは観測値によって置き換えられる。）しかし，我々は，これよりも単純なアプローチを試みる。即ち，と置いて，のみを推定するのである。その理由を述べよう。真のモデルが，数個のクラスターから成っていると仮定する。それらの全ては，一つを除いて，右肩上がりの等高線（正の相関）を持っているとしよう。例外のクラスターは左肩上がりの等高線（負の相関）を持っている。すると，カーネルもまた，右肩上がりと推定されるに違いない。このことは，ほとんどのクラスターの降下過程にとって，有利な背景となる。しかし，一つだけの例外のクラスターにとっては，問題を生じるだろう。降下過程を進むに従って，クラスターの本来の形が形成されてくる。これは，例外のクラスターの場合，もちろん，左肩上がりである。しかし，形成されつつあるクラスターと，残りのカーネルの間の相似測度の値は小さい。何故ならば，形が違うからである。この困難は，高次元では，更に大きくなる。高次元では，より多くの直交性が存在するのだから，カーネルが，統合されるべきクラスターに対して，直交に近く向いている可能性も大きくなる。カーネル推定のたいていの文献で，ρは推定されているが，その意味について，特に筋の通った解釈がなされているわけではない。どんな場合にρをモデルに含み，どんな場合に含めないかということは，もっと注意を向けられて良い問題である。
以上述べてきた推定法は，最も標準的なカーネル推定を，若干修正したものにすぎない。だから，カーネルについて良く知っている読者には，詳しく展開する必要は，あるまい。しかし，我々の降下procedureが，どんなところから始まるかを，明らかにするために，数学的な形を示そう。
今，N個の観測値があるとしよう。カーネルによって推定された確率分布は，で与えられる。は，観測値の「上に乗っている」カーネル，即ち，平均の2次元正規分布である。xとyの両方向に等しい分散は，hで表わされている。ベクトルxとをと書こう。我々は,に関する単純化の仮定を持っているから，包摂バンドの決定には，最小二乗クロス・バリデーション（Least Squares Cross-Validation, LSCV）が最も良く用いられる。LSCVは二乗誤差積分平均（Mean Integral Squared Error, MISE）の最小化を目指している。
ここには真の密度関数である。最後の項は最小化すべきパラメターを含んでいないから，最初の二項のみを問題にする。以下の式が，最初の二項の不偏推定量を与えることがわかっている。
LSCVの一項目は，2次元正規分布の積の積分を含んでいるが，これは以下の公式（Wand and Jones [14] 参照）によって簡単に求まる。
Phase3では，最適なクラスター数選択のために，各nにおいて選ばれたモデルの情報量基準を計算する。我々のアルゴリズムではAICとBICを計算し，そのモデル選択における得失については，次章で検討する。既に予告した通り，MEASURE1と2の違いは重要である。しかし，これはシミュレーション結果を見ることによって，最も効率的に説明できる。そこで，次章の後に独立の章を設けて議論する。