ピアソン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{②} \]

となっている。

(続く)

lognorm()関数の引数指定 (2017/04/12)

SciPyのlognorm()関数は引数の指定が分かりにくい。

例えば、LibreOffice Calcで適当なセルに「=LOGNORMDIST(4,3.5,1.2,1)」と入力すると「0.0390835557」と返される。

「4」は調べたい数値(対数にする必要はなく、そのままでよい)、「3.5」はデータを対数に変換した後の平均、「1.2」はデータを対数に変換した後の標準偏差、「1」は戻り値として累積分布を指定、という意味である。書式はExcelと同じである(Excelの場合は「LOGNORM.DIST」()関数)。

これと同じことをSciPyでやるには

import numpy as np
from scipy.stats import lognorm
lognorm.cdf(x=4, s=1.2, loc=0, scale=np.exp(3.5))

のようにする。

0.039083555706800471

つまり、「x」には調べたい数値、「s」にはデータを対数に変換した後の標準偏差、「loc」には0、「scale」にはデータを対数に変換した後の平均を指数に変換(ややこしい)した数値を指定するのである。

SciPyのnorm()関数では

norm.cdf(x=1.0, loc=0, scale=1)

とした場合、xは同じだが、「loc」は平均、「scale」は標準偏差を意味する。ここに注意しないと間違える。

Twitter botの作成② (2017/04/11)

サンプルプログラムの作成

Twitter APIには「REST API」、「Streaming API」、「Sample API」などがある

。REST APIは過去のツイートを検索できる機能もあって、よく使われていると思う。だが、検索にかかるのはほんの一部なので、検索結果を元に統計分析を行うのには適していない。

Streaming APIはリアルタイムですべてのツイートが流れてくる。だが、日本語での検索がうまくいかないようだ。

そこで今回はSample APIを使ってみる。

Sample APIはツイート全体の1%しか利用できない。だが、その1%が無作為に選ばれているのであれば(どのように選ばれているかは知らないが)、それなりに全体を代表していると言えるだろう。1%といってもデータは大量で、REST APIとは比較にならない。

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

「# キー、アクセストークンを設定する。」のブロックは「Twitter botの作成①」でメモしたものに書き換える。

# coding: utf-8

import sys
import time
import twitter
from datetime import datetime

# キー、アクセストークンを設定する。
#consumer_key = "*************************"
#consumer_secret = "**************************************************"
#access_token_key = "**********-***************************************"
#access_token_secret = "*********************************************"

if __name__ == "__main__":
    api = twitter.Api(
        consumer_key=consumer_key, consumer_secret=consumer_secret,
        access_token_key=access_token_key,
        access_token_secret=access_token_secret)
    argvs = sys.argv
    query = argvs[1]
    count = 0
    i = 0
    while True:
        # 回線不通の場合を処理する。
        try:
            for line in api.GetStreamSample():
                # 空のデータが流れてきた場合を処理する。
                try:
                    now = datetime.now()
                    count = count + line["text"].count(query)
                    interval = now.minute % 1 == 0 and now.second == 0
                    if i == 0 and interval:
                        print(
                            now.strftime("%Y-%m-%d %H:%M:%S"), query, " = ",
                            count, "回")
                        count = 0
                        i = 1
                    if interval == False:
                        i = 0
                except:
                    # 何もしたくないのでパスする。
                    pass
        except:
            print("回線不通です。")
            time.sleep(60)

使い方

記事を書いている時点で話題になっているのは浅田真央選手の引退だろうか。そこで「浅田真央」というキーワードが使われた回数を5分ほど数えてみた。

以下のコマンドを実行してTwitter botで「浅田真央」を検索する。

%run ~/py/test_bot.py 浅田真央


2017-04-11 23:04:00 浅田真央  =  3 回
2017-04-11 23:05:00 浅田真央  =  3 回
2017-04-11 23:06:00 浅田真央  =  5 回
2017-04-11 23:07:00 浅田真央  =  1 回
2017-04-11 23:08:00 浅田真央  =  1 回

「浅田真央」というキーワードは5分間で13回使われた。これが全体の1%だとすると、単純計算で1300回。1日に換算すると35万回を超える。かなりの人が浅田真央選手の引退に関心を持っているように見える。

参考リンク

python-twitterのドキュメント

https://python-twitter.readthedocs.org/en/latest/

Twitter botの作成① (2017/04/11)

アプリ作成者としての登録

Twitterでbotを使用するには先ずアプリ作成者として登録する必要がある。

①Twitterのアカウントを作成する。

アカウント作成の手順はあえて書く必要もないと思うので省略する。Twitterをやっている人なら、そのアカウントを使ってもいい。だが、もしそのアカウントに携帯電話番号を登録していない場合は登録する必要がある。

②以下のアドレスをクリックする。

https://apps.twitter.com/

③右上隅の「Sign in」をクリックする。

④「Twitterにログイン」でIDとパスワードを入力し、「保存する」にチェックを入れ、「ログイン」をクリックする。

⑤「Create New App」ボタンをクリックする。

⑥「Application Details」で「Name」にアプリの名前、「Description」にアプリの説明、「Website」にウェブサイトのアドレスを入力する。

Webアプリを作るのでなければ、ウェブサイトは実在しないアドレスでも構わないようだ。

⑦「Developer Agreement」で「Yes, I have read and agree to the Twitter Developer Agreement.」にチェックを入れる。

⑧「Create your Twitter application」ボタンをクリックする。

キーとアクセストークンの取得

上に続けて実行する。

①「Keys and Access Tokens」タブをクリックする。

②「Consumer Key (API Key)」と「Consumer Secret (API Secret)」をメモする。

③「Create my access token」をクリックする。(※)

④「Your Access Token」で「Access Token」と「Access Token Secret」をメモする。(※)

※ 私が取得したときはこのようであったが、以来、新しいアプリは作っていないので未確認。

python-twitterのインストール

python-twitterはPythonでTwitterのAPIを動かすライブラリである。他にも多くのライブラリがあって、どれがよいか分からない。とりあえず、更新頻度が高く、活発そうなのでこれにした。

以前、通常バージョンはPython3では正常に動かなかったので(今はどうか確認していない)、開発バージョンをインストールする。

$ pip install --upgrade git+git://github.com/bear/python-twitter.git

ファイル名が一部不明な場合のリネーム (2017/04/08)

ファイル名が一部不明な場合にリネームするコードをPythonで書いてみる。

例えば「text@@@@.txt」というファイルがあるとする。これを「text.txt」にリネームしたい。だが、ファイル名が「text」から始まっているのは知っているが、「@@@@」が何か分からない場合、どうするか。

ワイルドカードを使って、「text*.txt」としてもエラーになる。「*.txt」ならエラーにはならないが、それでは拡張子がtxtのすべてのファイルが対象になってしまう。簡単そうで、意外と難しい(私だけ?)。

あれこれ試行錯誤した結果、あまりすっきりはしないが、一応の方法を見つけた。例として、作業ディレクトリに「text1234.txt」ファイルがあるとして、これを「text.txt」にリネームする。だが、「1234」の部分が正確にはどうであるかは分かっていない。

import glob
import os
new_name = 'text.txt'
for old_name in glob.glob('text*.txt'):
    os.rename(old_name, new_name)


「text1234.txt」ファイルが「text.txt」ファイルにリネームされる。「なんだ、「text*.txt」でもできるじゃないか」と思われるかもしれない。だが、これはglobを噛ませているからできるのである。

単にos.rename('text*.txt', new_name)とやるとエラーになる。