1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
|
## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE,
comment = "",
fig.height = 6,
fig.width = 8,
fig.align = "center")
# out.width = "0.25\\textheight")
library(polynom)
setHook("plot.new",
list(las = function() par(las = 1),
pch = function() par(pch = 20)),
"append")
## -----------------------------------------------------------------------------
He <- list(polynomial(1), polynomial(0:1))
x <- polynomial()
for (n in 3:6) {
He[[n]] <- x * He[[n-1]] - (n-2) * He[[n-2]] ## R indices start from 1, not 0
}
He <- as.polylist(He)
plot(He)
plot(deriv(He))
plot(integral(He))
## -----------------------------------------------------------------------------
x <- c(0,1,2,4)
(op <- poly.orth(x))
(fop <- lapply(op, as.function))
(P <- sapply(fop, function(f) f(x)))
zapsmall(crossprod(P)) ### Verify orthonormality
## -----------------------------------------------------------------------------
(p1 <- poly.calc(1:6))
(p2 <- change.origin(p1, 3))
predict(p1, 0:7)
predict(p2, 0:7)
predict(p2, 0:7 - 3)
(p3 <- (p1 - 2 * p2)^2) # moderate arithmetic expression.
fp3 <- as.function(p3) # should have 1, 2, 3 as zeros
fp3(0:4)
## -----------------------------------------------------------------------------
x <- 80:89
y <- c(487, 370, 361, 313, 246, 234, 173, 128, 88, 83)
p <- poly.calc(x, y) ## leads to catastropic numerical failure!
predict(p, x) - y
p1 <- poly.calc(x - 84, y) ## changing origin fixes the problem
predict(p1, x - 84) - y
plot(p1, xlim = c(80, 89) - 84, xlab = "x - 84")
points(x - 84, y, col = "red", cex = 2)
#### Can we now write the polynomial in "raw" form?
p0 <- as.function(p1)(polynomial() - 84) ## attempt to change the origin back to zero
## leads to problems again
plot(p0, xlim = c(80, 89))
points(x, y, col = "red", cex = 2) ## major numerical errors due to finite precision
|
