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()