Search on the blog

2015年2月28日土曜日

matplotlibでポワソン分布のヒストグラムを描画

 最近matplotlibというPythonのグラフ描画ライブラリをインストールしたので、ちょっと遊んでみた。まずは簡単なものからということで、ポワソン分布のヒストグラムを描画してみた。
  1. numpy.random.poisson
  2. 自作のMCMC(メトロポリス法)
の2通りの方法で乱数を発生させてヒストグラムにしてみた。

自作MCMCはランダムウォークするときの次の数の選び方がかなりテキトーだったのだが、それっぽい結果になっていた。メトロポリス法はダメで、メトロポリス・ヘイスティングを使わないといけないというシナリオを期待していたのだけど...

グラフ描画の処理
公式ドキュメントのヒストグラム作成サンプルを参考に書いた。
rpoisson(lam, n)がポワソン分布に従う乱数を生成する関数。
hist()のビンの数、幅はうまく整数値がはまるように調整した。
import math
import matplotlib.pyplot as plt
from metropolis import rpoisson

# pdf of poisson distibution (mean = lam) at k
def p(lam, k):
    ret = math.exp(-lam)
    for i in xrange(k):
        ret *= 1.0 * lam / (i + 1)
    return ret

n = 1000000  # number of samples
lam = 10      # mean 
samples = rpoisson(lam, n)  # get n samples

maxv = max(samples)
plt.hist(samples, range=(-0.5, maxv+0.5), bins=maxv+1, normed=1, facecolor='blue', alpha=0.5)

y = [p(lam, x) for x in range(maxv+1)]
plt.plot(range(maxv+1), y, 'r')
plt.xlabel('x')
plt.ylabel('p(x)')
plt.title(r'Histogram of x ~ Poisson distribution($\lambda$ = %d)' % lam)

plt.subplots_adjust(left=0.15)
plt.show()

numpyの関数を使って乱数生成
import numpy as np

def rpoisson(lam, n):
    return np.random.poisson(lam, n)

lam = {1, 4, 10}と変化させて乱数を生成した。サンプル数は1,000,000。
以下がヒストグラム。赤色の線は理論値。



ポワソン分布にしたがって乱数が生成されていることが分かる。

自作MCMCを使って乱数生成
import math
from random import randint, random
from itertools import islice

def rpoisson(lam, n):
    
    def p(lam, k):
        if k < 0:
            return 0.0
        ret = math.exp(-lam)
        for i in xrange(k):
            ret *= 1.0 * lam / (i + 1)
        return ret

    def metropolis(lam):
        x = lam
        while True:
            yield x
            xt = (x-1) if randint(0, 1) == 0 else (x+1)
            alpha = p(lam, xt) / p(lam, x)
            if random() < alpha:
                x = xt

    return list(islice(metropolis(lam), n))
次の乱数候補には、50%ずつの確率で1大きい数 or 1小さい数を採用している。
おそらくこのような単純なやり方ではダメだろうと思っていたが、意外とそれらしい感じで乱数を生成できた。




2015年2月22日日曜日

キャッシュヒット率の測定

 perfというツールを使ってキャッシュヒット率を測定してみた。
キャッシュの影響なんてたかが知れていると甘くみていたら、痛い目るので要注意。

キャッシュの用語
知らない用語がちらほらあったので、簡単にまとめておく。

L1、L2、L3キャッシュ
Level1、Level2、Level3の略。
それぞれ速度、容量に違いがある。

速度L1 > L2 > L3
容量L1 < L2 < L3

高速なキャッシュほど容量が小さくヒット率は低い。
容量の大きいキャッシュはヒット率は高いが、低速。

というトレードオフがある。

プロセッサは以下の順序で処理実行に必要なデータを探す。
①L1キャッシュにアクセスし、ヒットすれば処理実行。ミスの場合は、L2へ。
②L2キャッシュにアクセスし、ヒットすれば処理実行。ミスの場合は、L3へ。
③L3キャッシュにアクセスし、ヒットすれば処理実行。ミスの場合はメインメモリへ。

L1dキャッシュとL1iキャッシュ
L1dはLevel 1 data cache。
L1iはLevel 1 instruction cache。

data cacheは、データをキャッシュするためのキャッシュ。

instruction cacheはプロセッサによって実行される命令をキャッシュするためのキャッシュ。

命令をキャッシュするとはどういうことか?

ストアドプログラム方式では実行されるプログラムもメモリに格納されているため、プロセッサは実行する命令をメモリから読まなければいけない。
よってプロセッサから見ると命令もある種のデータなのである。

キャッシュサイズを参照
Debian系マシンではlscpuコマンドでキャッシュサイズを確認できる。
L1-L3トータルで3MB以上ある。予想してたよりも容量大きい。

kenjih$ lscpu
Architecture:          i686
CPU 操作モード:   32-bit, 64-bit
Byte Order:            Little Endian
CPU(s):                4
On-line CPU(s) list:   0-3
コアあたりのスレッド数:2
ソケットあたりのコア数:2
Socket(s):             1
ベンダー ID:       GenuineIntel
CPU ファミリー:   6
モデル:             42
ステッピング:    7
CPU MHz:               800.000
BogoMIPS:              4988.44
仮想化:             VT-x
L1d キャッシュ:   32K
L1i キャッシュ:   32K
L2 キャッシュ:    256K
L3 キャッシュ:    3072K

キャッシュヒット率を計測
Linuxのパフォーマンス解析ツールperfを使って、キャッシュヒット率を計測してみた。 例として、以下の2つの方法で行列の掛け算を行ったときのキャッシュヒット率を計測する。

プログラムA (a.cpp)
#include <iostream>
#include <vector>

using namespace std;

typedef vector<int> vec;
typedef vector<vec> mat;

const int N = 1024;

mat multiply(const mat &x, const mat &y) {
    int r = x.size();
    int m = y.size();
    int c = y[0].size();

    mat z(r, vec(c));
    for (int i = 0; i < r; i++) {
        for (int j = 0; j < c; j++) {
            for (int k = 0; k < m; k++) {
                z[i][j] += x[i][k] * y[k][j];
            }       
        }
    }

    return z;
}

int main(int argc, char **argv) {
    
    mat x(N, vec(N));
    mat y(N, vec(N));

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            x[i][j] = i;
            y[i][j] = j;
        }
    }

    mat z = multiply(x, y);

    long long sum = 0;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            sum += z[i][j];
        }
    }

    cout << sum << endl;

    return 0;
}

プログラムB (b.cpp)
#include <iostream>
#include <vector>

using namespace std;

typedef vector<int> vec;
typedef vector<vec> mat;

const int N = 1024;

mat multiply(const mat &x, const mat &y) {
    int r = x.size();
    int m = y.size();
    int c = y[0].size();

    mat z(r, vec(c));
    for (int i = 0; i < r; i++) {
        for (int k = 0; k < m; k++) {
            for (int j = 0; j < c; j++) {
                z[i][j] += x[i][k] * y[k][j];
            }       
        }
    }

    return z;
}

int main(int argc, char **argv) {
    
    mat x(N, vec(N));
    mat y(N, vec(N));

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            x[i][j] = i;
            y[i][j] = j;
        }
    }

    mat z = multiply(x, y);

    long long sum = 0;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            sum += z[i][j];
        }
    }

    cout << sum << endl;

    return 0;
}

プログラムAとプログラムBの違いは、multiply関数のループの部分だけである。
プログラムAはi, j, kの順で、プログラムBはi, k, jの順でループを回す。

一見大した違いは無さそうであるが、キャッシュのヒット率によって大きな違いが生じる。

コンパイルは以下のように行った。

g++ -Wall -O2 a.cpp -std=c++0x -lm -o a

実行してみると、

プログラムA 9.260s
プログラムB 1.184s

となった。
約8倍ほどプログラムBの方が高速。

それぞれのプログラムでキャッシュ参照数、キャッシュミス数などを以下のコマンドで出力してみた。

perf stat -e cycles,instructions,cache-references,cache-misses ./a
結果は以下のとおり。

 プログラムA
 Performance counter stats for './a':

    27,363,421,784 cycles                    #    0.000 GHz                    
    11,880,422,836 instructions              #    0.43  insns per cycle        
        71,384,828 cache-references                                            
        57,242,356 cache-misses              #   80.188 % of all cache refs    

プログラムB
Performance counter stats for './b':

     3,618,482,791 cycles                    #    0.000 GHz                    
     6,495,956,182 instructions              #    1.80  insns per cycle        
         9,516,922 cache-references                                            
         4,136,614 cache-misses              #   43.466 % of all cache refs    

プログラムAでは、プログラムBの約13.8倍のキャッシュミスが発生していた。 つまりメモリアクセスの回数が13.8倍多いということ。 メモリアクセスがボトルネックとなって、1サイクルあたりの命令実行数にも差が生じていることが分かる。

元気が出る動画(3)シゴトでココロおどろう


「好きなことを仕事にする。」
「好きなことだけを続ける。」
「全ては後からついてくる。」
「計画を立てない。」
「目の前のことを全力でやってみてダメだったら次へ。」
「シゴトでココロおどろう。」

2015年2月15日日曜日

中国剰余定理のプログラム

 中国剰余定理のプログラム(連立合同式を解くプログラム)を書いてみた。

中国剰余定理とは?
2で割ると1余り、3で割ると2余り、5で割ると3余る数なーんだ?
というような問題を考える。

この問題を式で表すと、以下のような連立合同式になる。

x = 1 (mod 2)
x = 2 (mod 3)
x = 3 (mod 5)

中国剰余定理によると、
法が互いに素であれば、法の総積を法としてxは一意に定まる。

以下ではそのようなxを効率的に求めることを考える。

解法
法が互いに素な場合の解法は、以下のとおり。

x = a0 (mod m0)
x = a1 (mod m1)

を解きたいとする。
(m0, m1) = 1であるため、

p m0 + q m1 = 1

となるような、整数p, qが存在する。

x = a0 q m1 + a1 p m0

とすれば、xは上の方程式の解である。

これは以下のように確認できる。

x = a0 (1 - p m0) + a1 p m0 = a0 - a0 p m0 + a1 p m0 = a0 (mod m0)
x = a0 q m1 + a1 (1 - q m1) = a0 q m1 + a1 - a1 q m1 = a1 (mod m1)

次にm0, m1が互いに素ではない場合を考える。この場合は、m0、m1が互いに素になるように素因数を片方に寄せればよい。

g = (m0, m1)

とし、

g = p0^n0 p1^n1 p2^n2 .... pk^nk ...

と素因数分解されるとする。

素因数pkの個数を考える。m0の方がm1よりpkを多く持っていれば、pkはすべてm0に寄せる。そうでない場合は、pkはすべてm1に寄せる。

これをすべての素因数について行えば、最終的に(m0, m1) = 1となる。
よってm0, m1が互いに素な場合に帰着する。

プログラム
def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a%b)

# return [g, x, y]
# g = gcd(a, b)
# x, y satisfies a x + b y = g
def extgcd(a, b):
    if b == 0:
        return [a, 1, 0]
    g, x, y = extgcd(b, a%b)
    return [g, y, x - a/b * y]

# eq0: x = a0 (mod m0)
# eq1: x = a1 (mod m1)
# returns [xt, mod] such that x = xt + k mod for integer k.
def crt(eq0, eq1):
    a0, m0 = eq0
    a1, m1 = eq1

    g = gcd(m0, m1)

    if a0 % g != a1 % g:
        raise Exception("x doesn't exist.")

    if g > 1:
        m0 /= g
        m1 /= g

        while True:
            gt = gcd(m0, g)
            if gt == 1:
                break
            m0 *= gt
            g /= gt
        
        m1 *= g

        a0 %= m0
        a1 %= m1

    g, p, q = extgcd(m0, m1)
    
    x = a0 * q * m1 + a1 * p * m0
    mod = m0 * m1
    x = x % mod

    return [x, mod]

テスト
テストコード。
# for debug use only
def test():
    for x in range(1000):
        for m1 in range(1, 1000):
            for m2 in range(1, 1000):
                a1 = x % m1
                a2 = x % m2

                y, mod = crt((a1, m1), (a2, m2))
                
                assert x % mod == y % mod

    print "test success"

サンプル
標準入力から連立合同方程式を受け取って、解くプログラム。
解を10個出力するようにしている。
# read system of linear congruence
n = int(raw_input('input # of system:'))
eqs = []
for _ in xrange(n):
    a, m = map(int, raw_input().split())
    eqs.append((a, m))

# solve the system
x, mod = reduce(crt, eqs, (0, 1))

# print the answer
for _ in range(10):
    print x
    x += mod
まず、 連立合同式
x = 1 (mod 2)
x = 2 (mod 3)
x = 3 (mod 5)
を解いてみる。
input # of system:3
1 2
2 3
3 5
x = 
23
53
83
113
143
173
203
233
263
293

次に、法が互いに素ではない連立式を解いてみる。
input # of system:3
15471 102102
357813 504504
24691357 98765432
x = 
123456789
105883678906461
211767234356133
317650789805805
423534345255477
529417900705149
635301456154821
741185011604493
847068567054165
952952122503837

2015年2月13日金曜日

漸化式を解く5つの方法

 漸化式を解く5つの方法を纏めた以下のページがなかなかよかった。

How to Solve Recurrence Relations

漸化式の形が

1. an = an-1 + d
2. an = an-1 * r
3. an = an-1 + p(n),  p(n)はnの多項式
4. an = 前のk項の線形結合

の場合に、どのように一般項を導出すればよいかが記載されている。

1.と2.はarithmetic/geometric sequenceなので高校数学で習う。
3.は知らなかった。anはp(n)より1つ高次の多項式になるらしい。
4.は高校のとき習った。特性方程式を解けばいい。

最後に、
5. 母関数を使って一般項を出す方法
も紹介されている。

せっかくなので、4.と5.の方法を使ってフィボナッチ数列の一般項を求めてみた。
4. で解くほうが断然楽だった。

2015年2月11日水曜日

mod Nで等差数列/等比数列の和を計算

最近知ったちょいテクをメモっておく。

等差数列の和
等差数列の和をmod Nで計算することを考える。
初項a, 等差d, 項数nとすれば、求める和は、

sum = n (2a + (n-1) d) / 2

である。この値はmod Nではどうなるか?
gcd(2, N) = 1であれば、分子に2の逆元をかければいい。

Nが偶数の場合はどうすればいいか?一般に q | pのとき

p / q (mod N) = p (mod N * q) / q

が成り立つ。

よって、(分子 mod 2N) / 2を計算すればよい(※)。
ll sum(ll a, ll d, ll n, ll mod) {
    mod *= 2;
    ll ret = (2 * a + d * (n-1)) % mod;
    ret = ret * n % mod;
    return ret / 2;
}
※実は、n、2a + (n-1) dのいずれかは必ず2で割り切れるので、予め割れる方を割っておけば赤字のテクニックを知らなくても計算できる。

等比数列の和
初項a, 等比r, 項数nとする。求める和は、

sum = a(r^n - 1) / (r - 1)

である。法Nの値によっては、(r-1)の逆元が存在しない場合もある。
この場合も、赤字のテクニックを使うことで、和を計算できる。

または、n項の和をn/2項ずつの和に分割して計算することも出来る。
ll pow(ll x, ll p, ll mod) {
    ll ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

ll sum(ll a, ll r, ll n, ll mod) {
    if (n == 1)
        return a % mod;

    ll x = sum(a, r, n/2, mod);
    ll ret = (x + pow(r, n/2, mod) * x) % mod;
    if (n & 1)
        ret = (a + r * ret) % mod;
    return ret;
}

2015年2月10日火曜日

フィボナッチ数列は法nで周期的

 フィボナッチ数列は法nで周期的である。
例えば、フィボナッチ数の下1桁を表示すると

0 1 1 2 3 5 8 3 1 4 5 9 4 3 7 0 7 7 4 1 5 6 1 7 8 5 3 8 1 9 0 9 9 8 7 5 2 7 9 6 5 1 6 7 3 0 3 3 6 9 5 4 9 3 2 5 7 2 9 1 0 1 1 2 3 5 8 3 1 4 5 9 4 3 7 0 7 7 4 1 5 6 1 7 8 5 3 8 1 9 0 9 9 8 7 5 2 7 9 6 ....

のように周期60の数値が並ぶ。これはフィボナッチ数列が法10で周期的であることを意味している。同様に、他のnに対してもフィボナッチ数列は法nで周期的となる。

この理由について考えてみた。

まず、mod nの要素はn個しかない。よって、フィボナッチ数列の隣合う数(f_n, f_n+1)をmod n で考えると、pigeon hole principleにより、n*n+1回目以内に既出の点に戻ることになる。よって、f_n mod nは最終的にcycleに入る。

ここまでは簡単。
問題はただcycleに入るだけじゃなくて、(f_0, f_1) mod nに戻ってくるというのをどのように示せばよいかだ。

分からなかったので、ググってみると以下のような証明が見つかった。

Proof that Fibonacci Sequence modulo m is periodic?

なるほど。

途中からサイクルに入るような列をpreperiodic(前周期的)というらしい。

一般に、有限サイズの集合Sに対して、f: S → Sを定義すると、Sの要素に対して、fを繰り返し適用してえられる列はpreperiodicになる。
つまり、あるtが存在して、十分大きなnに対して、f(n + t) = f(n)が成り立つ。

ここで、fがbijectiveな場合は、上式の両辺にf^-1を繰り返し適用すると、いずれf(t) = f(0)となるので、(strictly) periodicになる。
フィボナッチ数列の式f(a, b) = (b, a+b)はbijectiveなので、(strictly) periodicということらしい。

フィボナッチ数列だけじゃなくて、bijectiveな関数なら何でもいいらしい。

ということで、ちょっと実験してみる。
mod 10でf(x) = ax+7を考えてみる。gcd(a, 10) = 1のときは、逆関数f^-1(x) = a^-1 (x-7)が存在するので、(strictly)periodicになることを確認する。
ついでに、開始点を0にすることで、gcd(a, 10) > 1のときは0に戻ってこないことを確認する。(∵ gcd(a, 10) > 0のときは0はfのimageには含まれないため。)

def gen(a, mod):                                                                
    x = 0                                                                       
    while True:                                                                 
        yield x                                                                 
        x = (a * x + 7) % mod                                                   
                                                                                
mod = 10                                                                        
for a in range(mod):                                                            
    print "a = %d:" % a,                                                        
    g = gen(a, mod)                                                             
    for _ in range(mod + 2):                                                    
        print g.next(),                                                         
    print

以下実行結果。

a = 0: 0 7 7 7 7 7 7 7 7 7 7 7                                                
a = 1: 0 7 4 1 8 5 2 9 6 3 0 7                                                
a = 2: 0 7 1 9 5 7 1 9 5 7 1 9                                                
a = 3: 0 7 8 1 0 7 8 1 0 7 8 1                                                
a = 4: 0 7 5 7 5 7 5 7 5 7 5 7                                                
a = 5: 0 7 2 7 2 7 2 7 2 7 2 7                                                
a = 6: 0 7 9 1 3 5 7 9 1 3 5 7                                                
a = 7: 0 7 6 9 0 7 6 9 0 7 6 9                                                
a = 8: 0 7 3 1 5 7 3 1 5 7 3 1                                                
a = 9: 0 7 0 7 0 7 0 7 0 7 0 7

予想通りになってた。

2015年2月5日木曜日

Pythonのデコレータで遊ぶ

 Pythonのデコレータで遊んでみた。Javaのアノテーションでアスペクトを挿入するようなイメージに近いのかなと思った。Pythonのデコレータの方がお手軽感がある。

関数の実行前後にログを出す
まず考えてみた例。関数の実行前と実行後に関数名 + "begin/end."と表示させるデコレータ。
def log(f):
    def _f():
        print "%s begin." % (f.__name__)
        f()
        print "%s end." % (f.__name__)
    return _f 

@log
def func():
    pass

if __name__ == '__main__':
    func()

関数をn回繰り返す
次に考えた例。
受け取った数だけ関数を繰り返し実行するデコレータを返す関数を使って、こんなこと出来るのだろうかと思ったら出来た。
def repeat(n):
    def _repeat(f):
        def _f():
            for _ in range(n):
                f()
        return _f
    
    return _repeat

@repeat(5)
def sayHello():
    print "Hello!"

if __name__ == '__main__':
    sayHello()

再帰関数をメモ化する
最後。可変長リストを使って、任意の再帰関数をメモ化するデコレータ。
def memorize(f):
    cache = {}
    def _f(*args):
        if args not in cache:
            cache[args] = f(*args)
        return cache[args]
    return _f

@memorize
def tarai(x, y, z):
    if x <= y:
        return y
    return tarai(tarai(x-1,y,z), tarai(y-1,z,x), tarai(z-1,x,y))

@memorize
def fib(n):
    if n < 2:
        return n
    return fib(n-1) + fib(n-2)

if __name__ == '__main__':
    print fib(100)
    print tarai(100, 50, 0)