[R言語] Rのお部屋::あーるのおへや[別館]

ライブラリnnetなのですが、ソースを見るとBLASを使っていなかったので、BLASを使うように修正してみました．
ニューロンの数にもよりますが、半分くらいの時間で計算できるようになりました．
くわしくはこちら（Cのソースもあります） ⇒ nnetのソースを修正して計算速度の改善をはかる

バグ報告等は、こちらのエントリーにどうぞ．↓

ATLAS関連の新記事を更新． ⇒BugTrack- Ｒ備忘録/63

ＡＴＬＡＳについては【Rで並列処理 - Atlasを使う】でも触れているのですが、gentoo LinuxのPortageを使ってインストールを行っているので自分でも何がインストールされているのかわかりにくい．なので、ソースからインストールしてみることにしてみた．

こういうチューニングを行ったのが、REvolutionRなのかな・・・

参考サイト ⇒ REvolutionR は連邦の新型か

「Rによるベイズ統計分析」を、Amazonにて、購入．

まだ、読んでません（笑

まさかこの期に及んで、またRserveを使う時がやってくるとは思わなかった・・・

JavaからRserveを使う
JavaからRserveを使う 例②（画像ファイルを扱う方法）
JavaからRserveを使う③

ご質問などはこちらのブログのほうへ、どうぞ

Rのforステートメントは激しく遅いので有名です．この繰り返しをpython側でやっちゃったらどうなるかなっと

# R内で直接まわす
now = datetime.datetime.now()
robjects.r('nn <- 100000;for( ii in 1:nn ){ x <- 4}')
print datetime.datetime.now() - now

# pythonでまわす ①
func = robjects.r('function(x) x <- 4')
now = datetime.datetime.now()
for ii in range(100000):
    tmp = func(4.0)

print datetime.datetime.now() - now

# pythonでまわす ②
func = robjects.r('func <- function(x) x <- 4')
print robjects.r.func
now = datetime.datetime.now()
for ii in range(100000):
    tmp = robjects.r.func(4.0)

print datetime.datetime.now() - now
>

何をしているのか、そのうち思い出します

import rpy2.rinterface as ri
ri.initr()
eval = ri.baseenv["eval"]
expression = ri.parse('x <- 4')

# pythonでまわす ①
now = datetime.datetime.now()
for ii in range(100000):
    tmp = eval(expression)

print datetime.datetime.now() - now
>

parseを使ってみる．

import rpy2.rinterface as ri
ri.initr()
expression = ri.parse('1 + 2')
len(expression)
len(expression[0])
ri.str_typeint(expression[0][0].typeof)
'SYMSXP'
tuple(expression[0][1])
(1.0,)
tuple(expression[0][2])
(2.0,)

>>> ri.baseenv["eval"]
>>> eval = ri.baseenv["eval"]
>>> eval(expression)
>>> print eval(expression)
>>> eval(expression)[0]
3.0
>>> print eval(expression)[0]
3.0
>>> print eval(expression)[1]
Traceback (most recent call last):
File "", line 1, in
IndexError: Index out of range.
>>> print eval(expression)[0]
3.0

>>> expression = ri.parse('3 + 2')
>>> print eval(expression)[0]
5.0
>>> print eval(expression)
<rpy2.rinterface.SexpVector - Python:0xb7fae1c0 / R:0x86da7b0>
>>> tuple(eval(expression))
(5.0,)

Pass-by-value paradigmという、内容が良くわからなかった．

An infamous example is when the column names for a matrix are changed, bringing a system to its knees when the matrix is very large, as the whole matrix ends up being copied

あまりにも悪名高い振る舞いが、matrixオブジェクトのカラム名を変更しようとするときである．Rはmatrixオブジェクトのすべてをコピーしようとするため、matrixが非常に大きい場合、システムが崩壊する．

以下のようにデータを付値すると、リストになってしまう．

>>> robjects.globalenv["a"] = [3,2,1]
>>> print robjects.r.a
[[1]]
[1] 3 

[[2]]
[1] 2

[[3]]
[1] 1

>>> print robjects.IntVector([3,2,1])
[1] 3 2 1

>>> robjects.globalenv["a"] = robjects.IntVector([3,2,1])
>>> print robjects.r.a
[1] 3 2 1

>>> print robjects.r('a')
[1] 3 2 1

>>> print robjects.r('a[1]')
[1] 3

>>> robjects.Vector([3,2,1]).r_repr()
'list(3L, 2L, 1L)'
>>> tmp = robjects.Vector([3,2,1]).r_repr()
>>> tmp
'list(3L, 2L, 1L)'
>>> robjects.r('a <- %s' % tmp)
>>> print robjects.r('a')
[[1]]
[1] 3

[[2]]
[1] 2

[[3]]
[1] 1

import datetime
testdata = range(100000,0,-1)
roInt = robjects.IntVector(testdata)
now = datetime.datetime.now()
tmp = roInt.r_repr()
robjects.r('a <- %s' % tmp)
print datetime.datetime.now() - now

now = datetime.datetime.now()
robjects.globalenv["a"] = roInt
print datetime.datetime.now() - now
>

RをPythonに組み込むことができる．それにインターフェースもちゃんとしてる．こいつはクールだぜ．

rpy2
ドキュメント ： rpy2.robjects

>>> import rpy2.robjects as robjects
>>> robjects.r
<rpy2.robjects.R object at 0xa189ccc>
>>> print robjects.r
<rpy2.robjects.R object at 0xa189ccc>
platform: i686-pc-linux-gnu
arch: i686
os: linux-gnu
system: i686, linux-gnu
status:
major: 2
minor: 15.1
year: 2012
month: 06
day: 22
svn rev: 59600
language: R
version.string: R version 2.15.1 (2012-06-22)
nickname: Roasted Marshmallows
>>> robjects.r.pi
<FloatVector - Python:0xa189e0c / R:0xaa2dff0>
[3.141593]
>>> robjects.r.pi[0]
3.141592653589793
>>> robjects.r.letters
<StrVector - Python:0xa19228c / R:0xa87e2a0>
['a', 'b', 'c', ..., 'x', 'y', 'z']
>>> robjects.r.letters[0:]
<StrVector - Python:0xa19290c / R:0xa87e3d0>
['a', 'b', 'c', ..., 'x', 'y', 'z']
>>> robjects.r.letters[0]
'a'
>>> robjects.r.letters[1]
'b'
>>> robjects.r('pi')
<FloatVector - Python:0xa192aec / R:0xaa2dff0>
[3.141593]
>>> robjects.r('x <- 4')
<FloatVector - Python:0xa192e6c / R:0xaa2ddd0>
[4.000000]
>>> robjects.r('x')
<FloatVector - Python:0xa192dec / R:0xaa2ddd0>
[4.000000]
>>> robjects.r['x']
<FloatVector - Python:0xa1929ec / R:0xaa2ddd0>
[4.000000]
>>> [ee for ee in robjects.r.letters]
['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z']

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

Rのお部屋で自分用のメモとして書いておいた記事（２次計画法（２次計画問題）のメモ）を、少しわかりやすくまとめてみることにしました．

自分用のメモのため、見に来てくれた方にはわかりにくいものとなっていましたので．

Rで二次計画法を行うには「quadprog」というライブラリを使います．

quadprogの説明（別ウインドウでPDFを開きます）にノーテーションを合わせます

min(-d^Tb + 1/2 b^TDb)
with constraints A^Tb ≧ b₀

ここで

b ： 求めたい変数のベクトル n×1の行列（要素n個のベクトル）
D ： n×nの対角行列
d ： n×1の行列（要素n個のベクトル）
A ： 制約条件の左辺の係数
b₀ ： 右辺の定数項

1/2 b^TDb

-d^Tb

制約条件は≧としかなっていませんが、等式制約と不等式制約両方扱えます（あたりまえ）．
ここでquadprogの制約条件は行列Aと制約条件の右辺のb₀で設定しますが
等式制約 → 不等式制約
の順でAとb₀を設定しなくてはなりません．これは後で説明します．

関数solve.QPのサンプルの例を実行してみましょう．最初のコメント部分に以下のように書いてあります．

##
## Assume we want to minimize: -(0 5 0) %*% b + 1/2 b^T b
## under the constraints: A^T b >= b0
## with b0 = (-8,2,0)^T
## and
##     (-4 2 0)
## A = (-3 1 -2)
##     ( 0 0 1)
## we can use solve.QP as follows:
##

これはまず求めたい変数の次元が３で、目的関数が

-(0 5 0) %*% b + 1/2 b^T b

制約条件はbを（b1、b2、b3）とすると、 b0 = (-8,2,0)^Tが右辺、Aが左辺ですから

-4×b1 ＋ -3×b2 ＋       ≧ -8
 2×b1 ＋  1×b2 ＋       ≧  2
          -2×b2 ＋ 1×b3 ≧  0

ここで2次の項を構成する対角行列のDがありませんが

1  0  0
0  1  0
0  0  1

Rのソースを記述します．DmatがD行列を表しています．Dは単位行列ですので以下でDmatを設定できます．

Dmat <- matrix(0,3,3)
diag(Dmat) <- 1
Dmat
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

次に目的関数の1次の項のdvecを設定します．ノーテーションではn×1の行列と書きましたが必ず1列なのでRでいうところのただのベクトルで設定します．

dvec <- c(0,5,0)

Amatは上記制約条件のとおりですが、以下の並びになるように設定します．制約条件の不等式の係数の並びとは行と列が逆になっていますので注意が必要です．

     (-4 2 0)
 A = (-3 1 -2)
     ( 0 0 1)

以下でAmatを設定します．

Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
Amat
     [,1] [,2] [,3]
[1,]   -4    2    0
[2,]   -3    1   -2
[3,]    0    0    1

bvecは制約条件の右辺ですので

bvec <- c(-8,2,0)

実行して

> solve.QP(Dmat,dvec,Amat,bvec=bvec)
$solution
[1] 0.4761905 1.0476190 2.0952381

$value
[1] -2.380952

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.0000000 0.2380952 2.0952381

$iact
[1] 3 2

等式制約がある場合
Rによる二次計画法 例

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

制約条件の順序は、等式制約がある場合前の記事で

等式制約 → 不等式制約

制約条件を以下のように変更してみます（1番目の制約条件を等式制約に変更しました）．

-4×b1 ＋ -3×b2 ＋       ＝ -8
 2×b1 ＋  1×b2 ＋       ≧  2
          -2×b2 ＋ 1×b3 ≧  0

この場合、不等式制約と同じようにRの変数を設定します．

Dmat <- matrix(0,3,3)
diag(Dmat) <- 1
Dmat
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
dvec <- c(0,5,0)
Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
Amat
     [,1] [,2] [,3]
[1,]   -4    2    0
[2,]   -3    1   -2
[3,]    0    0    1
bvec <- c(-8,2,0)

実行するときにmeqという引数に等式制約の数を設定します．今回は1個ですので「meq=1」を設定します．

solve.QP(Dmat, dvec, Amat, bvec=bvec, meq=1)
$solution
[1] 1.123596 1.168539 2.337079

$value
[1] -1.797753

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.2808989 0.0000000 2.3370787

$iact
[1] 3 1

これで等式制約を含んだ2次計画法を実行することができました．

-4 × 1.123596 + -3 × 1.168539 -> -8.000001

等式制約条件も満たしています．

Rによる二次計画法 例

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

2次計画法の典型的な例として、3資産を組み合わせてリスクを最小にするようなポートフォリオの各資産のウエイトを求めてみましょう．

３資産を資産A・資産B・資産Cとします．各資産の例えば月次のリターンを

資産A ：  3%、-1.5%、1%、 7%、4%
資産B ： -2%、 1.5%、7%、 8%、3%
資産C ：  4%、 8.5%、3%、-7%、6%

各資産のウエイトですからそれをW_A、W_B、W_Cと置き、また各資産は必ずロングポジションで保有することにすると、制約条件は

W_A＋W_B＋W_C ＝ 1.0
W_A         ≧ 0.0
    W_B     ≧ 0.0
        W_C ≧ 0.0

それではRで各変数を設定していきましょう．
目的関数はW_A、W_B、W_Cで構成されるポートフォリオのリスクですからまずは共分散行列を求めなければなりません．

retA <- c(0.03,-0.015,0.01,0.07,0.04) #資産Aのリターンの系列
retB <- c(-0.02,0.015,0.07,0.08,0.03) #資産Bのリターンの系列
retC <- c(0.04,0.085,0.03,-0.07,0.06) #資産Cのリターンの系列
ret <- cbind(retA,retB,retC)
C <- cov(ret)
C
           retA       retB       retC
retA  0.0010200  0.0004875 -0.0015475
retB  0.0004875  0.0016750 -0.0015750
retC -0.0015475 -0.0015750  0.0035050

Cをそのまま使っても良いのですがDmatに代入しておきます．

Dmat <- C

目的関数に1次の項はありませんので

dvec <- c(0,0,0)

制約条件の左辺から

## W_A＋W_B＋W_C  → (1,1,1)
## W_A          → (1,0,0)
##     W_B      → (0,1,0)
##         W_C  → (0,0,1)

Amat <- cbind(c(1,1,1), c(1,0,0), c(0,1,0), c(0,0,1))
Amat
     [,1] [,2] [,3] [,4]
[1,]    1    1    0    0
[2,]    1    0    1    0
[3,]    1    0    0    1

制約条件の右辺から

bvec <- c(1,0,0,0)

ここで1番目の制約条件だけが等式制約であることを思い出してください．ですのでmeq引数に1をセットしてsolve.QPを実行します．

res <- solve.QP(Dmat, dvec, Amat, bvec, meq=1, factorized=FALSE)

結果を見てみます．

res
$solution
[1] 0.4624505 0.2148457 0.3227038

$value
[1] 3.852634e-05

$unconstrained.solution
[1] 0 0 0

$iterations
[1] 2 0

$Lagrangian
[1] 7.705267e-05 0.000000e+00 0.000000e+00 0.000000e+00

$iact
[1] 1

$solutionが求めたい変数の結果ですので、これがウエイトになっています．合計すると１になっていて制約条件をちゃんと満たしていますね．

sum(res$solution)
[1] 1

これで、これらの資産のリスクが最小になるウエイトを求めることができました．

cov.wt(x, wt = rep(1/nrow(x), nrow(x)), cor = FALSE, center = TRUE, method = c("unbiased", "ML"))

method = "ML"

とすると、(n-1)で割った場合でなくnで割った値を算出．

z <- rbind(c(0,1), c(0,1), c(0,1), c(0,1), c(0,1))
z2 <- rbind(z, -z)
z2
      [,1] [,2]
 [1,]    0    1
 [2,]    0    1
 [3,]    0    1
 [4,]    0    1
 [5,]    0    1
 [6,]    0   -1
 [7,]    0   -1
 [8,]    0   -1
 [9,]    0   -1
[10,]    0   -1
w2 <- c(rep(1,5), rep(1,5))
cov.wt(z2, w2, method = "ML")

# --- これと同じ共分散マトリックス
z3 <- c(0,1)
z4 <- rbind(z3, -z)
z4
   [,1] [,2]
z3    0    1
      0   -1
      0   -1
      0   -1
      0   -1
      0   -1
w4 <- c(5, rep(1,5))
cov.wt(z4, w4, method = "ML")

確率を行として入れる場合

z5 <- rbind(c(0,0), c(0,0), c(0,1), c(0,1), c(0,1))
z6 <- rbind(z5, -z)
z6
      [,1] [,2]
 [1,]    0    0
 [2,]    0    0
 [3,]    0    1
 [4,]    0    1
 [5,]    0    1
 [6,]    0   -1
 [7,]    0   -1
 [8,]    0   -1
 [9,]    0   -1
[10,]    0   -1
w6 <- c(rep(1,5), rep(1,5))
cov.wt(z6, w6, method = "ML", center = F)
$cov
     [,1] [,2]
[1,]    0  0.0
[2,]    0  0.8

$center
[1] 0

$n.obs
[1] 10

$wt
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

# --- これと同じ共分散マトリックス
z7 <- c(0,sqrt(0.6))
z8 <- rbind(z7, -z)
z8
   [,1]       [,2]
z7    0  0.7745967
      0 -1.0000000
      0 -1.0000000
      0 -1.0000000
      0 -1.0000000
      0 -1.0000000
w8 <- c(5, rep(1,5))
cov.wt(z8, w8, method = "ML", center = F)
$cov
     [,1] [,2]
[1,]    0  0.0
[2,]    0  0.8

$center
[1] 0

$n.obs
[1] 6

$wt
[1] 0.5 0.1 0.1 0.1 0.1 0.1

各行が０or１のフラグだと考え、その確率を行として加えたい場合は、上記のようにする．
分散は平均周りの二次のモーメントであるため1行目から5行目までの確率0.6になるように

sqrt(0.6)

をセットし、そのときのウエイトを５とすることで、同じ共分散マトリックスが得られる．

ただし、

center = F

の場合に限る．

なぜかというと、「center=T」では、各列の平均を計算し平均で各列を基準化するのだが、z6とz8で平均が同じにならないため．

いまさらですが、多次元正規分布．

ライブラリmvtnormで、多次元正規分布の密度関数などが提供されています．

> dmvnorm(x=rbind(c(0,0),c(0,1)), mean=c(1,1))
[1] 0.05854983 0.09653235

マトリックスをxとして入力すると、各行の密度を返してくれます．