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

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rのグラフにスプライン曲線を描きこむ

TTFなどのフォントもそうですし、CADなどグラフィックス分野の方が使う事が多いと思いますが、間欠的なデータから連続した線を描画したいときは実用上は多々あります。ラグランジュ補間、ニュートン補間などバリエーションは幾つもありますが、Rではスプライン補間をする事が容易です。
R でグラフ中にプロットされてない値を予測してみる」を見ていたら簡単そうなので、試してみました。y=\frac{\sin(x)}{x}と言う簡単そうで積分するとトリッキーで大変なルベーグ積分不可の関数の極値をつないでみます。
まずは、y=\frac{\sin(x)}{x}を描いてみます。

次に、極値を求めます。

最後に、極大値と極小値をスプライン曲線で結びます。

意味が良く分からないものの、それらしいグラフになりました。ソースコードは以下になります。

# 関数を定義
f <- function(x){
  sin(x)/x
}


# fの一階微分
df <- function(x){
  (cos(x)*x-sin(x))/x^2
}


# fの二階微分
ddf <- function(x){
  (-sin(x)*x*2*x - (cos(x)*x-sin(x))*2*x)/x^4
}


# xの変分
dx <- 0.01
# xを生成(始点を僅かに正にすることでyが∞を回避)
x <- seq(1e-6, 30, dx)
# yを計算
y <- f(x)


getExtremum <- function(x){
  dy <- df(x)
  n <- length(dy)
  # 微分係数dy[i]*dy[i+1]がゼロ以下の場合は符号が逆転しており極値近傍と見なせる
  x[dy[-n]*dy[-1] <= 0]
}
extremum <- getExtremum(x)


# ddf(x[i])が負ならば、極大値
ulim <- c(extremum[ddf(extremum)<0])
# ddf(x[i])が正ならば、極小値
llim <- c(extremum[ddf(extremum)>0])


# スプライン関数を導出
sp_u <- smooth.spline(ulim, f(ulim))
sp_l <- smooth.spline(llim, f(llim))


# xとyのそれぞれ上下限
ylim <- c(-0.5, 1)
xlim <- c(min(x), max(x))


# 余白位置などを調整
par(mgp=c(1, 1, 0), mai=c(1,1,1,1), mar=c(4,3,3,1))


# 元データをプロット
curve(f(x), xlim=xlim, ylim=ylim, type="l", xlab="x", ylab=expression(sin(x)/x), axes=FALSE, main="スプライン関数のお試し")


# スプラインを描画/predict()関数で補間点を計算
par(new=T); plot(predict(sp_u, x), xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="", ylab="", axes=FALSE, col="red2")
par(new=T); plot(predict(sp_l, x), xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="", ylab="", axes=FALSE, col="red2")


# スプライン関数の元もプロットしておく
par(new=T); plot(ulim, f(ulim), xlim=xlim, ylim=ylim, type="p", lty="dotted", xlab="", ylab="", axes=FALSE, bg="gray")
par(new=T); plot(llim, f(llim), xlim=xlim, ylim=ylim, type="p", lty="dotted", xlab="", ylab="", axes=FALSE, bg="gray")


# 軸を描画
axis(1, at=seq(0, 30, 10), pos=ylim[1])
axis(2, at=c(seq(ylim[1], ylim[2], 0.25)), pos=0)

極大値と極小値を計算するほうが手間がかかっていますね。