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

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

今週末の数理統計研究所でのRの集い楽しみですね. 東京は遠いのでいかないですが, ニコ生で見る予定です.
R研究集会 「データ解析環境Rの整備と利用」2012 ~午前の部~ - ニコニコ生放送
R ニコ生