重回帰式の導出プログラム

jyuukaiki <- function(data) {
  attach(data)

  n <- length(data[[1]])
  p <- length(data)
  s <- array(1:(p-1)*(p-1), dim=c(p-1,p-1))
  syi <- array(1:p-1, dim=c(p-1,1))
  a <- array(1:p-1, dim=c(p-1,1))

  for(i in 1:p-1) {
    for(j in 1:p-1){
      Sij <- 0.0;
      for(k in 1:n){
        dif <- ( data[[i+1]][k] - mean(data[[i+1]]) ) * ( data[[j+1]][k] - mean(data[[j+1]]) )
        Sij <- Sij + dif
      }
      s[i,j] <- Sij
    }
  }
  for(j in 1:p-1){
    Syi <- 0.0
    for(k in 1:n){
      dif <- ( data[[1]][k] - mean(data[[1]]) ) * ( data[[j+1]][k] - mean(data[[j+1]]) )
      Syi <- Syi + dif
    }
    syi[j] <- Syi
  }
  Syy <- 0.0
  for(k in 1:n){
    dif <- ( data[[1]][k] - mean(data[[1]]) ) * ( data[[1]][k] - mean(data[[1]]) )
    Syy <- Syy + dif
  }
  syy <- Syy

  a <- solve(s,syi)

  x <- array(1:p-1, dim=c(1,p-1))
  for(i in 1:p-1) {
    x[i] <- a[i] * mean(data[[i+1]])
  }
  wa <- 0.0
  for(i in 1:length(x)){
    wa <- wa + x[i][1]
  }
  a0 <- mean(data[[1]]) - wa

  z <- c("")
  for(i in 1:length(a)){
    z <- paste(z, a[i][1], "X", i,"+")
  }
  z <- paste("y =", z, a0)
  detach(data)

  s
}




出力例
＜データ（教科書p.26営業所別の売上額）＞
y	x1	x2
8	5	6
9	5	8
12	7	10
11	5	12
13	8	12
17	12	12

出力結果
y =  0.862745098039216 X 1 + 0.46078431372549 X 2 + 1.01960784313725