x-means
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.19.3377
クラスタリングの方法は、色々あるけども、わたしの知る限り、ほとんどの方法では、事前にいくつのクラスターに分かれるか決める必要がある。クラスタリングしたい多くの状況で、いくつのクラスターに分かれるか事前には分からないというシチュエーションは山ほどあると思う。
ネットワークのクラスタリングには、Newman法とかあるらしいけど、x-meansはk-meansの最適なkを自動で決定するアルゴリズムであるらしい。やってることは単純で、k=2でk-meansを使って、どんどん分割していって、適切なとこまで分割できたら、そこで分割を打ち切る。こんなことは、最適なクラスター数が分からんという事態に直面したら誰でも考えるので、適切かどうかの判定をどうやってやるかが肝。論文ではBICないしAICというのを使えばいいよと書いてある。BICが何かは知らないけれど、対数尤度に謎の項を足したもの。多分大雑把には、k-meansは各クラスターがクラスター中心の周りに正規分布していると仮定しているようなものなので、このモデルがどれくらい適合しているかを測る指標があればよいという感じなんじゃないかと思う。まあ、BICが何か分からなくても、論文に具体的に式が書いてあるので実装はできる
OpenCVのサンプルに、k-meansを使った減色/画像分割があるので、これをx-meansで置き換えてみる
http://opencv.jp/sample/misc.html#color-sub
実際にやってみると、クラスター数が数百になったりして(画像サイズにもよるだろうけど)、相当遅い。ていうか、実装は合ってるのか怪しい。まあ、結論としては、意味なかった。思うに、共分散が単位行列に比例するとしているので、特定方向のみ分散が大きいようなケースでは分割しすぎるという問題があるんじゃなかろうか。あと、情報量基準でよいモデルを選択できるという根拠がいまいちよくわかってない。もっとroughに適合度検定とかじゃダメなんだろうか
実装。OpenCVの新しいPythonインターフェースは分からないことが多いので、とりあえず動くことだけを考えて作ってる
import sys sys.path.append("C:\\\\OpenCV2.1\\Python2.6\\Lib\\site-packages") from cv import * import numpy from numpy.linalg import norm from math import log from struct import * def xmeans(vectors , N0=1): R = len(vectors) M = len(vectors[0]) def IC(indices , K): eps = 1.0e-6 Rn = float(len(indices)) center= reduce(lambda x,y:x+y , [vectors[idx] for idx in indices] , numpy.zeros(M))/Rn variance = eps + sum([norm(vectors[idx]-center)**2 for idx in indices])/Rn #AIC #ic = Rn*log(2*numpy.pi) + Rn*M*log(variance) + (Rn-K) - 2*Rn*log(Rn/R) + 2*(M+1)*K #BIC ic = Rn*log(2*numpy.pi) + Rn*M*log(variance) + (Rn-K) - 2*Rn*log(Rn/R) + K*(M+1)*log(Rn) return ic def kmeans(indices , N): clusters = CreateMat(len(indices) , 1 , CV_32SC1) points = CreateMat(len(indices) , M , CV_32FC1) for i in xrange(points.rows): for j in xrange(points.cols): idx = indices[i] points[i,j] = vectors[idx][j] KMeans2(points , N , clusters , (CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 10, 1.0)) ret = [[] for _ in xrange(N)] for i in xrange(clusters.rows): idx = int(clusters[i,0]) ret[idx].append(indices[i]) return ret def split(indices): result = [] stack = [indices] while True: if len(stack)==0:return result targets = stack[-1] if len(targets)<=1: stack.remove(targets) result.append(targets) continue ic1 = IC(targets , 1) D1,D2 = kmeans(targets,2) if len(D1)==0 or len(D2)==0: stack.remove(targets) result.append(targets) continue ic2 = IC(D1,2)+IC(D2,2) if ic2 < ic1: stack.remove(targets) stack.append(D1) stack.append(D2) else: stack.remove(targets) result.append(targets) result = [] ls = kmeans(range(R) , N0) for indices in ls: result += split(indices) ret = numpy.zeros(R , numpy.uint32) for i,indices in enumerate(result): for idx in indices: ret[idx] = i return ret,len(result) if __name__=="__main__": src_img = LoadImage("test.jpg") sz = src_img.width*src_img.height dst_img = CloneImage(src_img) vectors = [] imgstr = src_img.tostring() for i in xrange(sz): c1,c2,c3 = unpack('BBB' , imgstr[i*3:i*3+3]) vectors.append( numpy.array([float(c1),float(c2),float(c3)]) ) print "start x-means" clusters,k = xmeans(vectors,2) print "%d clusters" % k count = numpy.zeros(k) r = numpy.zeros(k) g = numpy.zeros(k) b = numpy.zeros(k) for i in xrange(sz): idx = clusters[i] count[idx] += 1 (cb,cg,cr) = tuple(vectors[i]) r[idx] += cr g[idx] += cg b[idx] += cb for idx in range(k): if count[idx]>0.0: r[idx] /= count[idx] g[idx] /= count[idx] b[idx] /= count[idx] for h in xrange(dst_img.height): for w in xrange(dst_img.width): idx = clusters[h*dst_img.width+w] dst_img[h,w] = (b[idx] , g[idx] , r[idx]) NamedWindow ("result", CV_WINDOW_AUTOSIZE) ShowImage ("result", dst_img) WaitKey (0) DestroyWindow("result")
Zipfの法則を真面目に検定してみる
Zipfの法則は、「ユリシーズ」に現れる単語の出現頻度を多い順に並べると、出現頻度は概ね順位に反比例することを発見したのが始まりらしい。Zipfは“Human Behavior and the Principle of Least Effort”という本を書いて、都市の人口とかにも当てはまるとか言ったらしいけど、本は読んでないので、実際何が書いてあるかは知らない。Zipfの法則は、色々な現象で普遍的に観察される現象だということになってるらしい。
当たり前の話だけど、単調減少する関数fで、x->∞でf(x)->0になるようなものを考えると、fがある程度緩やかに減衰するなら、漸近的にf(x)〜c*x^aと書けて、両辺のlogとれば、log f= log(c)+a*log(x)なので、両対数グラフを最小二乗法でfittingすれば、とりあえず、指数aと係数cは求まってしまう。偶然aが-1付近になってしまうことも多々あるだろう。そんなわけで、きちんと検定してないものはZipf則に従ってると言えるか怪しい。
Zipfの法則を満たすと言われてるものは、色々あるので、1個ずつ見ていく
1)Webページへのアクセス頻度
これを最初に言い出したのが誰かは分からない。沢山の人が独立に同じことを言ってるのかもしれない。一例として
http://ci.nii.ac.jp/naid/110002937400/
を発見したけれど、単純にグラフにプロットして、大体成り立ってるよねと言ってるだけだった
2)苗字の頻度分布
Power-law Distribution of Family Names in Japanese Societies
http://arxiv.org/abs/cond-mat/9912035
日本人の名字の統計解析
http://ci.nii.ac.jp/naid/110003502725
どちらも適当にfittingしているだけの論文
3)都市の人口
これはZipf自身が既に言ってたことであるらしい。ので、関連研究も多く、全部を調べるとか流石にできない。ただ例えば、日本の場合、東京＞横浜＞大阪の順らしいけど、東京は、23区で測るらしい。それってどうなんだろうか。なんとなく作為的な匂いを感じなくもない。そもそも、都市といっても、RPGじゃないんだから町とフィールドがくっきり分かれてないわけで、どこを境界とするかは作為的にならざるをえない気がするけど。
4)遺伝子発現量
General Statistics of Stochastic Process of Gene Expression in Eukaryotic Cells
http://www.genetics.org/cgi/content/full/161/3/1321
Universality and flexibility in gene expression from bacteria to human
http://www.pnas.org/content/101/11/3765.abstract
5)地震の規模/粒径分布
地震の規模についてはGutenberg-Richter則と呼ばれ、粉体の粒径分布はRosin-Rammler則として、Zipfとかとは無関係に、それぞれの分野で経験則として確立しているらしい。Wikipediaには、これらもZipfの法則に入れられているが、従属変数がrankではなく、マグニチュードや粒子半径のような連続量なので、これをZipfの法則というのは、話を広げすぎだと思う。
というわけで、まともに検定もせず、適当にグラフ書いただけでZipf則とか言ってる人が沢山いるらしいことは分かった。でも、それだけではZipf則を否定する理由にもならないので、アメリカ/日本の都市人口と単語の出現頻度についてカイ二乗検定をしてみた。期待頻度は、rankをnとしたときにX/nとなると仮定している。Xをどのように決定するかは問題だけど、単純に最小二乗推定で決めた(というか、多分カイ二乗検定では、期待頻度の和と観測頻度の和は一致するのが原則なんじゃないかと思うので、それしかない)
単語の出現頻度は、
http://www.cs.cmu.edu/~cburch/words/top.html
から持ってきたデータ。都市の人口分布は、適当にぐぐって拾った。
結果としては、いずれもZipf則を満たすという仮説は棄却された(有意水準は5%)。そもそもカイ二乗検定で検定するのが適切かは疑問がなくもないが、実際のカイ二乗値を見ると誤差が相当大きい。なのでZipf則は、何か背景に深い理由があるというよりは、冒頭に書いたような理由によって、偶然そう見えるだけの産物と考えるのが妥当なのではないかと思う。特に、両対数グラフにプロットすると、誤差が見えにくくなるというのは一役買ってるように思う。傾きが-1前後あたりになりやすいというのも、傾きが-2とか-3になってくると、寡占・独占状態に近づき、逆に-0.5とか-0.4になると相当flatな感じに近づいて、いずれにせよ、-1から大きく離れると、人間が集計してみようという気になりづらいというのが一因ではないかと思う。
まあ、Zipfの法則も大雑把な経験則としては役に立つこともあるかもしれない
検証に使ったコード
from math import * def normsdist(u): d1=0.0498673470 d2=0.0211410061 d3=0.0032776263 d4=0.0000380036 d5=0.0000488906 d6=0.0000053830 if u > 0: temp=1+d1*u+d2*pow(u,2)+d3*pow(u,3)+d4*pow(u,4)+d5*pow(u,5)+d6*pow(u,6) return pow(temp,-16.0)/2 else: u=-u temp=1+d1*u+d2*pow(u,2)+d3*pow(u,3)+d4*pow(u,4)+d5*pow(u,5)+d6*pow(u,6) return 1.0 - pow(temp,-16.0)/2 def chidist(chi2 , f): if chi2 < 0.0: chi2 = 0.0 if f > 40: temp = (pow(chi2/float(f),1.0/3) - (1-2.0/9.0/float(f)))/sqrt(2.0/9.0/float(f)) return normsdist(temp) elif f%2==0: s = t = exp(-0.5 * chi2) for k in range(2, f, 2): t *= chi2 / k s += t return s else: chi = sqrt(chi2) if (f == 1): return 2 * normsdist(chi); s = t = chi * exp(-0.5 * chi2) / sqrt(2 * pi) for k in range(3, f, 2): t *= chi2 / k s += t return 2 * (normsdist(chi) + s) def zipftest(dat , alpha=0.05): ls = sorted(dat,reverse=True) X = sum(ls)/sum([1.0/n for n in xrange(1,len(ls)+1)]) chi2 = sum([(x - X/(n+1))**2/(X/(n+1)) for (n,x) in enumerate(ls)]) return chidist(chi2 , len(dat)-1)>alpha def main(): print "test:",zipftest([100.0/n for n in range(1,10)]) print "Americaの都市人口:",zipftest([8143197,3845541,2842518,2016582,1463281,1461575,1256509,1255540,1213825,912332]) print "日本の都市人口:",zipftest([8396594,3559867,2633819,2204496,1870170,1521362,1463941,1393659,1306,992,1146413]) if __name__=="__main__": main()
CiNiiの著者名検索の精度検証
http://cinii.jp/post/486298233/cinii-author-search
どんなもんかと。日本人研究者を対象に、なるべく分野が偏らないように集計してみた。集計自体は完全に正確ではないだろうけど、誤差数%程度なんじゃないかと思う。適合率/再現率は、該当する著者の文献が最も多いクラスターで測定した
北沢宏一 適合率100%再現率67.3%(312/463) 佐藤勝彦@宇宙論 適合率91.9%(57/62)再現率27.5%(57/207) 萩谷昌己 適合率100%再現率73.6%(95/129) 平木敬 適合率100%再現率79.2%(175/221) 木下佳樹 適合率100%再現率31.5%(6/19) 米澤明憲 適合率100%再現率73.8%(93/126) 林崎良英 適合率100%再現率33.9%(80/236) 横山茂之 適合率100%再現率27.3%(64/234) 西郷薫 適合率100%再現率33.3%(18/54) 牧野淳一郎 適合率100%再現率18.4%(12/65) 治部眞里 適合率100%再現率33.3%(9/27) 北浦和夫@FMO 適合率100%再現率15.2%(5/33) 柏原正樹 適合率100%再現率39.7%(25/63) 広田良吾 適合率100%再現率51.2%(87/170) 平均 適合率(略)再現率50.7%(1038/2047)
非日本人研究者はそもそも少ないのだろうけどと、Fedorov Dmitri(FMOとか研究してる人)とかは、全く名寄せされてなかった。
全体としては、適合率〜100%、再現率〜50%という感じだと思う。理情のセンセーだけ妙に再現率がいいので、そのおかげで、若干底上げされてる感もなきにしもあらず。平均的にはやる気なさすぎだし、わたしなら適合率〜100%、再現率〜90%くらいまでは確実にいけるので、NIIはお金出してくれんかなぁ
