ピアソンIV型分布の乱数生成 (2017/05/10)

PythonではピアソンIV型分布の乱数を生成するプログラムは私の知る限り、存在しない。そこで、前にピアソンIV型分布の乱数を生成するプログラムをモンテカルロ法で作成したが、あまりにも時間がかかるので結局使わなかった。今度は逆関数法でやってみる。

うまい具合にC言語で書かれたプログラムを掲載している論文があった。R言語のPearsonDSパッケージもこのプログラムを利用しているらしい(どこで見たかは忘れた)。

https://www-cdf.fnal.gov/physics/statistics/notes/cdf6820_pearson4.pdf

そこで今回はこのプログラムをPythonに直訳してみる。このプログラムでは乱数は1個ずつ生成するので遅い。配列で生成するプログラムを後日作成する。

少し腑に落ちないところがいくつかある。rpears4()関数の引数はm=尖度、nu=歪度、a=標準偏差、lam=平均であると思われる。また、正規分布の尖度は3とする場合と0とする場合がある。rpears4()関数の場合は3としているようだ。

すると、rpears4(3, 0, 1, 0)とすれば正規分布に従う乱数が生成するはずである。だが、確認してみると、歪度、平均は大体合っているが、標準偏差は小さく、尖度は大きい傾向がある。

また、グラフで最頻値である0前後で縦軸が0.8を超えるのは大きすぎるように思われる。縦軸は正規化されており、ビンごとの発生確率のはずである。だが、これはmatplotlibの問題だろう。

この辺りをチェックしながら記事を随時修正、追加するつもりである。

①以下のコードを「test.py」ファイルとして「~/py」フォルダーに保存する。

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

DBL_EPSILON = 1.0e-9
M_2_SQRTPI = 2 / np.sqrt(np.pi)
M_PI_2 = np.pi / 2

def gammar2(x, y):
    y2 = y * y  # const
    xmin = 2 * y2 if 2*y2 > 10.0 else 10.0  # const
    r = 1
    s = 1
    p = 1
    f = 0
    count = 0
    while x < xmin:
        if count == 0:
            t =  y / x  # const
            count = 1
        x += 1
        r *= 1 + t*t
    while p > s*DBL_EPSILON:
        p *= y2 + f*f
        f += 1
        p /= x * f
        x += 1
        s += p
    return 1.0 / (r*s)

def type4norm(m, nu, a):
    assert m > 0.5, 'error'
    return (0.5 * M_2_SQRTPI * gammar2(m, 0.5*nu)
            * np.exp(math.lgamma(m)-math.lgamma(m-0.5)) / a)

def ranu():
    return np.random.rand()

def rpears4(m, nu, a, lam):
    k = type4norm(m, nu, 1.0)  # const
    b = 2*m - 2  # const
    M = np.arctan(-nu / b)  # const
    cosM = b / np.sqrt(b*b + nu*nu)  # const
    r = b*np.log(cosM) - nu*M  # const
    rc = np.exp(-r) / k
    x = 0
    z = 0
    assert m > 1, 'error'
    while True:
        s = 0
        z = 0
        x = 4*ranu()  # have to check
        if x > 2:
            x -= 2
            s = 1
        if x > 1:
            z = np.log(x - 1)  # have to check
            x = 1 - z  # have to check
        x = M + rc*x if s else M - rc*x  # have to check
        if np.cos(x) < 0:  # added
            continue
        if not ((math.fabs(x) >= M_PI_2)
              | (z + np.log(ranu()) > b*np.log(np.cos(x)) - nu*x - r)):
            break
    return a*np.tan(x) + lam

if __name__ == '__main__':
    n = 100000
    rnd = np.zeros(n)
    for i in range(n):
        rnd[i] = rpears4(3, 0, 1, 0)
    mean = np.mean(rnd)
    std = np.std(rnd)
    skew = stats.skew(rnd)
    kurt = stats.kurtosis(rnd)
    print('mean = ', mean)
    print('std = ', std)
    print('skew = ', skew)
    print('kurt = ', kurt)
    plt.hist(rnd, bins=50, normed=True, range=(-4, 4))
    plt.savefig('pearson4.png', dpi=150)
    plt.show()
    plt.close()

②以下のコマンドをIPyhonコンソールに入力して「Enter」キーを押す。

%run -t ~/py/test.py


mean =  -0.00187551097649
std =  0.579023696444
skew =  0.12142983501094355
kurt =  6.15981877928405

論文とWikipediaの対応関係

指定した尖度、歪度、標準偏差、平均に従う乱数を生成したい。だが、rpears4()関数の引数であるm、nu、a、lamがどのようなものであるか、いまいちはっきりしない。

そこで、Wikipediaの記事を参考にしながら、はっきりさせようと思う。

論文p1によると、ピアソンIV型分布の確率密度関数は

\[ f(x)dx = k\left[1 + \left(\frac {x - \lambda} {a}\right)^{2} \right]^{-m}\exp\left[-v\tan^{-1}\left(\frac {x - \lambda} {a} \right)\right]dx \]

となっている。

一方、Wikipediaでは以下のようになっている。(https://en.wikipedia.org/wiki/Pearson_distribution)

\[ p(x) = \frac {\left|\frac {\Gamma\left(m + \frac {v} {2}i\right)} {\Gamma(m)}\right|^{2}} {\alpha B\left(m - \frac {1} {2}, \frac {1} {2}\right)}\left[1 + \left(\frac {x - \lambda} {\alpha}\right)^{2} \right]^{-m}\exp\left[-v\tan^{-1}\left(\frac {x - \lambda} {\alpha} \right)\right] \]

論文の式から両辺の\( dx \)を除き、\( f(x) \)を\( p(x) \)に、\( a \)を\( \alpha \)に置き換えると以下の式になる。

\[ p(x) = k\left[1 + \left(\frac {x - \lambda} {\alpha}\right)^{2} \right]^{-m}\exp\left[-v\tan^{-1}\left(\frac {x - \lambda} {\alpha} \right)\right] \]

すると、

\[ k = \frac {\left|\frac {\Gamma\left(m + \frac {v} {2}i\right)} {\Gamma(m)}\right|^{2}} {\alpha B\left(m - \frac {1} {2}, \frac {1} {2}\right)} \]

となる。

論文p5には

\[ k = \frac {\Gamma(m)} {\sqrt{\pi}a\Gamma(m-1/2)} {\left|\frac {\Gamma(m + iv/2)} {\Gamma(m)}\right|^{2}} \]

とある。すると、

\[ \frac {1} {\alpha B\left(m - \frac {1} {2}, \frac {1} {2}\right)} = \frac {\Gamma(m)} {\sqrt{\pi}a\Gamma(m-1/2)} \]

ということになる。

\( a \)と\( \alpha \)は同じだから変形すると、

\[ \frac {1} {B\left(m - \frac {1} {2}, \frac {1} {2}\right)} = \frac {\Gamma(m)} {\sqrt{\pi}\Gamma(m-1/2)} \] \[ B\left(m - \frac {1} {2}, \frac {1} {2}\right) = \frac {\sqrt{\pi}\Gamma(m-1/2)} {\Gamma(m)} \] \[ B\left(m - \frac {1} {2}, \frac {1} {2}\right) = \frac {\Gamma(m-1/2)\sqrt{\pi}} {\Gamma(m)} \]

となる。

ところで、Wikipediaによるとベータ関数とガンマ関数は以下のような関係にある。(https://ja.wikipedia.org/wiki/ベータ関数)

\[ B(x, y) = \frac {\Gamma(x)\Gamma(y)} {\Gamma(x + y)} \]

ここで、

\[ x = m - \frac {1} {2}, \ y = \frac {1} {2} \]

とすると、

\[ \Gamma\left(\frac {1} {2}\right) = \sqrt{\pi} \]

となる。

ガンマ関数と円周率との間にはこのような関係にあることがWikipediaに書かれている。(https://ja.wikipedia.org/wiki/ガンマ関数)

これで論文とWikipediaにおける数式の対応関係が分かった。では本題のm、nu、a、lamについて調べる。

mについて

先ず、\( m \)について調べる。「Pearson distribution」には

\[ m = \frac {1} {2b_{2}} \tag{①} \]

とある。では\( b_{2} \)とは何かというと、

\[ b_{2} = \frac {2\beta_{2} - 3\beta_{1} - 6} {10\beta_{2} - 12\beta_{1} - 18} \tag{②} \]

となっている。では\( \beta_{1} \)、\( \beta_{2} \)は何かというと、

\[ \beta_{1} = \gamma_{1}^{2} \tag{③} \] \[ \beta_{2} = \gamma_{2} + 3 \tag{④} \]

となっている。\( \gamma_{1} \)は歪度、\( \gamma_{2} \)は尖度と説明されている。

さて、③と④を①に代入する。

\[ b_{2} = \frac {2(\gamma_{2} + 3) - 3\gamma_{1}^{2} - 6} {10(\gamma_{2} + 3) - 12\gamma_{1}^{2} - 18} \] \[ b_{2} = \frac {2\gamma_{2} + 6 - 3\gamma_{1}^{2} - 6} {10\gamma_{2} + 30 - 12\gamma_{1}^{2} - 18} \] \[ b_{2} = \frac {2\gamma_{2} - 3\gamma_{1}^{2}} {10\gamma_{2} - 12\gamma_{1}^{2} + 12} \tag{⑤} \]

そして⑤を①に代入する。

\[ m = \frac {1} {2b_{2}} \] \[ m = \frac {1} {2\frac {2\gamma_{2} - 3\gamma_{1}^{2}} {10\gamma_{2} - 12\gamma_{1}^{2} + 12}} \] \[ m = \frac {1} {\frac {2\gamma_{2} - 3\gamma_{1}^{2}} {5\gamma_{2} - 6\gamma_{1}^{2} + 6}} \] \[ m = \frac {5\gamma_{2} - 6\gamma_{1}^{2} + 6} {2\gamma_{2} - 3\gamma_{1}^{2}} \]

つまり、\( m = 尖度\)ではないことが分かる。

vについて

次に\( v \)について調べる。「Pearson distribution」には

\[ v = -\frac {2b_{2}a+b_{1}} {2b_{2}^{2}\alpha} \tag{①} \]

とある。では\( a \)とは何かというと、

\[ a = \sqrt{\mu_{2}}\sqrt{\beta_{1}}\frac {\beta_{2}+3} {10\beta_{2} - 12\beta_{1} - 18} \tag{②} \]

となっている。

(続く)