R语言与函数估计学习笔记（函数展开）

函数估计

说到函数的估计我们可以肯定的一点是我们很难得到原模型的函数，不过我们可以找到一个不坏的函数去逼近它，所以我们的函数估计从函数展开开始说起。

函数展开

Taylor展开

首先不得不提的就是大名鼎鼎的Taylor展开，它告诉我们一个光滑的函数在x=t的一个邻域内有Taylor展式

它给我们的一个重要启示就是我们可以把我们感兴趣的函数拆解成若干个简单函数q0(x),q1(x)⋯,的线性组合。

那么还剩一个问题，就是qj(x)选什么。当然一个简单的选择就是qj(x)=xj,或者我们取t=x¯,qj(x)=(x−x¯)j。我们来看看这组函数基qj(x)=xj对标准正态密度函数的估计效果

x <- seq(-3, 3, by = 0.1)
y <- dnorm(x)
model <- lm(y ~ poly(x, 2))
plot(y, type = "l")
lines(fitted(model), col = 2)

summary(model)

##
## Call:
## lm(formula = y ~ poly(x, 2))
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -0.07901 -0.06035 -0.00363  0.05864  0.10760
##
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.64e-01   8.09e-03    20.2   <2e-16 ***
## poly(x, 2)1 -1.77e-16   6.32e-02     0.0        1    
## poly(x, 2)2 -9.79e-01   6.32e-02   -15.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0632 on 58 degrees of freedom
## Multiple R-squared:  0.805,  Adjusted R-squared:  0.799
## F-statistic:  120 on 2 and 58 DF,  p-value: <2e-16

从图像上来看，这个拟合不是很好，我们可以认为是p较小造成的，一个解决办法就是提高p的阶数，令p=10我们可以试试：

model1 <- lm(y ~ poly(x, 10))
x <- seq(-3, 3, by = 0.1)
y <- dnorm(x)
model <- lm(y ~ poly(x, 2))
plot(y, type = "l")
lines(fitted(model), col = 2)
lines(fitted(model1), col = 3)

summary(model1)

##
## Call:
## lm(formula = y ~ poly(x, 10))
##
## Residuals:
##       Min        1Q    Median        3Q       Max
## -3.86e-04 -2.03e-04  1.45e-05  1.83e-04  2.83e-04
##
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.64e-01   2.94e-05  5572.4   <2e-16 ***
## poly(x, 10)1  -1.92e-16   2.29e-04     0.0        1    
## poly(x, 10)2  -9.79e-01   2.29e-04 -4268.8   <2e-16 ***
## poly(x, 10)3   2.36e-16   2.29e-04     0.0        1    
## poly(x, 10)4   4.54e-01   2.29e-04  1979.0   <2e-16 ***
## poly(x, 10)5  -1.65e-16   2.29e-04     0.0        1    
## poly(x, 10)6  -1.54e-01   2.29e-04  -672.4   <2e-16 ***
## poly(x, 10)7   1.67e-17   2.29e-04     0.0        1    
## poly(x, 10)8   4.09e-02   2.29e-04   178.5   <2e-16 ***
## poly(x, 10)9   2.07e-16   2.29e-04     0.0        1    
## poly(x, 10)10 -8.85e-03   2.29e-04   -38.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.000229 on 50 degrees of freedom
## Multiple R-squared:     1,   Adjusted R-squared:     1
## F-statistic: 2.26e+06 on 10 and 50 DF,  p-value: <2e-16

从上图看到，拟合效果好了不少，这样看上去我们只需要提高基函数阶数就可以解决拟合优度的问题了。但是注意到随着阶数提高，可能出现设计阵降秩的情形，也有可能出现复共线性，这是我们不希望看到的。为了解决第一个问题，我们的做法是限制p的最大取值，如将p限制在5以下；对于第二个问题，我们的做法便是采用正交多项式基。

正交多项式展开

正交多项式的相关定义可以参阅wiki，这里就不在啰嗦了，我们这里列出Legendre多项式基与Hermite多项式基。
其中Legendre多项式基已经在wiki中给出了，其取值范围是[-1,1]，权函数是1，表达式为：

Legendre多项式基的递归表达式可以表达为：

Hermite多项式基的取值范围为(−∞,∞),对应的权函数是标准正态分布的核函数，表达式为：

Hermite多项式基的递归表达式可以表达为:

我们这里来看一个例子，假设真实模型为y=5xcos(5πx),我们一共做了10次试验，得到了10个观测，现在我们要找一个拟模型来近似这个真实模型。我们来看看多项式基的效果：

x <- seq(-1, 1, length = 20)
y <- 5 * x * cos(5 * pi * x)
f <- function(x) 5 * x * cos(5 * pi * x)
curve(f, -1, 1)
points(x, y)
A <- data.frame(x = seq(-1, 1, length = 1000))
model.linear <- lm(y ~ poly(x, 6))
lines(seq(-1, 1, length = 1000), predict(model.linear, A), col = 2)
model.linear1 <- lm(y ~ poly(x, 9))
lines(seq(-1, 1, length = 1000), predict(model.linear1, A), col = 3)
z <- matrix(rep(NA, 6 * length(x)), length(x), 6)
z[, 1] <- x
z[, 2] <- (3 * x^2 - 1)/2
z[, 3] <- (5 * x^3 - 3 * x)/2
z[, 4] <- (35 * x^4 - 30 * x^2 + 3)/8
z[, 5] <- (2 * 5 - 1)/5 * x * z[, 4] - 0.8 * z[, 3]
z[, 6] <- (2 * 6 - 1)/6 * x * z[, 5] - 5/6 * z[, 4]
model.linear2 <- lm(y ~ z)
x <- seq(-1, 1, len = 1000)
z <- matrix(rep(NA, 6 * length(x)), length(x), 6)
z[, 1] <- x
z[, 2] <- (3 * x^2 - 1)/2
z[, 3] <- (5 * x^3 - 3 * x)/2
z[, 4] <- (35 * x^4 - 30 * x^2 + 3)/8
z[, 5] <- (2 * 5 - 1)/5 * x * z[, 4] - 0.8 * z[, 3]
z[, 6] <- (2 * 6 - 1)/6 * x * z[, 5] - 5/6 * z[, 4]
B <- as.data.frame(z)
lines(x, predict(model.linear2, B), col = 4)
letters <- c("orignal model", "6 order poly-reg", "9 order poly-reg", "6 order orth-reg")
legend("bottomright", legend = letters, lty = 1, col = 1:4, cex = 0.5)

Fourier展开

这里我们就可以看到，多项式拟合对于这种含周期的问题的解决效果是很不好的，正交多项式完全不行，可见问题并不是出在复共线性上，对于含周期的函数的逼近我们可以引入Fourier基：

我们来看看拟合效果：

x <- seq(-1, 1, length = 10)
y <- 5 * x * cos(5 * pi * x)
f <- function(x) 5 * x * cos(5 * pi * x)
curve(f, -1, 1, ylim = c(-15.5, 15.5))
points(x, y)
model.linear <- lm(y ~ poly(x, 7))
A <- data.frame(x = seq(-1, 1, length = 1000))
lines(seq(-1, 1, len = 1000), predict(model.linear, A), col = 2)
model.linear1 <- lm(y ~ poly(x, 9))
lines(seq(-1, 1, len = 1000), predict(model.linear1, A), col = 3)
z <- matrix(rep(NA, 6 * length(x)), length(x), 6)
z[, 1] <- cos(2 * pi * x)
z[, 2] <- sin(2 * pi * x)
z[, 3] <- cos(4 * pi * x)
z[, 4] <- sin(4 * pi * x)
z[, 5] <- cos(6 * pi * x)
z[, 6] <- sin(6 * pi * x)
model.linear2 <- lm(y ~ z)
x <- seq(-1, 1, len = 1000)
z <- matrix(rep(NA, 6 * length(x)), length(x), 6)
z[, 1] <- cos(2 * pi * x)
z[, 2] <- sin(2 * pi * x)
z[, 3] <- cos(4 * pi * x)
z[, 4] <- sin(4 * pi * x)
z[, 5] <- cos(6 * pi * x)
z[, 6] <- sin(6 * pi * x)
B <- as.data.frame(z)
lines(x, predict(model.linear2, B), col = 4)
letters <- c("orignal model", "7 order poly-reg", "9 order poly-reg", "Fourier-reg")
legend("bottomright", legend = letters, lty = 1, col = 1:4, cex = 0.5)

可见Fourier基对周期函数的拟合还是很好的。但是这必须是不含趋势的结果，含趋势的只能在局部有个不错的拟合，如果我们把上面的模型换为5x+cos(5πx),可以看到Fourier基拟合的效果是十分糟糕的。

x <- seq(-1, 1, length = 10)
y <- 5 * x + cos(5 * pi * x)
f <- function(x) 5 * x + cos(5 * pi * x)
curve(f, -1, 1)
points(x, y)
model.linear <- lm(y ~ poly(x, 7))
A <- data.frame(x = seq(-1, 1, length = 1000))
lines(seq(-1, 1, len = 1000), predict(model.linear, A), col = 2)
model.linear1 <- lm(y ~ poly(x, 9))
lines(seq(-1, 1, len = 1000), predict(model.linear1, A), col = 3)
z <- matrix(rep(NA, 6 * length(x)), length(x), 6)
z[, 1] <- cos(2 * pi * x)
z[, 2] <- sin(2 * pi * x)
z[, 3] <- cos(4 * pi * x)
z[, 4] <- sin(4 * pi * x)
z[, 5] <- cos(6 * pi * x)
z[, 6] <- sin(6 * pi * x)
model.linear2 <- lm(y ~ z)
x <- seq(-1, 1, len = 1000)
z <- matrix(rep(NA, 6 * length(x)), length(x), 6)
z[, 1] <- cos(2 * pi * x)
z[, 2] <- sin(2 * pi * x)
z[, 3] <- cos(4 * pi * x)
z[, 4] <- sin(4 * pi * x)
z[, 5] <- cos(6 * pi * x)
z[, 6] <- sin(6 * pi * x)
B <- as.data.frame(z)
lines(x, predict(model.linear2, B), col = 4)
letters <- c("orignal model", "7 order poly-reg", "9 order poly-reg", "Fourier-reg")
legend("bottomright", legend = letters, lty = 1, col = 1:4, cex = 0.5)

样条基展开

有些时候我们对全局的拟合是有缺陷的，所以可以进行分段的拟合，一旦确定了分段的临界点，我们就可以进行局部的回归，样条基本上就借鉴了这样一个思想。
为了增加局部的拟合优度，我们在原来的函数基1,x,x2,⋯,xp上加上其中，knot表示节点，函数(x−knoti)+表示函数(x−knoti)取值为正时取函数值，否则取0.

x <- seq(-1, 1, length = 20)
y <- 5 * x * cos(5 * pi * x)
f <- function(x) 5 * x * cos(5 * pi * x)
curve(f, -1, 1)
points(x, y)
model.reg <- lm(y ~ poly(x, 5))
A <- data.frame(x = seq(-1, 1, length = 1000))
lines(seq(-1, 1, len = 1000), predict(model.reg, A), col = 2)
ndat <- length(x)
knots <- seq(-1, 1, length = 10)
f <- function(x, y) ifelse(y > x, (y - x)^3, 0)
X <- matrix(rep(NA, length(x) * (3 + length(knots))), length(x), (3 + length(knots)))
for (i in 1:3) X[, i] <- x^i
for (i in 4:(length(knots) + 3)) X[, i] <- f(knots[(i - 3)], x)
model.cubic <- lm(y ~ X)
x <- seq(-1, 1, length = 1000)
X <- matrix(rep(NA, length(x) * (3 + length(knots))), length(x), (3 + length(knots)))
for (i in 1:3) X[, i] <- x^i
for (i in 4:(length(knots) + 3)) X[, i] <- f(knots[(i - 3)], x)
A <- as.data.frame(X)
lines(seq(-1, 1, len = 1000), predict(model.cubic, A), col = 3)

从上图中我们可以看到加上样条基后，拟合效果瞬间提高了不少，三阶样条基就可以匹敌5~6阶的多项式基了。R中的splines包中提供了polyspline函数，来做样条拟合，我们可以看看在这个例子中它几乎就是原函数的“复制”。

x <- seq(-1, 1, length = 20)
y <- 5 * x * cos(5 * pi * x)
library(splines)
model <- polySpline(interpSpline(y ~ x))
# print(model)
plot(model, col = 2)
f <- function(x) 5 * x * cos(5 * pi * x)
curve(f, -1, 1, ylim = c(-15.5, 15.5), add = T)
points(x, y)

本节的最后，我们最后来看看函数展开的相关内容，如果说我们已经知道了函数f(x)的表达式，想求解一个近似的函数展开式的系数，我们只需要将f(x)拆解为f(x)=g(x)p(x),其中p(x)为密度函数，那么展开式系数可以近似的表示为其中x1,⋯,xn是由p(x)产生的随机数。