読者です 読者をやめる 読者になる 読者になる

J's blog

趣味で統計•データ解析をしています

RでLPを解く

R

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

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

R

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

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

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

R

とりあえず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.")))
> # 成功

Rのプロキシ設定(Mac用)

R

研究室などのプロキシ環境下において、ネットワークに接続することができないことがあります。
そういう時は.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")
}


参考:Macのネットワーク環境に合わせてHTTP_PROXYを切り替えるシェルスクリプト - Qiita

manipulate関数によるインタラクティブな描画

R

メモランダム.


RStudio導入時におそらく自動でインストールされている(?),manipulateパッケージのmanipulate関数を使用することでインタラクティブな図が描ける.

使い方はここ→http://www.rstudio.com/ide/docs/advanced/manipulateを見ればよ〜〜くわかる.

ちなみにggplotでもちゃんと使える.



その他参考:http://geocities.yahoo.co.jp/gl/snatool/view/20110616/1308176519