ちょうど待ち行列MAP/G/1の系内客数分布を数値計算しようとしていて, どうせなので汎用性高い物にしようと思ったの動機.
まぁありきたりな話, 新しいライブラリを使うと大抵最初はまる.
・numpy のインストール
と思ったらすでにnumpy は入っていた. バージョンは1.2.
古いようなので, easy-install でupgrade する.
% sudo easy_install --upgrade numpyとここで, 鬼のようにwarning と error っぽいメッセージが. 萎える.
・numpy を使う
matrixってのが行列の定義のようだ.
>>> C_matrix = [
... [-0.2,0.2],
... [0.0, -0.2]
... ]
>>> C_matrix = matrix(C_matrix)
>>> D_matrix = [
... [0,0],
... [0.2,0]
... ]
>>>
>>> C_matrix = matrix(C_matrix)
>>> D_matrix = matrix(D_matrix)
>>> C_matrix
matrix([[-0.2, 0.2],
[ 0. , -0.2]])
>>> D_matrix
matrix([[ 0. , 0. ],
[ 0.2, 0. ]])
のように定義できた. array ってのもあるが, これは完全に転置っぽいので直感的でないため, やめる.
(注) matrix で[[1,2],[3,4]] は array では [[1,3],[2,4]] と列で定義するとうこと.
演算もらくらくで便利そうだ.
>>> E_VEC = matrix([[1],]*2)
>>> E_VEC
matrix([[1],
[1]])
>>> C_matrix + D_matrix
matrix([[-0.2, 0.2],
[ 0.2, -0.2]])
>>> (C_matrix + D_matrix) * E_VEC
matrix([[ 0.],
[ 0.]])
>>> dot((C_matrix + D_matrix), E_VEC)
matrix([[ 0.],
[ 0.]])
* は各成分の掛け算になるらしい. dot()で行列のかけざん.
で, solve ってのが連立方程式の解のよう.
とりあえず, 行ベクトルPI * (C + D) = 0, PI*e = 0 を解いてCおよびDで構成される, (C+Dをジェネレータとして持つ)マルコフ連鎖の定常状態確率を解くことにする.
ここではまった.
solve は正方行列でないと扱えないらしく, おそらくフルランクでないと無理だろう.
とりあえず, C+Dは明らかにフルランクではなく, PI*e = 0と組み合わせることでようやく解ける.
どうしようかと思ったけどとりあえず暫定的に, C+Dの一列を抜いてeへ替える. そして右辺は[0,...,0,1]とおく.
以下のようになる.
from numpy.linalg import solve
def calculate_pi():
sumC_D = C_matrix + D_matrix
pi_solver = []
for column in sumC_D.tolist():
column.append(1)
pi_solver.append(column[1:])
zeros_plus1 = [0] * (M_SIZE -1) + [1]
pi_solver = mat(pi_solver)
zeros = mat(zeros_plus1)
あまりに煩雑すぎるが, 特にメソッドもなさそうなので暫定的に...
pi_solverは最左行を省く. (まぁこれがのちのち問題だと思う)
zeros はさっきかいた[0,...,0, 1]の行ベクトル.
で, 早速solve.
>>> pipi = solve(t,zeros) Traceback (most recent call last): File "", line 1, in File "/Library/Python/2.6/site-packages/numpy-1.5.1-py2.6-macosx-10.6-universal.egg/numpy/linalg/linalg.py", line 316, in solve raise LinAlgError, 'Incompatible dimensions' numpy.linalg.linalg.LinAlgError: Incompatible dimensions
はーなんでやーと思って, 色々やってみると, 同じことを matrix のかわりに arrayなら通る.
で, もしかたしらとちょっと不安だったが, どうやら
行* 行列 = 行
はだめで
行列 * 列 = 列
しかやってくれないぽい.
ふざけんなよーと思って, またもや暫定的に
def calculate_pi():
sumC_D = C_matrix + D_matrix
pi_solver = []
for column in sumC_D.tolist():
column.append(1)
pi_solver.append(column[1:])
zeros_plus1 = [0] * (M_SIZE -1) + [1]
pi_solver = mat(pi_solver)
zeros_plus1 = mat(zeros_plus1)
pi_vec = solve(pi_solver.T,zeros_plus1.T)
return pi_vec.T
とした. Tで転置らしい. これは先が思いやられるな...
