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を使うのがよさそうです(参考)。

縦と横にしか動けない世界で(1994年東京大学入学試験理系数学第6問)


1994年の東京大学の入学試験、理系数学第6問は次のようなものでした。

平面上の2点P, Qに対し、PとQをx軸またはy軸に平行な線分からなる折れ線で結ぶときの経路の長さの最小値をd(P, Q)で表す。

(1) 原点O(0, 0)と点A(1, 1)に対し、
d(O, P)=d(P, A)を満たす点P(x, y)の範囲をxy平面上に図示せよ。

(2) 実数a>=0に対し、点Q(a, a^2+1)を考える。
次の条件(*)を満足する点P(x, y)の範囲をxy平面上に図示せよ。
(*) 原点O(0, 0)に対し、d(O, P)=d(P, Q)となるようなa>=0が存在する。

この問題は、d(O, P)=Abs[x]+Abs[y]であることがわかれば解けます。d(O, P)=Sqrt[x^2+y^2]ではありません。

Mathematicaで試します(UMMでも動きます)。

(*1*)
d[p_, q_] := Total[Abs[p - q]]
expr = d[{0, 0}, {x, y}] == d[{x, y}, {1, 1}];
cond = Reduce[expr, {x, y}, Reals];
reg = ImplicitRegion[cond, {x, y}];
RegionPlot[reg, PlotRange -> {{-2, 2}, {-2, 2}}]

(*2*)
expr = Exists[a, a >= 0,
   d[{0, 0}, {x, y}] == d[{x, y}, {a, a^2 + 1}]];
cond = Reduce[expr, {x, y}, Reals];
reg = ImplicitRegion[cond, {x, y}];
RegionPlot[reg, PlotRange -> {{-2, 2}, {-2, 2}}]

描画領域の境界にも実線が引かれている、という問題がありますが、とりあえずはこれでいいでしょう。

数学ガールのミルカさんが「こうなる」とひと言発する間にやっていること(ウラムの螺旋)


4797374152結城浩『数学ガールの秘密ノート 整数で遊ぼう』(SBクリエイティブ, 2013)(索引あり・参考文献リストあり)にウラムの螺旋が出てきます。

100くらいまで描いたところで、

「あ、あたしっ、この先にも興味があります。ぐるぐるぐるぐるぐるぐるぐる……とずっと続けたらどんな図形が表れるんでしょう」

と言うテトラちゃんに応えて、

「こうなる」

と言ってミルカさんがかなり大きなウラムの螺旋を見せています。(p.70)

「ウラムの螺旋」で検索すればぱっと出せそうな気もしますが(インタラクティブなものなんかも)、本書で描かれているのは1から始まるよくある螺旋ではなく、0から始まる螺旋なので、そうでもないかもしれません。そもそもミルカさんは、「こうなる」といいながらウェブ検索をするキャラではありません。

ですから彼女は、「こうなる」とひと言発するわずかの時間でウラムの螺旋を描くスクリプトを書いているんだと思うのです。(彼女はノートPCを持ち歩いている設定でしたっけ?)

そんなことができるのかと思って最初に思いつくのはこんなスクリプトです(Mathematica)。UMMでも動きます。

spiral[size_, start_] := With[{center = Ceiling[(2 size - 1)/2]},
  Module[{
    m = Table[start, {2 size - 1}, {2 size - 1}],
    i = start},
   Do[
    Do[m[[r, center + k]] = ++i, {r, center + k - 1, center - k, -1}];
    Do[m[[center - k, c]] = ++i, {c, center + k - 1, center - k, -1}];
    Do[m[[r, center - k]] = ++i, {r, center - k + 1, center + k, 1}];
    Do[m[[center + k, c]] = ++i, {c, center - k + 1, center + k, 1}],
    {k, 1, size - 1}];
   ArrayPlot[PrimeQ[m] /. {True -> 1, False -> 0}]]]

これで「spiral[4, 0]」などとすればウラムの螺旋を描けますが、会話のテンポはかなり遅くなるでしょう。

会話のテンポを自然なものに保つためには、もっとコンパクトなスクリプトを書けなければなりませんが、すぐには思いつきません。(かつて書いたスクリプトがノートPCに保存してあった?)

というわけで、彼女は優秀な高校生だという物語の設定は、そのとおりだと思いました。

会話形式で書かれている数学ガールには、会話のペースに合わせてできるものだろうかと考えながら読む楽しみがあります。

せっかくなので、サイズを大きくしたときの変化を見てみましょう。

Animate[spiral[size, 1], {size, 1, 100, 1}]

サイズを大きくしたときの変化

中心の数を変化させたときの様子も見てみましょう。

Animate[spiral[50, start], {start, 0, 100, 1}]

中心の数を0から100まで変化させた様子

Wolfram CDF Playerがインストールされていれば、インタラクティブに調べられます。


MathematicaのNMinimize, NMinValueのバグ


Ver.11で修正されました。

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

a2 + c2 == 1, b2 + d2 == 1という条件の下で、次のような関数の最小値を求めたいとします。

Mathematicaには最大・最小を求めるためのさまざまな関数が用意されていますが、NMinimizeMinimizeNMinValueMinValueは、頭にNをつけるかどうかで、数値的な方法と解析的な方法を切り替えられて便利です。最小値とその時の各変数の値を知りたいときは[N]Minimize、最小値だけを知りたいときは[N]MinValueを使います。

しかし、ちょっとうまくいかない例に遭遇しました(UMMでも試せます)。バグが修正されたバージョンを受け取るために、プレミアユーザになっています。

問題その1

f1 = Sqrt[
 1 + 4 Sqrt[3] b c d + d^2 - 5 c^2 (-1 + d^2) + 
  2 a (-Sqrt[2] b d + Sqrt[6] c (-1 + d^2))];

NMinimize[{f1, a^2 + c^2 == 1 && b^2 + d^2 == 1}, {a, b, c, d}]

このように実行すると、

関数の値4.35271 +1.45672\ Iは{a,b,c,d} = {0.79784,1.52261,-0.444508,0.634251}において実数ではありません.

という警告とともに結果{1.22618, {a -> -0.82781, b -> 0.492467, c -> 0.224354, d -> -0.333602}}が返ります。このとき、a2 + c2は約0.74、b2 + d2は約0.35なので、与えた制約条件から大きくずれています。

問題その2

f2の最小値を求めようとしたときにも問題が起こります。

f2 = Sqrt[
  1 - 4 Sqrt[3] b c d + d^2 - 5 c^2 (-1 + d^2) + 
   2 a (-Sqrt[2] b d - Sqrt[6] c (-1 + d^2))];

NMinimize[{f2, a^2 + c^2 == 1 && b^2 + d^2 == 1}, {a, b, c, d}]

このように実行しても、入力がそのまま返ってくるだけで、まとも(?)な結果が得られません(間違った答えが返ってくるよりはましかもしれません)。

この問題は、WolframAlphaでも起こります(2015年5月17日確認)。

f1に関しては間違った結果が返ってきます(下はNMinValueの場合。NMinimizeでも同様)。

f2に関しては、珍しいことに、WolframAlphaは完全に沈黙します(下はNMinValueの場合。NMinimizeは違う結果になります)。

MathematicaからWolframAlphaに問い合わせたときはまた違う結果で、f2は解釈できないとのことでした。ほとんど同じf1は解釈できたので、そんなことはないと思うのですが。

そもそも、こういう関数の最小値を何の工夫もせずに得ようとするのが間違いなのかもしれませんが、あえて工夫をしないのにも理由があるのです。それについては別の機会に書きます。

MathematicaのAgglomerateのバグ


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

Mathematicaには、階層的クラスタリングができる関数が3つ用意されています。FindClustersAgglomerateClusteringComponentsです。3つもあるのが問題な気もしますが、とりあえずそれはよしとして・・・

FindClustersにはバグがあることを以前報告しました。

ですから、階層的クラスタリングをしたいときは、AgglomerateかClusteringComponentsを使わなければなりませんが、残念なことに、Agglomerateにもバグがあります。FindClustersの場合と同様、このデータと次のコードで再現できます(RSSリーダーでは表示されないかもしれません)。Assertによる警告は表示されませんが、UMMでも動きます。

クラスタリング結果にはすべての要素(この例では63個)が入っていなければなりませんが、結果を数えてみると足りません(この例では38個)。このバグは報告済みですが、解決方法はわかっていません。バグが修正されたバージョンを受け取るために、プレミアユーザになっています。

FindClustersとAgglomerateが使えないとなると、残るはClusteringComponentsだけになるわけですが、これにもちょっと問題があります(MathematicaのClusteringComponentsの困ったところ)。

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