10.1 一个简单的例子

这是一个简单的例子，它用来计算双样本的 t-统计量，并且显示“所有运算步骤”。这是一个人为的例子，当然还有其他更简单的办法得到一样的结果。

函数定义如下：

> twosam <- function(y1, y2) {
    n1  <- length(y1); n2  <- length(y2)
    yb1 <- mean(y1);   yb2 <- mean(y2)
    s1  <- var(y1);    s2  <- var(y2)
    s <- ((n1-1)*s1 + (n2-1)*s2)/(n1+n2-2)
    tst <- (yb1 - yb2)/sqrt(s*(1/n1 + 1/n2))
    tst
  }

通过这个函数，你可以利用下面的命令实现双样本 t-检验：

> tstat <- twosam(data$male, data$female); tstat

第二个例子是仿效 Matlab 里面的反斜杠命令。它是用来返回向量 y 正交投影到 X 列空间上面的系数。(这常常被称为回归系数的最小二乘法估计。) 这可以用函数 qr() 来实现；但是直接使用这个函数有时需要一点技巧。下面提供了一个简单而又安全的 函数。

给定一个 n×1 的向量 $y$ 和一个 n×p 的矩阵 $X$，因此 $X$ \ $y$ 可以定义如下$(X′X)^{-1}X′y$, 其中 $(X′X)^{-1}$ 这就是 $X′X$ 的广义逆矩阵(generalized inverse)。

> bslash <- function(X, y) {
  X <- qr(X)
  qr.coef(X, y)
}

当这个对象创建后，它可用于这样的命令中：

> regcoeff <- bslash(Xmat, yvar)

经典的 R 函数 lsfit() 可以很快地完成这个功能和做其他一些相关的事情¹。它 以一种有点有悖直觉的方式依次用函数 qr() 和 qr.coef() 去完成这部分计算。如果一部分代码常常被使用，我们可以把这部分代码单独列出来成为函数使用。如果是这样考虑，我们可能期望矩阵的二元操作(binary operator)能以一种更为便利的方式进行。

¹. 参见R中的统计模型中描述的方法 ↩