Rのggplot2でグラフを作りたい

Rおじさんだよ

描画が綺麗にできるRパッケージ, ggplot2を使いたいです.
紹介のために, ローレンツ方程式の数値解を描画します.
ローレンツ方程式とは.
http://en.wikipedia.org/wiki/Lorenz_system
ローレンツ方程式はカオス的な振る舞いをする有名な問題です.
\begin{align*}
\frac{\mathrm{d}x}{\mathrm{dt}} &= ax + by \\
\frac{\mathrm{d}y}{\mathrm{dt}} &= b(y-z) \\
\frac{\mathrm{d}z}{\mathrm{dt}} &= -xy + cy -z
\end{align*}
パラメータを神調整すればカオスアトラクタが見られます.

Rで常微分方程式を解く

RのdeSolveパッケージは, 適当な微分方程式を指定した解法で解いてくれます.
Rスクリプトなので遅いです. 簡単な問題ならこれで十分ですが, 実用は無理かと.
この話はまた今度.

# install.packages("deSolve")
library(deSolve) 

# deSolve用に解きたい微分方程式を書きます
# ローレンツ方程式を書きました
Lorenz<-function(t, state, parameters) {
    with(as.list(c(state, parameters)),{
        # 微分方程式を書く (xの時間微分) = f(x)の形
        dX <- a*X + Y*Z
        dY <- b * (Y-Z)
        dZ <- -X*Y + c*Y - Z   
        # listに入れる 
        list(c(dX, dY, dZ))
        }
    ) 
}

deSolve用に微分方程式を書きました.
あとは条件を指定すれば勝手に解いてくれます.

# 時間と初期値, パラメータを決めます
parameters <- c(a = -8/3, b = -10, c = 28) 
# こうするとストレンジアトラクタです
state <- c(X = 1, Y = 1, Z = 1)
times <- seq(0, 100, by = 0.001)


# 常微分方程式を解いてくれる関数odeにお願いします
# odeの引数は...見た目通り.

out <- ode(y = state, times = times, parms = parameters,
           func = Lorenz, method = "rk4")

# method = "" で解法を指定できます.
# デフォルトではlsodeが選択されます.

outに常微分方程式を解いた結果が入っています.
スキームは常微分方程式の特徴にあったものを選びます.
問題がstiffがどうか, などに気をつけます. 授業でやりましたね.

デフォルトで指定されるのはlsodeです.
stiffであるかどうかを反映してステップ幅や解法を勝手に決めます.
何も分からなければデフォルトでいいはずです.
可変ステップなので, ステップ幅を固定したい場合はRK4などに指定します.

次に描画です.

ggplot2による描画

Rにはデフォルトのplot関数とそのお友達がいます.
機能は充実しているのですが, 見栄えが悪いです.
ダサいので, ちゃんとしたグラフを描くときにはMatlabなどを使っていました.
しかしRでもかっこ良くプロットできます.
ggplot2をインストールしましょう.

# install.packages("ggplot2")
library(ggplot2)

out <- data.frame("time"=out[,1], "X" = out[,2], "Y" = out[,3])

# ふつうのplot
pdf("dasaiplot.pd")
plot(out[,2],out[,3], "l", 
     xlab="X", ylab="Y", main = "Lorentz Equation",
     col = "blue") 
dev.off()

ふつうのplotをしました.
f:id:gg_hatano:20130704202640p:plain
ダサいですね.

# ggplotを使ったplot
pdf("Kakkoiiplot.pdf")
p <- ggplot(out, aes(x=X, y=Y)) +   # plotに使うdataframeを指定
      geom_path(aes(color=time)) +  # 点の移動を描画したいときにはgeom_path
      ggtitle("Lorentz Equation")   # タイトルを指定
plot(p)
dev.off()