R語(yǔ)言與函數(shù)估計(jì)學(xué)習(xí)筆記（函數(shù)展開(kāi)）

函數(shù)估計(jì)

說(shuō)到函數(shù)的估計(jì)我們可以肯定的一點(diǎn)是我們很難得到原模型的函數(shù)，不過(guò)我們可以找到一個(gè)不壞的函數(shù)去逼近它，所以我們的函數(shù)估計(jì)從函數(shù)展開(kāi)開(kāi)始說(shuō)起。

函數(shù)展開(kāi)

Taylor展開(kāi)

首先不得不提的就是大名鼎鼎的Taylor展開(kāi)，它告訴我們一個(gè)光滑的函數(shù)在x=t的一個(gè)鄰域內(nèi)有Taylor展式

它給我們的一個(gè)重要啟示就是我們可以把我們感興趣的函數(shù)拆解成若干個(gè)簡(jiǎn)單函數(shù)q0(x),q1(x)?,的線性組合。

那么還剩一個(gè)問(wèn)題，就是qj(x)選什么。當(dāng)然一個(gè)簡(jiǎn)單的選擇就是qj(x)=xj,或者我們?nèi)=xˉ,qj(x)=(x?xˉ)j。我們來(lái)看看這組函數(shù)基qj(x)=xj對(duì)標(biāo)準(zhǔn)正態(tài)密度函數(shù)的估計(jì)效果

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

從圖像上來(lái)看，這個(gè)擬合不是很好，我們可以認(rèn)為是p較小造成的，一個(gè)解決辦法就是提高p的階數(shù)，令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

從上圖看到，擬合效果好了不少，這樣看上去我們只需要提高基函數(shù)階數(shù)就可以解決擬合優(yōu)度的問(wèn)題了。但是注意到隨著階數(shù)提高，可能出現(xiàn)設(shè)計(jì)陣降秩的情形，也有可能出現(xiàn)復(fù)共線性，這是我們不希望看到的。為了解決第一個(gè)問(wèn)題，我們的做法是限制p的最大取值，如將p限制在5以下；對(duì)于第二個(gè)問(wèn)題，我們的做法便是采用正交多項(xiàng)式基。

正交多項(xiàng)式展開(kāi)

正交多項(xiàng)式的相關(guān)定義可以參閱wiki，這里就不在啰嗦了，我們這里列出Legendre多項(xiàng)式基與Hermite多項(xiàng)式基。
其中Legendre多項(xiàng)式基已經(jīng)在wiki中給出了，其取值范圍是[-1,1]，權(quán)函數(shù)是1，表達(dá)式為：

Legendre多項(xiàng)式基的遞歸表達(dá)式可以表達(dá)為：

Hermite多項(xiàng)式基的取值范圍為(?∞,∞),對(duì)應(yīng)的權(quán)函數(shù)是標(biāo)準(zhǔn)正態(tài)分布的核函數(shù)，表達(dá)式為：

Hermite多項(xiàng)式基的遞歸表達(dá)式可以表達(dá)為:

我們這里來(lái)看一個(gè)例子，假設(shè)真實(shí)模型為y=5xcos(5πx),我們一共做了10次試驗(yàn)，得到了10個(gè)觀測(cè)，現(xiàn)在我們要找一個(gè)擬模型來(lái)近似這個(gè)真實(shí)模型。我們來(lái)看看多項(xiàng)式基的效果：

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展開(kāi)

這里我們就可以看到，多項(xiàng)式擬合對(duì)于這種含周期的問(wèn)題的解決效果是很不好的，正交多項(xiàng)式完全不行，可見(jiàn)問(wèn)題并不是出在復(fù)共線性上，對(duì)于含周期的函數(shù)的逼近我們可以引入Fourier基：

我們來(lái)看看擬合效果：

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)

可見(jiàn)Fourier基對(duì)周期函數(shù)的擬合還是很好的。但是這必須是不含趨勢(shì)的結(jié)果，含趨勢(shì)的只能在局部有個(gè)不錯(cuò)的擬合，如果我們把上面的模型換為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)

樣條基展開(kāi)

有些時(shí)候我們對(duì)全局的擬合是有缺陷的，所以可以進(jìn)行分段的擬合，一旦確定了分段的臨界點(diǎn)，我們就可以進(jìn)行局部的回歸，樣條基本上就借鑒了這樣一個(gè)思想。
為了增加局部的擬合優(yōu)度，我們?cè)谠瓉?lái)的函數(shù)基1,x,x2,?,xp上加上其中，knot表示節(jié)點(diǎn)，函數(shù)(x?knoti)+表示函數(shù)(x?knoti)取值為正時(shí)取函數(shù)值，否則取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階的多項(xiàng)式基了。R中的splines包中提供了polyspline函數(shù)，來(lái)做樣條擬合，我們可以看看在這個(gè)例子中它幾乎就是原函數(shù)的“復(fù)制”。

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)

本節(jié)的最后，我們最后來(lái)看看函數(shù)展開(kāi)的相關(guān)內(nèi)容，如果說(shuō)我們已經(jīng)知道了函數(shù)f(x)的表達(dá)式，想求解一個(gè)近似的函數(shù)展開(kāi)式的系數(shù)，我們只需要將f(x)拆解為f(x)=g(x)p(x),其中p(x)為密度函數(shù)，那么展開(kāi)式系數(shù)可以近似的表示為其中x1,?,xn是由p(x)產(chǎn)生的隨機(jī)數(shù)。