milk.gas = milk[milk$v4 == "gasline",1:3] milk.gas milk.diesel = milk[milk$v4 == "diesel",1:3] milk.diesel n1 = dim(milk.gas)[1] n1 n2 = dim(milk.diesel)[1] n2 p = 3 y.bar.g = colMeans(milk.gas) y.bar.d = colMeans(milk.diesel) s.g = cov(milk.gas) s.d = cov(milk.diesel) s.pl = (s.g*(n1-1)+s.d*(n2-1))/(n1+n2-2) s.pl T.sq = (n1*n2)/(n1+n2)*t(y.bar.g-y.bar.d)%*%solve(s.pl)%*%(y.bar.g-y.bar.d) T.sq cut.off = (n1+n2-2)*p/(n1+n2-p-1)*qf(1-0.01,p,n1+n2-p-1) cut.off把上述R语言代码用自定义函数function表示

可以将上述R语言代码用自定义函数function表示如下：

calculate_Tsq <- function(data) {
  milk.gas <- data[data$v4 == "gasline", 1:3]
  milk.diesel <- data[data$v4 == "diesel", 1:3]
  n1 <- dim(milk.gas)[1]
  n2 <- dim(milk.diesel)[1]
  p <- 3
  y.bar.g <- colMeans(milk.gas)
  y.bar.d <- colMeans(milk.diesel)
  s.g <- cov(milk.gas)
  s.d <- cov(milk.diesel)
  s.pl <- (s.g*(n1-1)+s.d*(n2-1))/(n1+n2-2)
  T.sq <- (n1*n2)/(n1+n2)*t(y.bar.g-y.bar.d)%*%solve(s.pl)%*%(y.bar.g-y.bar.d)
  cut.off <- (n1+n2-2)*p/(n1+n2-p-1)*qf(1-0.01,p,n1+n2-p-1)
  return(list(T.sq=T.sq, cut.off=cut.off))
}

# example usage:
data <- read.csv("milk.csv")
result <- calculate_Tsq(data)
print(result$T.sq)
print(result$cut.off)