安定分布に従う乱数の生成

安定分布を生成させたくなったので、やり方を調べた。Chambers-Mallows-Stuckが提案した方法を使うのが一般的らしい(scipy,GLSLはこの方法に基づいた実装のようだった)。一様分布Uと指数分布E(平均が1となるように。これは手軽には、区間(0,1]上の一様分布Vに対して、-log(V)を計算すればいい)を組み合わせて、γ=1,δ=0の標準安定分布X(α,β)を作って、最後にスケールγと位置δの調整をする。基本的には、逆関数法であるけど、計算は結構大変。まぁ、何も分かってなくても、式を丸写しすれば実装出来てしまうけど...

#!/usr/bin/env python
import math,random

"""
J.M. Chambers, C.L. Mallows and B.W. Stuck (1976) 
"A Method for Simulating Stable Random Variables"  
JASA, Vol. 71, No. 354. pages 340-344  

"""
def random_stable(alpha,beta,gamma=1.0,delta=0.0):
   assert(alpha>0 and alpha<=2),("invalid argument:alpha=%f" % alpha)
   assert(beta>=-1 and beta<1),("invalid argument:beta=%f" % beta)
   assert(gamma > 0),("invalid argument:gamma=%f" % gamma)
   #-- 理屈上は、generic caseだけあれば十分
   if alpha==2:                         #-- Gaussian distribution
      X = random.gauss(0,math.sqrt(2))
   elif alpha==1.0 and beta==0.0:       #-- Cauchy distribution
      X = math.tan( random.uniform(-math.pi/2 , math.pi/2) )
   elif alpha==0.5 and fabs(beta)==1.0: #-- Levy distribution
      X = beta/(random.gauss(0,1)**2)
   elif alpha!=1.0 and beta==0.0:       #-- Symmetric alpha-stable distribution(except Cauchy case)
      U = random.uniform(-math.pi/2 , math.pi/2)
      E = random.expovariate(1.0)       #-- or -math.log(random.random())
      X = math.sin(alpha*U)/(math.pow(math.cos(U) , 1/alpha)) * math.pow(math.cos(U*(1-alpha))/E,(1-alpha)/alpha)
   elif alpha==1.0:                     #-- generic case for alpha==1
      U = random.uniform(-math.pi/2 , math.pi/2)
      E = random.expovariate(1.0)
      S = math.pi/2 + beta*U
      X = (2/math.pi)*(S*math.tan(U) - beta*math.log((math.pi*E*math.cos(U)/2)/S))
   else:                                #-- generic case for alpha!=1
      U = random.uniform(-math.pi/2 , math.pi/2)
      E = random.expovariate(1.0)
      t = beta*math.tan(math.pi*alpha/2)
      B = math.atan( t )
      S = math.pow(1+t*t , 1/(2*alpha))
      X = S*math.sin(alpha*U+B)/math.pow(math.cos(U),1/alpha) * math.pow(math.cos(U*(1-alpha)-B)/E,(1-alpha)/alpha)
   #-- scale and location
   if alpha!=1.0:
      return (gamma*X + delta)
   else:
      return (gamma*X + (2/math.pi)*beta*gamma*math.log(gamma) + delta)

実装しただけではアレなので、昔から気になってた実験。安定分布は、α<2の時、分散が発散し、α<1の時、平均も発散する(α<2では、αより大きい次数のモーメントは発散する)。のだけど、標本平均、標本分散は当然∞にはならないので、実際にどういう風に動くのか見てみたかった

def test_stbl(alpha,N=10):
   XS = []
   for n in xrange(N):
     XS = XS+[random_stable(alpha , 0.0 , 1.0) for _ in xrange(10000)]
     m = sum(XS)/len(XS)
     s = math.sqrt( sum([(x-m)*(x-m) for x in XS])/len(XS) )
     print (n+1)*10000,m,s

実際に、β=0,γ=1,δ=0で試してみると、
・αが小さいほど、同一データ数に於ける標本分散は大きくなる
・α<2の時、データ数を増やすと、大きい標本分散が出やすくなる
・α>=1なら、データ数が増えるにつれてちゃんと平均は0に収束していく(これは当然だが)
・α<1なら、データ数が増えるにつれて、平均は負の値になったり正の値になったりしつつ、大きい絶対値が出やすくなる
という傾向が見えた。歪度も(3次のモーメントだから)発散するので、実際に、標本歪度を計算すると、標本分散と同様、突然、正負が入れ替わって、収束しない(安定分布のパラメータβも歪み度とか言われるので、紛らわしいけど)

ついでに、β=0,γ=1,δ=0の対称安定分布を、こないだ実装したshapiro-wilk検定に放り込んでみる。

def test_stbl2(alpha,L,N=1000):
   cnt = 0
   for _ in xrange(N):
       xs = [random_stable(alpha , 0.0 , 1.0) for _ in xrange(L)]
       if swtest(xs)>0.05:cnt+=1
   print "%2.2f(%%) passed" % (100*float(cnt)/N)

当たり前だけど、αが2に近いほど、沢山誤判定する。とはいえ、α=1.8で、データ数25くらいでも、テストに通る(P値の閾値は0.05として)のは70%程度なので、まぁ割と頑張ってると言える？

Mandelbrotが綿花の価格変動の分布は、正規分布よりα=1.7の安定分布に近いと言ったらしいけど、α=1.7で、P-valueの閾値を0.05とすると
データ数25:Shapiro-Wilk検定を通る確率は60%程度
データ数50:Shapiro-Wilk検定を通る確率は40%程度
データ数100:Shapiro-Wilk検定を通る確率は17~8%程度
データ数250:Shapir-Wilk検定を通る確率は1%程度
という感じになった。