RでLagrange polynomial interpolation - サボタージュ禁止のおさぼり日記

私はまだ知らない. 尊敬している人の鼻毛がでているときどう伝えればいいのか.

あっ, どーも僕です.

Rでラグランジュ補間多項式

研究の息抜きにラグランジュ補間多項式のRでの仕方を調べてみました. まぁ息抜きというか現実逃避ですかね.

・・・

ラグランジュ補間多項式自体は簡単に説明してあるのでそれをどうぞ.
Rで補間法 - サボタージュ禁止のおさぼり日記

このときは, レポートのためにラグランジュ補間を計算したので, 自分で関数を組みました. しかし, Rならきっとパッケージがあるだろうし, 自作関数の検証もしたいということで今回探してみました.

まぁ探すといってもsosパッケージを使ったらであっというまでした( ｀・∀・´)ﾉ
こちらが, "Lagrange poly"で検索した結果です.
f:id:aaaazzzz036:20121129191646p:plain
何個かありそうですが, とりあえずpracmaパッケージのlagrangeInterpを使ってみることにしました.

使い方は,

library(pracma)
# X, Y : データ点
# x : 内挿したい点
lagrangeInterp(X, Y, x)

だけでとっても簡単ですね.

これを利用していつか自分が作ったものと比べてみました.
結果がこちら.
f:id:aaaazzzz036:20121129193659j:plain

自作関数の結果とパッケージでのぱっと見の結果が同じでバッチリですね！値はどうだろうということでこちら.

#検証
sum(round(abs(MyLag-PracmaLag), 6))
[1] 0

バッチリですね.

最後に, 前とほとんど同じですが一応スクリプトはこちらです.

#************************************************************
#ラグランジュ補間補間多項式によるsin(x)の近似
#データ
X   <- 0.5 + 0:6
Y   <- sin(X)

#内挿したい点
x <- seq(from = 0.5, to = 6.5, by = 0.1)

#オリジナル
#補間多項式の要素
lag <- function(i, x)  prod(x - X[-i])/prod(X[i]-X[-i])*Y[i]
#補間多項式
Lag <- function(x) sum(sapply(1:length(X), lag, x))
#結果
MyLag <- sapply(x, Lag)

#パッケージ
library(pracma)
PracmaLag <- lagrangeInterp(X, Y, x)

#グラフづくり
par(mai = c(1.5, 1.5, 0.5,0.5))
curve(sin, 0.5, 6.5, col = "#f39800", lwd = 3, type = "n",
      bg = "white", tcl = 0.5,
      cex = 1.3, cex.lab = 1.5, cex.axis = 1.3)
points(x, MyLag, pch = 19, col = "#0068b7", cex = 2.5)
points(x, PracmaLag, pch = 19, col = "#f39800", cex = 1.5)
legend("bottomleft",  lwd = 2, col = c("#0068b7","#f39800"),
       legend = c("MyLagrange", "pracmaLagrange"), 
       bty = "n", cex = 2, lty = c(-1, -1), pch = c(19, 19),
       pt.cex = 3)

#検証
sum(round(abs(MyLag-PracmaLag), 6))
[1] 0