MathematicaのClusteringComponentsの困ったところ


Mathematica 9.0, 10.0, 10.1, 10.2, 10.3, 10.4.1, 11 for Microsoft Windows (64-bit)と10.0.0 for Linux ARM (32-bit)でのことです。

Mathematicaには、階層的クラスタリングができる関数が3つ用意されています。FindClustersAgglomerateClusteringComponentsです。

FindClustersにはバグがあります。

Agglomerateにもバグがあります。

バグではありませんが、ClusteringComponentsにも困ったところがあります。データをn個のクラスタに分けたいと思ってClusteringComponents[array,n]としても、できるクラスタがnより少ないことがあるのです。マニュアルには「最高でn個のクラスタを求める」とあるので、nより少ないのはバグでは無いのですが、ちょうどn個のクラスタを作りたいときに使えないのは困ります。

次のコードで再現できます(UMM3でも動きます)。

data = Import["https://gist.github.com/taroyabuki/4996086/raw/be3b2d537a51b803790fa1149cc714663a8b6ee9/clustering_test_data2.csv"];

Length[Union[ClusteringComponents[data, 13, 1, DistanceFunction -> EuclideanDistance, Method -> "Optimize"]]]
(* 12 *)

13個のクラスタを作りたかったのですが、できたクラスタは12個でした。

データをシャッフルしてからならうまくいきます。

Length[Union[ClusteringComponents[RandomSample@data, 13, 1, DistanceFunction -> EuclideanDistance, Method -> "Optimize"]]]
(* 13 *)

というわけで、階層的クラスタリングをしたいときはRを使うのがよさそうです(参考)。

数式を使わないデータマイニング入門


4334033555岡嶋裕史『数式を使わないデータマイニング入門 隠れた法則を発見する』(光文社, 2006)

最初の3章はデータマイニングとは何か、ビジネスでの利用例、データの準備・整理法といった総論的な話。その後はデータマイニング手法を一つずつ紹介する各論となる。紹介されている手法は、回帰分析・決定木・クラスタ分析・自己組織化マップ・連関規則・ニューラルネット。タイトルどおり、数式は使われておらず、トレーニングを受けていない人でも読めるようになっている。コーヒーチェーン店の自己組織化マップなど、データとアルゴリズムで再現できない例があるのが残念。

参考文献リストあり・索引なし

一様乱数から正規乱数を作る方法?


ソフトウェア関係の雑誌を読んでいたら、次のようなことが書かれていました。

「0以上1未満の一様乱数を12回足して6を引くと正規乱数になる」というテクニックを覚えておけば、一様乱数しか提供されていない開発環境で大いに役立つはずだ

これは中心極限定理のいい例ですね。実際に乱数を大量(下の絵では100万個)に発生させてそのヒストグラムを描いてみると、正規分布になっているように見えます(太線は正規分布)。これで十分な場合も多いでしょう。

しかし、正規性の検定をすると、正規分布であるという仮説を棄却したくなるようなp値が得られます。Mathematicaで試します(参考:DistributionFitTest)。

In[1]:= randomGaussian[] := Sum[RandomReal[], {12}] - 6

In[2]:= data = Table[randomGaussian[], {1000000}];

In[3]:= DistributionFitTest[data, NormalDistribution[0, 1]]

                  -8
Out[3]= 8.33717 10

この方法で生成される乱数の分布は、足す数(上の例では12)を増やせば正規分布に収束しますが、正規分布ではありません。ですから、「正規乱数になる」という表現はいけません。(例えば「randomGaussian[] := (Sum[RandomReal[], {24}] - 12)/Sqrt[2]」ならもっと正規分布に近づきますが、依然として正規分布ではありません。)

一様乱数しか提供されていない環境では、ボックス=ミュラー法などを使うのがいいと思います。Excelでは「=NORM.INV(RAND(),0,1)」が簡単です。

4874085601Pressほか『ニューメリカルレシピ・イン・シー』(技術評論社, 1993)にこんな記述がありました。

悪いrand()のため結果が疑わしいすべての科学論文を図書館の棚から消せば、各棚にこぶしが入るほどの隙間ができるであろう。(p.203)

「悪いrand()」ではなく「悪い乱数」を対象にすると、隙間はどのくらいになるのでしょう。

Pythonで階層的クラスタリング


Pythonで階層的クラスタリングをするときは、scipy.cluster.hierarchyを使うのが簡単なのですが、先日ちょっと間違えてしまいまして。

準備

scipyやmatplotlibを使います(インストールは済んでいるとして)。

import scipy
import numpy
import scipy.spatial.distance as distance
from matplotlib.pyplot import show
from scipy.cluster.hierarchy import linkage, dendrogram
from random import random

特徴ベクトルの長さがdimのデータをn個ランダムに作り、実験用のデータとします。

n = 100
dim = 10
data = [[random() for i in range(dim)] for i in range(n)]

方法1

非類似度の指標とクラスタ連結法を指定して、クラスタリングを実行します。

result1 = linkage(data, metric = 'chebyshev', method = 'average')

デンドログラムを描きます。

dendrogram(result1)
show()

このように、データをそのまま使ってクラスタリングするのが簡単です。

方法2

データ間の距離を与えてクラスタリングすることもできます。ただし、データ間の距離は、n行n列の距離行列ではなく、pdist()で作られるベクトルで表現します。

dArray1 = distance.pdist(data, metric = 'chebyshev')
result2 = linkage(dArray1, method = 'average')

assert (result1 == result2).all()#同じ結果

補足

pdist()で作られるベクトルは、次のようなものです。

dArray2 = []
for i in range(n - 1):
  for j in range(i + 1, n):
    dArray2.append(distance.chebyshev(data[i], data[j]));

assert (dArray1 == dArray2).all()#同じもの

このベクトルは、n行n列の距離行列からsquareform()で作ることもできます。

dMatrix1 = numpy.zeros([n, n])
for i in range(n):
  for j in range(n):
    dMatrix1[i, j] = distance.chebyshev(data[i], data[j])

dArray3 = distance.squareform(dMatrix1)
assert (dArray1 == dArray3).all()#同じもの

逆に、pdist()の結果からn行n列の距離行列を作ることもできます。

dMatrix2 = distance.squareform(dArray1)
assert (dMatrix1 == dMatrix2).all()#同じもの

まとめ

  • 階層的クラスタリングは、linkage(data)あるいはlinkage(距離のベクトル表現)で行う。(ちょっとわかりにくい仕様ですね。)
  • 距離のベクトル表現は、pdist(data)あるいはsquareform(n行n列の距離行列)で作れる。
  • n行n列の距離行列は、squareform(距離のベクトル表現)で作れる。

間違い

linkage(n行n列の距離行列)でクラスタリングするのは間違いです。こうすると、行列の各行がクラスタリングの対象データになってしまいます。

しかし、距離行列の各行は、元データのそれぞれが、他のデータとどの程度似ていないかを表現したものなので、これ自体を特徴ベクトルと見なすこともできます。ですから、次のような「間違った」方法も、自覚して使うなら間違いではありません。

result3 = linkage(dMatrix1, metric = 'chebyshev', method = 'average')

さらに、非類似度の指標としてチェビシェフ距離を使うなら、「正しい」方法と同じ結果になります。

assert (result1 == result3).all()#同じ結果