2014-07-12

RでLPを解く

RでもLP程度なら「不等式制約付き最適化問題」であっても
数値計算で求められることがわかりました。
用いる関数は、optim()のラッパー関数であるconstrOptim()を使います。

例えば、以下のような(最大化)LPを考えます。

目的関数:
　 $3x_1 + 2x_2$
制約条件:
　 $3x_1 + x_2 \le 9$
　 $2.5x_1 + 2x_2 \le 12.5$
　 $x_1 + 2x_2 \le 8$
　 $x_1, x_2 \ge 0$

(解答は, $(x_1, x_2)=(2, 3)$ のとき, 最適値12)

わかりやすさのためにまず、constrOptim()の引数のために目的関数f、制約条件uiとciを定義しておきます。fに関しては、optim()と同様に扱えます。

> # 目的関数
> f <- function(par) {
+   return(3*par[1] + 2*par[2])
+ }

uiとciに関しては、「 $ui$ %*% $par \ge ci$ 」を満たすように構成します。
parは推定するパラメータです。

> # 制約条件その1
> (ui <- matrix(c(-3, -1, -2.5, -2, -1, -2, 1, 0, 0, 1), nc=2, byrow=TRUE))
     [,1] [,2]
[1,] -3.0   -1
[2,] -2.5   -2
[3,] -1.0   -2
[4,]  1.0    0
[5,]  0.0    1
> # 制約条件その2
> (ci <- c(-9, -12.5, -8, 0, 0))
[1]  -9.0 -12.5  -8.0   0.0   0.0

さて、目的関数と制約条件の定義が済んだので、実行に入ります。

引数thetaには推定パラメータの初期値を与え、
またこの関数はデフォルトで目的関数の最小化になっているので、control=list(fnscale=-1)と加えることで、最大化にしています。
引数gradはfのグラディエント関数で、これにNULLを与えるとNelder-Mead法、関数を与えると準ニュートン法による計算となります。

> # 実行（ Nelder-Mead法 ）
> # constrOptim(theta=初期値, f=目的関数, grad=fのグラディエント関数, ui, ci, control=optim()へ渡すパラメータ)
> constrOptim(c(1, 1), f, grad=NULL, ui=ui, ci=ci, control=list(fnscale=-1))
$par
[1] 2 3

$value
[1] 12

$counts
function gradient 
     538       NA 

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 5

$barrier.value
[1] 0.002829033

結果は、 $(x_1, x_2) = (2, 3)$ が最適解であり、最適値12が得られました。

目的関数fの引数はベクトルでなければならない点と、
初期値が0だとうまくいかない点(理由はよくわからない)には注意が必要です。

参考: 線形不等式制約付きの最適化関数 constrOptim - RjpWiki

2014-07-05

軸の目盛りを文字列にする

plotの軸の目盛りをいじりたい。しかしどうすれば...？
そんなときはaxisを使います。例えば以下です。

axis(side=1, at=1:10, labels=letters[1:10])

これなら、x軸の1から10番目の目盛りがa~jになります。

こんな感じに部分的にも使えます。

plot(trees$Volume, xaxt="n")
axis(1, at=11:20, labels=LETTERS[2:11])

f:id:jundoll:20140705234144p:plain

これで時間を目盛りにして見たい時も簡単。

2014-06-28

catによる出力を非表示にする

とりあえずinvisibleを使ってみるとこうなります。

> invisible(cat("I want to suppress output in R."))
I want to suppress output in R.
> # 失敗

ちなみにsuppressMessagesを使ってもダメでした。

invisible関数で非表示にできるのは、オブジェクトであるため、capture.output関数で出力を文字列として扱えるようにすれば良いようです。

> invisible(capture.output(cat("I want to suppress output in R.")))
> # 成功

2014-05-24

Rのプロキシ設定（Mac用）

研究室などのプロキシ環境下において、ネットワークに接続することができないことがあります。
そういう時は.Rprofileなどに以下のコードを入れておきましょう。

Sys.setenv("http_proxy"="http://proxy.hogehoge.jp:hogehoge2")
# hogehoge = host
# hogehoge2 = port

また、ネットワーク環境による切り替えが可能な状態であれば、以下のコードで自動的に切り替えられます。

if(system("echo `networksetup -getcurrentlocation`", intern=TRUE) == NetworkEnv.) {  
  # NetworkEnv. = "プロキシ設定のされたネットワーク環境名"
  Sys.setenv("http_proxy"="http://proxy.hogehoge.ac.jp:hogehoge2") 
  cat("=== Proxy setting changed. ===\n\n")
}