2015/02/18

継承は必ずしもis-a関係ではない(OOP) - Inheritance is not subtyping問題

Inheritance is not subtyping」は, Cookらによる1989年の論文で, そのままですが, オブジェクト指向における継承≠派生型(is-aの関係)ということのようです.

 派生型(Subtyping)については, is-a - Wikipediaにその説明が載っていますが, 確かにここでも, 継承が必ずis-a関係を構築するとは書かれてはいません.

 数理科学的バグ撲滅方法論のすすめ - 第4回 関数型言語とオブジェクト指向,およびOCamlの"O"について - ITpro に, その反例を示したソースコードがありますが, OCamlなのであまり実感わきません. オブジェクト指向言語は, JavaかPythonくらいしか使ったことがないので, Pythonでとりあえず書き直してみました.

 Pointクラスを継承したColoredPointクラス.
class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def equals(self, another):
        return self.x == another.x and self.y == another.y

class ColoredPoint(Point):
    def __init__(self, x, y, c):
        self.x = x
        self.y = y
        self.c = c
    def equals(self, another):
        xeq = self.x == another.x
        yeq = self.y == another.y
        ceq = self.c == another.c
        return xeq and yeq and ceq
 この例では, 確かにPointクラスが継承されていますが, equalsがオーバーライドされています. 一方で, Pointクラス間ではequals関数で比較できるように実装されています.

 PointクラスとColoredPointクラスのインスタンスがis-a関係ならば, 当然, Pointクラスのインスタンス⇔ColoredPointクラスのインスタンス間で, インスタンスのequalityの比較が可能なはずです.
>>> ColoredPoint(1,2,3).equals(ColoredPoint(1,2,3))
True
>>> ColoredPoint(1,2,3).equals(Point(1,2))
Traceback (most recent call last):
  File "", line 1, in 
  File "notisa.py", line 16, in equals
    ceq = self.c == another.c
AttributeError: Point instance has no attribute 'c'
 で, 比較してみると, エラーが発生する. という現象が反例としてあげられるため, 継承≠is-aということのようです.

 面白いのは, 継承しなくても, is-aの関係が成立してしまうケースが存在することで, 次のFruitとOarngeクラスの例が考えられます.
class Fruit:
    def __init__(self, price):
        self.price = price
    def equals(self, another):
        return self.price == another.price

class Orange:
    def __init__(self, price):
        self.price = price
    def equals(self, another):
        return self.price == another.price
厳密な意味(?)でis-a関係が成立しているわけではないですが, equalsによる比較ができてしまい, 一種の部分型(is-a関係)のような振る舞いが可能になります. OCamlなどでは, この辺を型推論により自動的に型付けしてくれるようで, このような型を構造的部分型(structual subtyping)と呼ぶようです(前述のITProの記事より).

2015/02/14

quilを使って巡回セールスマン問題の解の経路を表示する(Clojure) + 矢印の表示


 は, ProcessingのClojureラッパーです. グラフィック系のライブラリで, Clojureから, ProcessingのAPI(メソッド等)を簡単に触ることができます. もともと, Processing自体が非プログラマを対象としたヴィジュアルアート向きの, Javaベースのプログラミング言語なので, その機能をClojureで使えるようにしてくれます. Clojure/Swingでやるには少し煩雑なグラフィック処理を, Processing風の書き方で書くことができます. マウスやキーボードの情報の取得も可能なので, ゲームも作れそうです.  以下がAPI(関数)の一覧.


 少し前に, 巡回セールスマン問題を遺伝的アルゴリズムを使って解いたのですが, その時の解(となる経路)の可視化がいまいちだったので, quilを使ってやり直してみました. もともとは, Processing同様, quilも 次のquil.infoにあるようなヴィジュアルアートを作成することを目的としているようですが, 容易に絵が描かけるという点には変わりありません.

 ちなみに, quilのコードを記述する部分はほぼ100%(?)手続き型言語風の書き方になります. map等でループを書くことはできますが.

 他のClojureのプラグイン同様, プロジェクトを作成しているなら, project.cljの :dependenciesに[quil "2.2.5"]を追加(quilでは試していないため保証はできませんが, 一応これ) でOKのはずです.
 lein repl 単体ならば, "~/.lein/profiles.clj" (または, それに相当するファイル)の:pluginsに対応するリストへ  [quil "2.2.5"] を加えれば利用できるようになります.

 実際にプロットさせると次のようになりました. (ソースコードは一番下)

 上記の表示は,
(def sample-points
  {:a [20 70] :b [30 20] :c [50 50] :d [60 55] :e [40 30]
   :f [20 50] :g [80 60] :h [10 30] :i [65 10] :j [40 80]})
のようなXY上に点がばらついているときに, 次のような解が生成された時に, 表示されるルートです.
(def route [:a :f :j :d :g :i :c :e :b :h])
 単に線を引いて, 丸を描いて, 文字を打ち込んでいるだけです. なので非常に簡単でした. しかし, Swingではなく, quil独特(Processing特有?)の関数名がいくつかあり, 目的の関数を探すのに時間がかかりました.

 API一覧はあるのですが, 例えば, 丸を描く関数名がellipseだったり(SwingだとfillOvalとかですが), その丸を描く関数を使うときに円の内側の色を指定する関数名が, fillだったりします. 線の太さは, stroke-weightで設定します.

 一方で, Processingのラッパーだけあっって, ネット上にあるProcessingのサンプルをClojure/quilに読み替えれば, 様々なことができるので, その点は非常に良い点です. quilのAPI一覧を見てよくわからない場合は, ProcessingのAPIやサンプルを見れば良いわけですから.

 ウィンドウのみを表示させたければ, 次のように書けます.
user> (require '[quil.core :as q])
nil
user> (q/defsketch skt2)
#'user/skt2
という形で表示できます. このdefsketch関数のオプションで, 描画の設定等を行います.
(defn plot-route []
  (q/defsketch skt1
    :title "routes"
    :setup (fn [] (q/smooth)) ;; anti-aliased
    :draw (draw-route sample-points route)
    :size [400 400]))
のようにdefsketchを定義します. タイトル等はオプションです. :setupは, ウィンドウ生成時に, 実行される関数, :drawは, 定期的に実行される関数ということでしょう. smooth関数でアンチエイリアスが指定できます. 上記の例を見ればわかるかとは思いますが, 滑らかな表示に切り替わります. :sizeのパラメータは, ウィンドウのサイズ.

 こういう経路情報の表示で意外と面倒なのが, 矢印の表示です 点(x1, y1)と点(x2, y2)間に線を引いた時, 頂点の片方に矢印の三角形を付け加えて矢印にする作業です. そして, quilにも多分矢印を表示する関数は無かったと思います.

 仮想的なXY空間を用意し, そこに矢印の先となる三角形を想定した3つの点をおいて, 回転行列/アフィン変換等で適当な角度に回転させます. 次にその3つの点に必要なdrawPolygonメソッド等で描画するというステップが必要になります.

 色々ググって調べてみたのですが, Swing/Awtでは, やはり上記の手法が一般的なようでした. 回転の処理を直接書くのは面倒なのですが, Processingには, matrix stackというオブジェクトがあり, これを使うことで簡単に三角形のオブジェクトを回転させることができることができるそうです.

 これで回転の処理の問題が解決したので, あとは与えられた座標から矢印の先端を表す三角形を回転させる角度がほしいのですが, これはatan(アークタンジェント)で解決します. しかし, ゼロ除算の問題でatanをそのまま使うよりも, atan2という(x, y)の値を渡すとその角度(単位はラジアン)を返してくれる関数の方が有用です. 知らなかったのですが, atan2は, MayaのMELエクセルとか, C++にも標準でついてくるような一般的な関数のようです.
(defn draw-arrow-head [xy1 xy2 top-mergin bot-mergin]
      (let [x1 (first xy1) x2 (first xy2)
            y1 (second xy1) y2 (second xy2)
            size (+ 25 (- bot-mergin))]
        (set-white)
        (q/push-matrix)
        (q/translate x1 y1)
        (q/rotate (+ (q/radians 270) (q/atan2 (- y2 y1) (- x2 x1))))
        (q/triangle 0 (+ top-mergin 0) 5 size -5 size)
        (q/pop-matrix)))
上記のような実装で, 矢印の頭の部分が書けます. 矢印の本体の線は, line関数で書けばよいので, 特に問題はないかと思います.

 rotateは, 時計回りにオブジェクトを回転させるようです. GUIプログラミングではお馴染みですが, Y軸が反転しているので, その反転分の270度を足すことで, 全体として, 適切に傾いた矢印が描けます. 以下の図が参考.
 atan2は, x軸に対する(x, y)座標の角度θを返します. translateには, オブジェクトの原点に相当する点の位置を与えます. 矢印の三角形を原点から少し離しているのは, 原点のまわりに各点を表すノードを表す円(ellipse)を表示させるためです. こんな感じで, ノード間の矢印を表しました.

 quilを使っていて少し困ったのは, エラーメッセージがreplに出力されないことです, 一度defsketchでオブジェクトを生成されてしまうと, :setup/:drawに指定した関数内で実行時エラーが発生した場合, エラーメッセージが出力されず, エラーが発生しているのかすら, 分かりませんでした. これはおそらく, defsketchを実行した時点で別プロセス(スレッド)でその処理が実行されるため, replのあるスレッド上ではエラーが発生していないことに由来するようですが, デバッグが大変なので, 対策を考える必要がありそうです.

 全ソースコードは以下から.
uid0130 / plot-graph.clj - Gist

参考文献

2015/02/10

ClojureでリストのHylomorphismの変換マクロ

 前に書いた記事(Clojureでunfoldとhylomorphism) の続き.

 上の記事では, ClojureにおいてList版のHylomorphismを中間のデータ構造を生成しない純粋な再帰呼び出しで書きなおすことで高速化しました. しかし, 可読性が最悪で, はっきりいって読めない...
(defn correct-answer? [arrangement]
  (reduce andf
    (map not-conflict?
      (take-while first (iterate rest arrangement)))))

(defn correct-answer?-v2 [arrangement]
  (hylol andf true #(not (first %)) not-conflict? rest arrangement))
 上の関数が, リストを生成しながら計算するHylomorphism(に相当するもの), 下の関数が, リスト(中間のデータ構造)なしで同じ計算をするものです. (iterateで)リストを生成せずにプログラムが実行される分, 下の関数の方が高速であるという結果になったという話でした.

 hylolとはどいういった関数かというと, 次のような関数です.
(defn hylol
  [f e p g h b]
  (if (p b)
    e
    (f (g b) (hylol2 f e p g h (h b)))))
 ピカソの絵画ような抽象的な定義ですが, しかし抽象的であるがために, ある程度普遍性があります(様々なプログラムに出現する可能性が高いという意味です).

 できれば可読性を維持したまま(つまり, reduceとかmapとかが書かれた状態のまま) Hylomorphismに相当するプログラムを純粋な再帰に変換したいですね. というわけで, 関数でhylolを実装するのではなく, マクロで実装すれば良いのでは? という結論に至りました.
(hylol-expand
  (reduce + 0
       (take-while (partial not= 0)
                    (map dec (iterate dec 100)))))
とか
(hylol-expand
 (->> (iterate dec 100)
      (take-while (partial not= 0))
      (reduce + 0)))
みたいな書き方なら, 可読性に関しては問題ありません.

 というわけで上記のような書き方を維持したまま, Hylomorphismを再帰呼び出しに書き戻すマクロを作ってみました. 以下がそのマクロの実装. hylol-formで展開後の構造を定義し, hylol-expandで記法のパターンマッチを行い, それに対応するマクロフォームに展開します. 素朴な実装なので, 特に柔軟性などはありません.
(use '[clojure.core.match :only (match)])

(defn hylol-form
  "returns quoted hylol form"
  ([f-in-reduce init-in-reduce pred use-map-option f-in-map
   f-in-iterate init-in-iterate]
     `(letfn [(~'f [~'b]
                (if (~pred ~'b)
                  (~f-in-reduce
                   ~(if use-map-option (list f-in-map 'b) 'b)
                   (~'f (~f-in-iterate ~'b)))
                  ~init-in-reduce))]
        (~'f ~init-in-iterate))))

(defmacro hylol-expand
  "hylomorphism of lists without generating list"
  [hylol-exp]
  (match [hylol-exp]
    [(['reduce f-in-reduce init-in-reduce ;; no map & natural
       (['map f-in-map
         (['take-while pred
           (['iterate f-in-iterate init-in-iterate] :seq)]
            :seq)] :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred true f-in-map
                f-in-iterate init-in-iterate)
    [(['reduce f-in-reduce init-in-reduce ;; no map & natural
       (['take-while pred
         (['iterate f-in-iterate init-in-iterate] :seq)]
            :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred false nil
                f-in-iterate init-in-iterate)
    [(['->> (['iterate f-in-iterate init-in-iterate] :seq) ;; map & thread
       (['take-while pred] :seq)
       (['map f-in-map] :seq)
       (['reduce f-in-reduce init-in-reduce] :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred true f-in-map
                f-in-iterate init-in-iterate)
    [(['->> (['iterate f-in-iterate init-in-iterate] :seq) ;; no map & thread
       (['take-while pred] :seq)
       (['reduce f-in-reduce init-in-reduce] :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred false nil
                f-in-iterate init-in-iterate)
     [the-other-cases]
     (throw (Exception. (str "malformed hylol : " the-other-cases)))))
 上記のようなマクロ定義で, 次のような4通りの記述が可能になります.
(println (hylol-expand
          (reduce + 0
                  (map dec
                       (take-while (partial not= 0)
                                   (iterate dec 100))))))

(println (hylol-expand
          (reduce + 0
                  (take-while (partial not= 0) (iterate dec 100)))))

(println (hylol-expand
          (->> (iterate dec 100)
               (take-while (partial not= 0))
               (map inc)
               (reduce + 0))))

(println (hylol-expand
          (->> (iterate dec 100)
               (take-while (partial not= 0))
               (reduce + 0))))
 reduce, map,  take-while, iterateの位置は固定で, これ以外の書き方はできません. 引数となる「関数/(繰り返しの)初期値」のみが変更可能です. 上の二例は, 引数渡しで普通に書けるケース, 下の二例は, スレッディングマクロで書かれたものを再帰に書きなおすものです.

 簡単なベンチマークテストのため, 前の記事で作成した8queen問題でテストしてみました. 前回書いたのは, 次のような書き方(関数コール)でした.
(defn correct-answer?-v2 [arrangement]
  (hylol andf true #(not (first %)) not-conflict? rest arrangement))
 hylolの変換マクロの場合, 同じ意味のプログラムが,
(defn correct-answer?-mbt [arrangement]
  (hylol-expand ;; = identify
   (->> (iterate rest arrangement)
        (take-while first)
        (map not-conflict?)
        (reduce andf true))))
で書けます. 見た目は, 大分改善されました.

 以下が簡単なベンチマーク結果. correct-answers-with-hylol-macroがhylolの変換のマクロ実装, correct-answers-with-hylol-functionは, 前に書いた記事のcorrect-answers-v2, correct-answersは, ナイーブな(hylolの変換が入らない)実装でこれも前回書いた記事と同じ.
user> (time (count correct-answers-with-hylol-macro))
"Elapsed time: 590.56785 msecs"
92
user> (time (count correct-answers-with-hylol-function))
"Elapsed time: 594.688247 msecs"
92
user> (time (count correct-answers))
"Elapsed time: 739.69738 msecs"
92
 これで, 可読性を維持したまま, Hylomorphismのマクロ展開により高速化が可能になりました.

2015/02/07

Clojure, incanterを使ってクラスタリング, Single Linkage(単連結法) / Complete Linkage(完全連結法)

 Clojureのincanter(統計解析用のライブラリ)が使えるようになったので, とりあえず, クラスタリングをやってみました. Single Linkage (単連結法)および, Complete Linkage(完全連結法)は, 階層的クラスタリングの手法の一種で, 距離空間に散らばった多数の点をその位置からそれぞれグループ(クラスタ)ごとに分類する手法です.

 例えば, 下の図は, 5つのガウス分布に従ってプロットされた二次元上の点です (各ガウス分布ごとに色で分類されています).

 これをComplete Linkageによって, 5つのクラスタへ分類すると, 次のようになります.
 元のガウス分布にいくらか従いながら, 近い(近そうな)点どうしが, それぞれのクラスタとして, ひとまとめにされていることが観察できると思います.

 クラスタリングでは, 空間上の点の位置を元にして同じようなどうしを一つにまとめ, 元の分布を予測したり, グループごとに分類します.

データクラスタの可視化

 incanterでは, 二次元空間上のグラフのプロットや行列の処理が簡単に使えるので, このようなクラスタリングが, 簡単に実装できます.

(require '[incanter.core :as core]
         '[incanter.stats :as st]
         '[incanter.charts :as ch]
         '[clojure.math.numeric-tower :as tower])

;; ------------------------- ------------------------- -------------------------
;; show clusters
;; ------------------------- ------------------------- -------------------------

(defn plot-with-clusters [clusters] ;; clusters = [[[x1 y1] [x2 y2] ...] ...]
  (let [head (first clusters)
        init-plot (ch/scatter-plot (map first head) (map second head))]
    (reduce #(ch/add-points %1 (map first %2) (map second %2))
            init-plot (rest clusters))))
 名前の衝突を避けるためにrequireでincanterをロードします. incanterではscatter-plotで点を二次元座標上にプロットできます(一番目のクラスタをscatter-plotでプロット). add-pointsでchartオブジェクトに点を付加, reduceで各クラスタを一つにまとめます.

サンプルとなるデータ点のリスト

 これに, ランダムなパラメータを持つ複数のガウス分布に従うデータセットを作ります.
(defn generate-random-parameter []
  {:sd (+ 0.15 (rand 0.1)) :x (dec (rand 2)) :y (dec (rand 2))
   :size (+ 20 (rand-int 10))})

(defn generate-sample-set [parts]
  (let [generate-by #(st/sample-normal (:size %1) :mean (%2 %1) :sd (:sd %1))]
    (map (fn [p] (map list (generate-by p :x) (generate-by p :y)))
         (take parts (repeatedly generate-random-parameter)))))
 sample-normalは, ガウス分布に従うサンプルデータを生成します. 標準偏差等のパラメータの範囲は適当です. generate-sample-setでサンプルのデータ点(x, y)のリストを生成します.

 データ的には, [ [[0.1 0.2] [-0.3 0.1] ......] [[0.3 -0.02] ......] ......]のように, 各クラスタごとに(x, y)座標がひとつのリストに.

 (coreの)view関数で, (scatter)プロットされたオブジェクトを表示できます.
user> (core/view (plot-with-clusters (generate-sample-set 3)))
 このリストを(apply concat ls)で, 結合してサンプルのデータ点のリストを作ります. 何か面白いデータがあれば, それでクラスタリングができるのですが...... このデータ列が本来クラスタリングにより分析されるデータになります.

距離

 距離空間には, マンハッタン距離やマハラノビス距離などありますが, 今回は典型的なユークリッド距離を使います. 簡単にいえば, 視覚的な距離がそのまま, データ間の距離になります.
(defn distance [x1 y1 x2 y2]
  (double (tower/sqrt (+ (core/sq (- x1 x2)) (core/sq (- y1 y2))))))

(defn compare-pairs [xy1 xy2]
  (distance (first xy1) (second xy1) (first xy2) (second xy2)))

Single Linkage, Complete Linkageと階層的クラスタリング

 Single/Complete Linkageは, 階層的クラスタリングの一種です. 階層的クラスタリングとは, クラスタ群の中から2つクラスタを選びその2つを結合していき, 目標となるクラスタ数に達するまで結合を続けます. クラスタ群の初期状態では, データ点一つに対し, 一つクラスタを定義します. 階層的クラスタリングの欠点の一つは, 最終的なクラスタ数を指示する必要があることです. 最終的なクラスタ数の適切な値がわからない時にはこの設定が問題となります. しかし, 裏を返せば, 好きなようにクラスタ数を決められるということでもあります.

 Single LinkageとComplete Linkageの違いは, クラスタ間の距離の定義です. Single Linkageでは, 2つのクラスタ間どうしで最も近い要素間の距離がそのクラスタ間の距離になります. 一方, Complete Linkageは, 最も遠い要素間の距離がクラスタ間の距離とします.
(defn single-linkage "cluster a & cluster b -> distance" [a b]
  (reduce min (map #(reduce min (map (partial compare-pairs %) b)) a)))

(defn complete-linkage "cluster a & cluster b -> distance" [a b]
  (reduce max (map #(reduce max (map (partial compare-pairs %) b)) a)))
 プログラム上の違いはmin/maxだけですが, Single Linkageの考え方としては, 2つのクラスタが同じクラスタなら要素間に繋がりがあるはずだということで, Complete Linkageの方は, 同じクラスタに含まれるということは要素がそれなりに凝集しているはずだ, ということでしょう.

IncanterのMatricesと距離行列

 incanterを使っていて, せっかくなので距離行列(クラスタ間の距離のデータを保存している行列)もincanterの行列を使いました.

 Clojureのシーケンス(リスト)から, matrixオブジェクトを生成できます. 同じクラスタ群の距離を保存するための行列なので, 対象行列になり, 対格要素は0になります. ここでの計算において, 対格要素は使いませんが, わかりにくいので対格要素は無限大としました.
 こんな感じ.
user> (generate-init-matrix (map (fn [x] [x]) (apply concat (generate-sample-set 5))))
 A 118x118 matrix
 -----------------
 Infinity  1.79e-01  3.13e-01  .  1.96e+00  2.09e+00  2.26e+00 
 1.79e-01  Infinity  1.78e-01  .  1.85e+00  1.99e+00  2.14e+00 
 3.13e-01  1.78e-01  Infinity  .  1.93e+00  2.08e+00  2.21e+00 
 ... 
 1.96e+00  1.85e+00  1.93e+00  .  Infinity  2.07e-01  3.39e-01 
 2.09e+00  1.99e+00  2.08e+00  .  2.07e-01  Infinity  3.69e-01 
 2.26e+00  2.14e+00  2.21e+00  .  3.39e-01  3.69e-01  Infinity
 matrixオブジェクトは, デフォルトでtoStringメソッドが実装されているため, replでは, 間を省略した行列が表示されます. 上の例は, 対格要素をPOSITIVE_INFINITYにした118×118の行列. 巨大な行列でも間引かれているので表示されるときはコンパクトになります.

(def inf Double/POSITIVE_INFINITY)

(defn distances-xy-xys [xy xys]
  (cons inf (map (fn [xy1] (compare-pairs (first xy) (first xy1))) (rest xys))))

(defn generate-init-matrix [xys]
  (let [dist-vec (map distances-xy-xys (drop-last xys) (iterate rest xys))]
    (core/symmetric-matrix
     (apply concat (concat dist-vec [[inf]])) :lower false)))

(defn find-conj-pair [dist-matrix width]
  (loop [x (dec width) y (dec width) conj-pair [width width] min-value inf]
    (cond
     (< y 0) conj-pair
     (< x 0) (recur y (dec y) conj-pair min-value)
     :else (let [value (core/sel dist-matrix :cols x :rows y)]
             (if (< value min-value)
               (recur (dec x) y [x y] value) ;; update pair
               (recur (dec x) y conj-pair min-value)))))) ;; not update pair

(defn update-dists
  "updates distance matrix
   dist-matrix : distance matrix
   i, j : delete target
   new-column : newly generated distances from i,j"
  [dist-matrix i j new-column]
  (->> (core/sel dist-matrix :except-rows [i j] :except-cols [i j])
       (core/bind-rows new-column)
       (core/bind-columns (cons inf new-column))))
 上2つの関数は, 初期状態の距離行列を作るためのものです. find-conj-pairで距離行列から結合する2つのクラスタのペアをインデックスで返します. update-distsで距離行列を更新します. 新しく生成されたクラスタと他のクラスタ間の距離を行列の上の行(と左側の列)に付け加えて更新します.

階層的クラスタリング

 クラスタ群は, 簡単のため, i,jのインデックスで扱います.
(defn remove-clusters [i j clusters]
  (let [;; supposed to be c_0 c_1 ... c_n ... c_m ... c_max
        n (min i j) m (max i j)
        ;; [[c_0  ... c_n ... c_m-1] [c_m ... c_max]]
        first-m-last (split-at m clusters)
        ;; [[c_0 ... c_n-1] [c_n ... c_m-1]]
        first-n-m (split-at n (first first-m-last))]
    ;; (concat [c_0 ... c_n-1] [c_n+1 ... c_m-1] [c_m+1 ... c_max])
    (concat (first first-n-m) (rest (second first-n-m)) (rest (second first-m-last)))))
 remove-clustersでクラスタ群のi,j番目のクラスタを除去します.

 以下は, クラスタリングのループです. 一番近いクラスタを探し出し, 結合して更新という処理を指定されたクラスタ数になるまで続けます.
(defn clustering [calc-dists xys target-size]
  (loop [clusters (map (fn [x] [x]) xys)
         dist-matrix (generate-init-matrix clusters)
         cluster-size (count xys)]
    (if (<= cluster-size target-size)
      clusters
      (let [conj-pair (find-conj-pair dist-matrix cluster-size)
            i (first conj-pair) j (second conj-pair)
            new-cluster (concat (nth clusters i) (nth clusters j))
            other-clusters (remove-clusters i j clusters)]
       (recur (cons new-cluster other-clusters)
              (update-dists dist-matrix i j (calc-dists new-cluster other-clusters))
              (dec cluster-size))))))

(defn link-with [dist-f new-cluster other-clusters]
  (map (fn [cluster] (dist-f new-cluster cluster)) other-clusters))

(defn clustering-with-single-linkage [xys number-of-clusters]
  (clustering (partial link-with single-linkage) xys number-of-clusters))

(defn clustering-with-complete-linkage [xys number-of-clusters]
  (clustering (partial link-with complete-linkage) xys number-of-clusters))

実行結果

こんな感じでサンプルのデータ点を生成して,
user> (def A (generate-sample-set 5))
user> (core/view (plot-with-clusters A))

こんな感じのガウス分布になっているものに, 最終的なクラスタ数を5に設定してSingle Linkageをつかうと,
user> (core/view (plot-with-clusters (clustering-with-single-linkage (apply concat A) 5)))

となり, 同様に最終的なクラスタ数を5に設定したComplete Linkageでは,
user> (core/view (plot-with-clusters (clustering-with-complete-linkage (apply concat A) 5)))
になります.

 Complete Linkageの方が, うまく(?)クラスタリングされているように見えます. 元の分布を完全に予測しているわけではありませんが, 生成されたときの分布4に近いかたちで分類できています. Single Linkageは, 孤立した点が一つのクラスタに分類されてしまう傾向があるようです. 全体がほぼ2つのクラスタに分類されてしまっています.

参考文献
 全ソースコード.
(require '[incanter.core :as core]
         '[incanter.stats :as st]
         '[incanter.charts :as ch]
         '[clojure.math.numeric-tower :as tower])

;; ------------------------- ------------------------- -------------------------
;; show clusters
;; ------------------------- ------------------------- -------------------------

(defn plot-with-clusters [clusters] ;; clusters = [[[x1 y1] [x2 y2] ...] ...]
  (let [head (first clusters)
        init-plot (ch/scatter-plot (map first head) (map second head))]
    (reduce #(ch/add-points %1 (map first %2) (map second %2))
            init-plot (rest clusters))))

;; ------------------------- ------------------------- -------------------------
;; generate data set
;; ------------------------- ------------------------- -------------------------

(defn generate-random-parameter []
  {:sd (+ 0.15 (rand 0.1)) :x (dec (rand 2)) :y (dec (rand 2))
   :size (+ 20 (rand-int 10))})

(defn generate-sample-set [parts]
  (let [generate-by #(st/sample-normal (:size %1) :mean (%2 %1) :sd (:sd %1))]
    (map (fn [p] (map list (generate-by p :x) (generate-by p :y)))
         (take parts (repeatedly generate-random-parameter)))))

;; ------------------------- ------------------------- -------------------------
;; clustering
;; ------------------------- ------------------------- -------------------------

(defn distance [x1 y1 x2 y2]
  (double (tower/sqrt (+ (core/sq (- x1 x2)) (core/sq (- y1 y2))))))

(defn compare-pairs [xy1 xy2]
  (distance (first xy1) (second xy1) (first xy2) (second xy2)))

(defn single-linkage "cluster a & cluster b -> distance" [a b]
  (reduce min (map #(reduce min (map (partial compare-pairs %) b)) a)))

(defn complete-linkage "cluster a & cluster b -> distance" [a b]
  (reduce max (map #(reduce max (map (partial compare-pairs %) b)) a)))

(def inf Double/POSITIVE_INFINITY)

(defn distances-xy-xys [xy xys]
  (cons inf (map (fn [xy1] (compare-pairs (first xy) (first xy1))) (rest xys))))

(defn generate-init-matrix [xys]
  (let [dist-vec (map distances-xy-xys (drop-last xys) (iterate rest xys))]
    (core/symmetric-matrix
     (apply concat (concat dist-vec [[inf]])) :lower false)))

(defn find-conj-pair [dist-matrix width]
  (loop [x (dec width) y (dec width) conj-pair [width width] min-value inf]
    (cond
     (< y 0) conj-pair
     (< x 0) (recur y (dec y) conj-pair min-value)
     :else (let [value (core/sel dist-matrix :cols x :rows y)]
             (if (< value min-value)
               (recur (dec x) y [x y] value) ;; update pair
               (recur (dec x) y conj-pair min-value)))))) ;; not update pair

(defn update-dists
  "updates distance matrix
   dist-matrix : distance matrix
   i, j : delete target
   new-column : newly generated distances from i,j"
  [dist-matrix i j new-column]
  (->> (core/sel dist-matrix :except-rows [i j] :except-cols [i j])
       (core/bind-rows new-column)
       (core/bind-columns (cons inf new-column))))

(defn remove-clusters [i j clusters]
  (let [;; supposed to be c_0 c_1 ... c_n ... c_m ... c_max
        n (min i j) m (max i j)
        ;; [[c_0  ... c_n ... c_m-1] [c_m ... c_max]]
        first-m-last (split-at m clusters)
        ;; [[c_0 ... c_n-1] [c_n ... c_m-1]]
        first-n-m (split-at n (first first-m-last))]
    ;; (concat [c_0 ... c_n-1] [c_n+1 ... c_m-1] [c_m+1 ... c_max])
    (concat (first first-n-m) (rest (second first-n-m)) (rest (second first-m-last)))))

(defn clustering [calc-dists xys target-size]
  (loop [clusters (map (fn [x] [x]) xys)
         dist-matrix (generate-init-matrix clusters)
         cluster-size (count xys)]
    (if (<= cluster-size target-size)
      clusters
      (let [conj-pair (find-conj-pair dist-matrix cluster-size)
            i (first conj-pair) j (second conj-pair)
            new-cluster (concat (nth clusters i) (nth clusters j))
            other-clusters (remove-clusters i j clusters)]
        (recur (cons new-cluster other-clusters)
               (update-dists dist-matrix i j (calc-dists new-cluster other-clusters))
               (dec cluster-size))))))

(defn link-with [dist-f new-cluster other-clusters]
  (map (fn [cluster] (dist-f new-cluster cluster)) other-clusters))

(defn clustering-with-single-linkage [xys number-of-clusters]
  (clustering (partial link-with single-linkage) xys number-of-clusters))

(defn clustering-with-complete-linkage [xys number-of-clusters]
  (clustering (partial link-with complete-linkage) xys number-of-clusters))

2015/02/03

incanterでhistgramをプロットする

 Clojureのincanterについては, ずいぶん前から知っていたのですが, 使う機会があまりなかった為, 特に触っていませんでした. 最近, Clojureで生成したデータからグラフをプロットするという作業が必要になったこともあって, incanterを弄っているので, メモです.


 とりあえず, histgramを使うところまで.

 以下のコマンドは, REPLに
user> (use '(incanter core stats charts pdf))
を叩き込んだ前提で. (2015/02/05, pdf追加)
 名前的に他のライブラリ等の関数名と名前が重複する可能性がある(特にnumeric-towerなど, な)ので, prefixなしで直にロードするというのはあまりよくありませんが, 弄ってみるだけなら, これでも問題ない筈です.

 公式のチュートリアル(Github)にも載っていますが, ヒストグラムの表示は,
user> (view (histogram (sample-normal 10000)))
でいけます.

 ヒストグラムは, n個のデータサンプルをカテゴリ別に分類し, そのデータの分布を可視化するものなので, カテゴリ数を調整したくなります. カテゴリ数は, :nbinsオプションで指定できます. カテゴリ数を1000にする場合, 以下のように書きます.
user> (view (histogram (sample-normal 10000) :nbins 1000))
これを実行すると, 次のようになります.

 (sample-normal 10000)は, 意味的には,
user> (import 'java.util.Random)
user> (def r (Random .))
user> (map (fn [x] (.nextGaussian r)) (range 10000))
と一緒のようです. 平均0, 標準偏差1の正規分布に従い乱数を生成するものです.

 sample-normalが生成するのは, 単なる数値型(double)のリストです.
user> (sample-normal 3)
(1.2821074070154617 1.638783357428357 -0.2801943679292413)
 サンプルのリストをconcatでつなげれば, 別々のサンプルを合成できます. :sdは標準偏差(Standard Deviation), :meanは平均のオプション.
user> (view (histogram (concat (sample-normal 10000 :mean -5 :sd 4)
                               (sample-normal 10000 :mean 4)) :nbins 1000))

 ヒストグラムが, 正規分布に従っているかどうかを調べるために, 正規分布の曲線を重ねるのは, histgramで生成したオブジェクトに, add-line関数で重ねます. pdf-normalは, 正規分布の関数. histgramの:densityオプションは, 密度(Density)か頻度(Frequency)のオプション.
user> (view (doto (histogram (sample-normal 10000) :nbins 100 :density true)
                  (add-lines (range -4 4 0.01) (pdf-normal (range -4 4 0.01)))))
で, これが次のようにプロットされます.