PythonでIsotonic regression

Isotonic regressionの応用を考えているが、Rでは少し問題があるのでPythonでやってみようと思って挫折していたのだがなんか急にやる気が出てきたのでやった。
環境はUbuntu 12.04 64bit, Intel Corei5-2540M CPU 2.60 GHz * 4.
　
まずはインストール。

sudo apt-get install libatlas3gf-base
sudo apt-get install libc6
sudo apt-get install libdsdp-5.8gf
sudo apt-get install libfftw3-3
sudo apt-get install libglpk0 
sudo apt-get install libgsl0ldbl
sudo apt-get install python
sudo apt-get install python-central
sudo apt-get install python-cvxopt

Windowsとかサーバーへのインストールめんどくさすぎ…さすがUbuntu…
Windowsはここからインストーラーが入手できる。
Numpyに依存しているらしいので、同じページからNumpyインストーラーを入手して同じく。
Mac OSの場合。
　
入力はとひたすら並んだものを想定している。
20130131追記
上限と下限を設定できるように改良した。
最後のifとforの部分がmapで書けたらよかったのだが。

def IR_py(vec, n_mutant, n_concentration, max0 = float("inf"), min0 = -float("inf")):
	from cvxopt import matrix, solvers
	n_mt = n_mutant #変異株の数
	n_conc = n_concentration #観測濃度点の数
	#1次の項のベクトルは入力
	p1 = map(lambda x: -float(x), vec)
	p = matrix(p1)
	
	#解くべきベクトルの長さは細胞株数*濃度点
	#2次の項の行列を作る
	q1 = []
	for i in range(len(p1)):
		q2 = [0.0]*len(p1)
		q2[i] = 1.0
		q1 += [q2]
	
	Q = matrix(q1)
	
	#不等号制限行列
	wt = []
	for i in range(1, n_conc):
		w1 = [.0]*len(p1)
		w1[i] = 1.0
		w1[i-1] = -1.0
		wt += [w1]
	
	mt = []
	for i in range(n_mt):
		for j in range(1, n_conc):
			w2 = [.0]*len(p1)
			w2[(i+1)*n_conc + j] = 1.0
			w2[(i+1)*n_conc + j - 1] = -1.0
			mt += [w2]
	
	wt_mt = []
	for i in range(n_mt):
		for j in range(n_conc):
			w3  = [.0]*len(p1)
			w3[j] = -1.0
			w3[j + (i+1)*n_mt*n_conc] = 1.0
			wt_mt += [w3]
	
	tr = map(list, zip(*wt + mt + wt_mt))
	G = matrix(tr)
	#G = matrix(wt + mt + wt_mt, (3*n_conc-2, len(p1)))
	#G = matrix(wt + mt + wt_mt)
	#不等式の値
	h = matrix([.0]*(3*n_conc-2))
	#等式はないのでAとbは空
	A = matrix([.0]*len(p1), (1, len(p1)))
	b = matrix(1.0)
	
	res = solvers.qp(Q, p, G, h)["x"]
	for i in range(len(res)):
		if res[i] > max0:
			res[i] = max0
		
		if res[i] < min0:
			res[i] = min0
	
	return list(res)

データ2種類、観測点それぞれ3つのデータを入力してみる。

v = [100.18, 96.94, -0.7577, 96.14, 99.88, 0.0219]

IR_py(v, n_mutant = 1, n_concentration = 3)

#結果
     pcost       dcost       gap    pres   dres
 0: -1.7190e+04 -1.7905e+04  1e+03  4e+00  6e-01
 1: -1.9268e+04 -1.9853e+04  6e+02  4e-02  6e-03
 2: -1.9313e+04 -1.9335e+04  2e+01  6e-04  9e-05
 3: -1.9321e+04 -1.9323e+04  2e+00  3e-06  5e-07
 4: -1.9322e+04 -1.9322e+04  6e-02  5e-14  4e-17
 5: -1.9322e+04 -1.9322e+04  6e-04  3e-14  4e-17
Optimal solution found.
[100.18013202524466, 97.65324701419549, -0.3678801712936925, 97.65338473198302, 97.65323631645902, -0.3679199165885018]

データが0~100までだとすると

IR_py(v, n_mutant = 1, n_concentration = 3, max0 = 100, min0 = 0)

     pcost       dcost       gap    pres   dres
 0: -1.7190e+04 -1.7905e+04  1e+03  4e+00  6e-01
 1: -1.9268e+04 -1.9853e+04  6e+02  4e-02  6e-03
 2: -1.9313e+04 -1.9335e+04  2e+01  6e-04  9e-05
 3: -1.9321e+04 -1.9323e+04  2e+00  3e-06  5e-07
 4: -1.9322e+04 -1.9322e+04  6e-02  5e-14  4e-17
 5: -1.9322e+04 -1.9322e+04  6e-04  3e-14  4e-17
Optimal solution found.
[100.0, 97.65324701419549, 0.0, 97.65338473198302, 97.65323631645902, 0.0]

Pythonで二次計画法その1
Pythonで二次計画法その2
cvxoptの二次計画法の説明文。超分かりやすい。
Pythonの転置