2012-03-31

5. 補間へのイントロダクション：ぱらぱらめくる『A Short Course on Approximation Theory』

x軸上の異なるn個の値があって、それぞれにyの実数値が与えられるとxy平面上のn点を通るn次多項式がある
- Lagrange's 補間(Wiki)
- Newton's 補間((Wiki)
  - 傾きを考慮

# x=1,2,3を解とする多項式
p <- poly(c(1, 2, 3))
fp <- function(x) polyval(p, x)

x <- runif(4); y <- fp(x)
xx <- linspace(0, 4, 51)
yL <- lagrangeInterp(x, y, xx)
yN <- newtonInterp(x, y, xx)
## Not run: 
ezplot(fp, 0, 4)
points(xx, yL)
points(xx, yN,col=2)
## End(Not run)

- Hermite 補間(これも傾きを考慮)(Wiki)
  - 0次、1次、2次、…の傾きが合うようにする
  - pchip()関数
  - ３次までを合わせる

## Not run: 
x <- c(1, 2, 3, 4, 5, 6)
y <- c(16, 18, 21, 17, 15, 12)
plot(x, y, col="red", xlim=c(0,7), ylim=c(10,22),
     main = "Spline and 'pchip' Interpolation Polynomials")
grid()

xs <- seq(1, 6, len=51)
ys <- interp1(x, y, xs, "spline")
lines(xs, ys, col="cyan")
yp <- pchip(x, y, xs)
lines(xs, yp, col = "magenta")
## End(Not run)

ベジェ曲線・バーンスタイン多項式
- ベジェ曲線はコンピュータグラフィクスなどで点と点との間を滑らかにつなぐ曲線(Wiki)
- ベジェ曲線を引くときに大切なのがバーンスタイン多項式
- 『バーンスタイン多項式（-たこうしき、Bernstein polynomial）は、バーンスタイン基底関数 (Bernstein basis polynomials) の線形結合で与えられるバーンスタイン形式の多項式。数値的に安定』(Wikiより)
- Rのパッケージ pencopula やDPpackageにバーンスタイン多項式を使ったり求めたりする関数があるようだが、詳細不明
余談
- x軸上の異なるn個の値( $x_i;i=1,2,...,n$ )があるとき、 $(x_i)$ のすべてを通る多項式は

library(pracma)
x<-c(10,50,70)
# f(x)=0となるような多項式の係数を返す関数poly()
p<-poly(x)
# pなる係数を持つ多項式にxを代入した値
polyval(p,x) # 0になる
X<-seq(from=-1,to=101,length=10000)
plot(X,polyval(p,X),type="l")
abline(h=0)
abline(v=x)

色々な補間
- １次
  - constant: constant between points
  - linear: linear interpolation (default)
  - nearest: nearest neighbor interpolation
  - spline: cubic spline interpolation
    - スプライン曲線(Wiki)
  - cubic: cubic Hermite interpolation

library(pracma)
x <- c(0.8, 0.3, 0.1, 0.6, 0.9, 0.5, 0.2, 0.0, 0.7, 1.0, 0.4)
y <- x^2
xi <- seq(0, 1, len = 81)
yl <- interp1(x, y, xi, method = "linear")
yn <- interp1(x, y, xi, method = "nearest")
ys <- interp1(x, y, xi, method = "spline")

## Not run: 
plot(x, y); grid()
lines(xi, yl, col="blue", lwd = 2)
lines(xi, yn, col="black", lty = 2)
lines(xi, ys, col="red")

## End(Not run)

- ２次

library(pracma)
## Not run: 
    x <- linspace(-1, 1, 11)
    y <- linspace(-1, 1, 11)
    mgrid <- meshgrid(x, y)
    Z <- mgrid$X^2 + mgrid$Y^2
    xp <- yp <- linspace(-1, 1, 101)

    method <- "linear"
    zp <- interp2(x, y, Z, xp, yp, method)
    plot(xp, zp, type = "l", col = "blue")

    method = "nearest"
    zp <- interp2(x, y, Z, xp, yp, method)
    lines(xp, zp, col = "red")
    grid()
## End(Not run)

- ２次元(barycentric Lagrange)

##  Example from R-help
xn <- c(4.05, 4.10, 4.15, 4.20, 4.25, 4.30, 4.35)
yn <- c(60.0, 67.5, 75.0, 82.5, 90.0)
foo <- matrix(c(
        -137.8379, -158.8240, -165.4389, -166.4026, -166.2593,
        -152.1720, -167.3145, -171.1368, -170.9200, -170.4605,
        -162.2264, -172.5862, -174.1460, -172.9923, -172.2861,
        -168.7746, -175.2218, -174.9667, -173.0803, -172.1853,
        -172.4453, -175.7163, -174.0223, -171.5739, -170.5384,
        -173.7736, -174.4891, -171.6713, -168.8025, -167.6662,
        -173.2124, -171.8940, -168.2149, -165.0431, -163.8390),
            nrow = 7, ncol = 5, byrow = TRUE)
xf <- c(4.075, 4.1)
yf <- c(63.75, 67.25)
barylag2d(foo, xn, yn, xf, yf)
#  -156.7964 -163.1753
#  -161.7495 -167.0424

# Find the minimum of the underlying function
bar <- function(xy) barylag2d(foo, xn, yn, xy[1], xy[2])
optim(c(4.25, 67.5), bar)  # "Nelder-Mead"
# $par
# 4.230547 68.522747
# $value
# -175.7959

## Not run: 
# Image and contour plots
image(xn, yn, foo)
contour(xn, yn, foo, col="white", add = TRUE)
xs <- seq(4.05, 4.35, length.out = 51)
ys <- seq(60.0, 90.0, length.out = 51)
zz <- barylag2d(foo, xn, yn, xs, ys)
contour(xs, ys, zz, nlevels = 20, add = TRUE)
contour(xs, ys, zz, levels=c(-175, -175.5), add = TRUE)
points(4.23, 68.52)
## End(Not run)