データをプロットしたグラフに、推定された予測値や、推定されたパラメーターを数式ごと書き込みたいときはありませんか？
そう難しい話では無いのですが、幾つかコツがあるので例を作ってみました。
1. データセット
グラフを描くにはデータセットが要るので、今回は以下の数字を
# 奈良女子大学のページにあった実験観測データを拝借
t <- c(0,5,10,15,20,25,30) # 時間
Ts <- c(56.5, 32.1, 21.9, 17.0, 14.6, 11.9, 10.9) # 物体の温度
Tm <- 0 # 外気温（冷却材は氷でゼロ度だが、室温はもっと高い模様。なお、9.660583度を仮定すると、もっとも当てはまりが良い）
2. 推定を行う
推定を行い、予測値を計算してみましょう。モデルの制約によりpredict関数は使えません。
# ニュートンの冷却の法則を想定し、定数を推定する
dT <- (Ts - Tm)/(Ts[1]- Tm)
# OLSで推定。
r <- lm(log(dT) ~ t + 0)
# 予測値を計算
prdct <- Tm + (Ts[1] - Tm)*exp(coef(r)[1]*t)
4. データをプロットする
まずは観測値をプロットしてみます。
par(mar=c(5, 3.5, 1, 1), mgp=c(2, 1, 0))
plot(Ts ~ t, xlab="経過時間", ylab="", axes=FALSE, xlim=xlim, ylim=ylim)
x_at <- seq(xlim[1], xlim[2], 10)
y_at <- seq(ylim[1], ylim[2], 25)
x_labels <- sprintf("%d分", x_at)
y_labels <- sprintf("%d℃", y_at)
axis(1, at=x_at, labels=x_labels, pos=ylim[1])
axis(2, at=y_at, labels=y_labels, pos=xlim[1], las=2)
5. 推定値を重ね書きする
推定値をグラフに書き込みます。ここはRの紹介で良く触れられますね。
par(new=T)
plot(prdct ~ t, type="l", xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim, lty=2)
なお、こんな手間隙かけずにlines(t, prdct, lty=2)するべきだろと突っ込まれました。おっしゃる通り。
直線なのが気に入らない人は、plotの代わりに以下のようにcurveを使ってください。
par(new=T)
curve(Tm + (Ts[1] - Tm)*exp(coef(r)[1]*t), xname="t", xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim, lty=2)
6. 推定された式を凡例として表示する
まず、数式をsprintf関数で作って、parse関数で表現式に変えます。expression関数は、代入が手間なので使いません。使える数式の記号はhelp plotmathをすると一覧が見られます。
expr <- parse(text = sprintf("%.3f + (%.3f-%.3f)*e^{%1.3f*t}", Tm, Ts[1], Tm, coef(r)[1]))
次に、legend関数で凡例を表示します。推定値はlty、観測値はpchになりますが、使わない方には-1を入れていることに注意してください。
legend("topright", lty=c(2, -1), pch=c(-1, 1), legend=c(expr, "観測値"), bty="n")
7. データの出所を明記する
データを拝借しているので、出所を明記しておきましょう。
mtext("出所) 奈良女子大 - 「第11回の授業 (07/07/09)」の実験データから編集" , side=1, line=3.5, adj=0)
8. 完成図
相関係数などを書いても良かったかも知れません。
9. 補足
文字が化ける場合は、以下のようにpar()関数で指定するグラフィックス・パラメーターでフォントを指定してください。フォントはMeiryoを指定していますが、お使いのOSにあわせましょう。
windowsFonts(GraphFont = windowsFont("Meiryo"))
par(mar=c(5, 3.5, 1, 1), mgp=c(2, 1, 0), family="GraphFont")
画面のグラフを保存したい場合は、以下のような感じで。
dev.copy(png, "ニュートンの冷却の法則の実験結果.png", width=640, height=480, bg="white")
dev.off() # 画面のグラフを消す
dev.off() # PNGファイルを閉じる
*1:推定値のプロットにpar(new=T)とplot()を使うからズレるのであって、低水準関数のlines()を使えば座標の範囲を決めなくてもいいと言うツッコミがありましたが、ここを決めておかないと軸目盛などの位置もズレると言うか、あわせるのが困難になるので、やはり決めておいた方がいいです。
