Rで自作のラグランジェ補間法を

橋本愛さんに踵落ししてもらいたい系男子です.

あっどーも, 僕です.

こないだLagrange補間多項式のレポートがでました. 授業ではエクセルで教わりました. ですが, ご存知の方もいらっしゃると思いますが Lagrange補間多項式をエクセルでやるのは データ点が少なくても結構たいへんです.

実際レポートでは7点だったのですが正直しんどかったです. なのでエクセルでは諦めてRでやってみたら結構簡単にできたので 載せてみたいと思います.

N+1個の点, を通過する 曲線を表すLagrangeの補間多項式はこちらです.

また, Lagrangeの補間多項式をn+1次の多項式

を用いて表すと次式となります.

Rのコードはこちら.

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

#補間多項式の要素
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))

#結果
x <- seq(from = 0.5, to = 6, by = 0.1)
y <- sapply(x, Lag)

#グラフづくり
def.par <- par(no.readonly = TRUE) 
par(mai = c(1.5, 1.5, 0.5,0.5))
curve(sin, 0.5, 6, col = "#f39800", lwd = 3,
      bg = "white", cex = 1.3, cex.lab = 1.5, cex.axis = 1.3)
points(x, y, pch = 19, col = "#0068b7", cex = 1.3)
legend("bottomleft",  lwd = 2, col = c("#f39800", "#0068b7"),
       legend = c(expression(sin(x)), "Lagrange"), 
       bty = "n", cex = 1.3, lty = c(1, -1), pch = c(-1, 19),
       pt.cex = 2)
par(def.par)　　　　　　　　　　　

結果がこちらです. データ点が多いと役立たずって聞きましたが, これだと なかなか良い感じにできているように見えます. つーかこれ合っているのかな？ まぁもうレポート出してしまったからいいんだけどさ.

後でRでのスプラインの仕方を調べたら竹澤 先生の著書がわかりやすかったです.

Rによるノンパラメトリック回帰の入門講義

追記

ラグランジュ補間法のパッケージ調べました. こちらです.