MathematicaのClusteringComponentsの困ったところ

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

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

FindClustersにはバグがありました。(FindClustersのバグは11で解決)

Agglomerateにはバグがあります。

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

次のコードで再現できます。

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()#同じ結果

つながりやすさ、No.1へ?

スマートフォンのつながりやすさに関する調査で、ソフトバンクが1位になったという広告が大量に流れているようです。

なんでも、NTTドコモとau、ソフトバンクのそれぞれ3400, 3700, 5300人に、月間15回、7時から23時までの時間帯でランダムに発信した場合の通話接続率が、順に0.98, 0.981, 0.983だったとのこと。1位になったのはこれが初めてではないようですが、その瞬間にはCMの準備ができていなかったのでしょう。

違いが微少なので、誤差の範囲では?とも思うのですが、ちょっとRで計算してみると、そういうわけではなさそうです。

> trials <- c(3400, 3700, 5300) * 15
> successes <- trials * c(0.98, 0.981, 0.983)
> prop.test(successes, trials, correct=F)

        3-sample test for equality of proportions without continuity
        correction

data:  successes out of trials 
X-squared = 16.9408, df = 2, p-value = 0.0002096
alternative hypothesis: two.sided 
sample estimates:
prop 1 prop 2 prop 3 
 0.980  0.981  0.983 

4863540930石田基広『R言語逆引きハンドブック』(シーアンドアール研究所, 2012)初版p. 386のサンプルは、成功回数と試行回数が逆になっていますね。

とりあえず誤差ではないとしても、調査方法でよくわからないところもあって、なんとも言えない感じですが、「ちょっとうれしいんで、言わせてやってください」と言われると・・・

「つながりやすいソフトバンク」はスマートフォン限定(Slashdot)ソフトバンクの「つながりやすさ、No.1へ」に SNS からは厳しい声も(japan.internet.com)など、言わせてやってあげられない方面からはいろいろツッコミもあるようですが、それらの中にもツッコミどころがあったりして(「本来、このような品質調査の場合は、比較する3社でサンプル数を揃えて違いを探るものである」とか。そもそも「SNS から」っていうのもおかしい)、広告主にとってはかなりうれしい展開になっているようです。

このブログもつながりやすさの向上を目指しています。