征罚性回归zz
1.主成分回归（Principal Component Regression，PCR）
主成分回归是预测变量在主成分上的回归。你可以手动地计算它：
#计算特征值和和特征向量
e <- eigen(cor(x))
#计算主成分
enx <- scale(x) %*% e$vectors
#在新的向量（主成分）上做回归
lm(y ~ enx)
或者更直接计算：
lm(y~princomp(x)$scores)
如果只想要前三个主成分作回归：
lm(y~princomp(x)$scores[,1:3])
解释新的坐标（主成分）很麻烦，但是命今“biplot(princomp(x))”可以帮上忙：
你也可以用独立成分分析（Independant Component Analysis）代替主成分分析(Principal Component Analysis ,PCA),得到独立成分回归（“Independant Component Regression“）。如果你的数据不服从高斯分布，可以参考《Independent variable selection: Application of independent component analysis to forecasting a stock index》
然而，这儿还是存在一点问题：当我们在选择新的基的时候，我们仅仅考虑了自变量x,并没有考虑因变量y，这样一些自变量在PCA（主成分分析）中表面上是无关的，但它有可能和y联系得较紧密；一些自变量在PCA中是很关键的，但它可能和因变量无关。
2.偏最小二乘（Partial Least Squares,PLS)
偏最小二乘主要用来处理这类问题：它是在新基上的回归，类似考虑了因变量的主成分分析：
library(help=pls)
library(pls)
?simpls.fit
在互联网上，你可以从NIPALS（Non-Linearly Iterated Partial Least Squares,非线性偏最小二乘迭代算法）中找到如下的PCA（主成分分析）的算法：
(t=scores, p=loadings)
begin loop
select a proxy t-vector (any large column of X)
small = 10^(-6) (e.g.)
begin loop
p = X'*t
p is normalised to length one: abs(p)=1
t = X*p
end loop if convergence: (t_new - t_old) << small
remove the modelled information from X: X_new = X_old - t*p'
end loop (if all PCs are calculated)
http://www.bitjungle.com/~mvartools/muvanl/
你可以认为PLS是介于经典的多元回归和主成分回归之间，如果经典的线性回归是“0”，主成分回归是“1”，那么PLS就是“0.5”。
PCR和PLS之间的不同，就是因子得分的计算方法不同
PCR：因子从X'X中得出
PLS：因子从Y'XX'Y中得出
下面是PLS包的下载地址：
http://cran.r-project.org/src/contrib/Devel/ http://www.gsm.uci.edu/~joelwest/SEM/PLS.html
library(help="pls.pcr")
3.惩罚最小二乘（Penalized Least Squares）
最小二乘就是构造了一个函数f，并使下式最小化：
但是这样容易产生病态函数。为了使其在一个“好”的范围内求解，我们可以给它加上一项
这里的J(f)增加了函数了复杂度（naughtiness），λ是一个由用户选择的正实数，为了让函数曲线不是很高，可以选择
4.岭回归(Ridge regression)
经典的回归是计算
存在的问题是当数据存在复共线性（multicolinear）的时候，矩阵
几乎是奇异（singular）的，会导致最后估计结果的方差非常高。你可以通过变换矩阵的特征值来规避这个问题：
这样会减少估计的方差，但这样的估计是有偏的（biased）,如果减少的方差远远高于增加的估计偏差，MSE（Mean Square Error）减少了，那么新的方法将是有效的。
你可以将岭回归看作是一个惩罚回归，用最小化
代替最小化
，这样就避免很大的方差。
加入限制条件
问题就变成了选择恰当的k(或s)，使方差较小估计的偏差也较小。 你可以通过探索法，作图，证明来选择恰当的参数。
实际上，我们首先可以中心化和标准化X的各列（X中不应当有常数项），你可以自己编程完成：
my.lm.ridge <- function (y, x, k=.1)
{
my <- mean(y)
y <- y - my
mx <- apply(x,2,mean)
x <- x - matrix(mx, nr=dim(x)[1], nc=dim(x)[2], byrow=T)
sx <- apply(x,2,sd)
x <- x/matrix(sx, nr=dim(x)[1], nc=dim(x)[2], byrow=T)
b <- solve( t(x) %*% x + diag(k, dim(x)[2]), t(x) %*% y) #计算
c(my - sum(b*mx/sx) , b/sx)
}
或者用惩罚最小二乘的定义：
my.ridge <- function (y,x,k=0)
{
xm <- apply(x,2,mean)
ym <- mean(y)
y <- y - ym
x <- t( t(x) - xm ) #中心化
ss <- function (b)
{
t( y - x %*% b ) %*% ( y - x %*% b ) + k * t(b) %*% b
}
b <- nlm(ss, rep(0,dim(x)[2]))$estimate #非线性最小化
c(ym-t(b)%*%xm, b)
}
my.ridge.test <- function (n=20, s=.1)
{
x <- rnorm(n)
x1 <- x + s*rnorm(n)
x2 <- x + s*rnorm(n)
x <- cbind(x1,x2)
y <- x1 + x2 + 1 + rnorm(n)
lambda <- c(0, .001, .01, .1, .2, .5, 1, 2, 5, 10)
b <- matrix(nr=length(lambda), nc=1+dim(x)[2])
for (i in 1:length(lambda))
{
b[i,] <- my.ridge(y,x,lambda[i])
}
plot(b[,2], b[,3],
type="b",
xlim=range(c(b[,2],1)), ylim=range(c(b[,3],1)))
text(b[,2], b[,3], lambda, adj=c(-.2,-.2), col="blue")
points(1,1,pch="+", cex=3, lwd=3)
points(b[8,2],b[8,3],pch=15)
}
my.ridge.test()
op <- par(mfrow=c(3,3))
for (i in 1:9)
{
my.ridge.test()
}
par(op)
op <- par(mfrow=c(3,3))
for (i in 1:9)
{
my.ridge.test(20,10)
}
par(op)
在R的"MASS"包中有个函数“lm.ridge”，我们可以使用它。
library(MASS)
my.sample <- function (n=20)
{
x <- rnorm(n)
x1 <- x + .1*rnorm(n)
x2 <- x + .1*rnorm(n)
y <- 0 + x1 - x2 + rnorm(n)
cbind(y, x1, x2)
}
n <- 500
r <- matrix(NA, nr=n, nc=3)
s <- matrix(NA, nr=n, nc=3)
for (i in 1:n)
{
m <- my.sample()
r[i,] <- lm(m[,1]~m[,-1])$coef
s[i,2:3] <- lm.ridge(m[,1]~m[,-1], lambda=.1)$coef
s[i,1] <- mean(m[,1])
}
plot( density(r[,1]), xlim=c(-3,3),
main="Multicolinearity: high variance")
abline(v=0, lty=3)
lines(density(r[,2]),col='red')
lines(density(s[,2]),col='red',lty=2)
abline(v=1,col='red',lty=3)
lines(density(r[,3]),col='blue')
lines(density(s[,3]),col='blue',lty=2)
abline(v=-1,col='blue',lty=3)
生成样本数据，默认为n=20行，3列
x1,x2有很强的共线性
生成空矩阵r,s，用来保存后面的计算结果。
r 中每行是作经典线性回归的回归系数
s中每行是作岭回归的回归系数
岭回归中常数项的估计值为
绘制出经典回归和岭回归系数的核密度估计曲线,可以看到两种方法计算出来的回归系数的差别。
下面我们将演示岭回归是有偏的。
......续上面的程序......
evaluate.density <- function (d,x, eps=1e-6)
{
density(d, from=x-eps, to=x+2*eps, n=4)$y[2]
}
x<-mean(r[,2]); points( x, evaluate.density(r[,2],x) )
x<-mean(s[,2]); points( x, evaluate.density(s[,2],x) )
x<-mean(r[,3]); points( x, evaluate.density(r[,3],x) )
x<-mean(s[,3]); points( x, evaluate.density(s[,3],x) )
legend( par("usr")[1], par("usr")[4], yjust=1,
c("intercept", "x1", "x2"),
lwd=1, lty=1,
col=c(par('fg'), 'red','blue') )
legend( par("usr")[2], par("usr")[4], yjust=1, xjust=1,
c("classical regression", "ridge regression"),
lwd=1,lty=c(1,2),
col=par('fg'))
在x附近用4个点估计密度函数。
作出每个估计参数的均值点。
左侧的图例
par("usr")是作图的边界，yjust=1是指相对于y轴右对齐。
c()是说明文字,lwd=1宽度为1,lty指线型
col=当然是颜色了，par('fg')是指默认前景色
我们还可以计算估计量的均值和方差：
#LS regression
> apply(r,2,mean)
[1] 0.01031310 0.91142520 -0.92335369
> apply(r,2,var)
[1] 0.05631858 3.12950890 3.18404170
#Ridge regression
> apply(s,2,mean)
[1] 0.01116668 0.56577249 -0.57728805
> apply(s,2,var)
[1] 0.04901199 1.14735701 1.17649852
但是我们不能用它们来比较两种方法的优劣，我们更宁愿用均方误差的均值（Mean Square Error,MSE）:它是估计的参数的真实值之间差的平方的均值（它与方差的定义差不多，将真实值换为均值就变成了方差。）
> v <- matrix(c(0,1,-1), nr=n, nc=3, byrow=T)
> apply( (r-v)^2, 2, mean )
[1] 0.0563123 3.1310954 3.1835483
> apply((s-v)^2,2,mean)
[1] 0.04903866 1.33361583 1.35283092
我们也可以画出MSE随k变化而变化的过程。
