R语言基础安装鸢尾花数据集iris多组均值检验

鸢尾花数据集iris包含四个属性：'尊片长度' (Sepal.Length)，'尊片宽度' (Sepal.width)，'花瓣长度' (Petal.Length)，'花瓣宽度' (Petal.width)，数据集含有150个样本，分为三类：'Setosa'、'Versicolor'、'Virginica'。

本案例将对这三类鸢尾花进行多组均值检验，检验其四个属性的均值是否相等。

检验步骤

第一类鸢尾花Setosa的四个属性的均值记为p1，同时记Mo=(5.0,3.4,1.5,0.2)'，检验H0:Pi=Po-H1:mfH。
- 计算样本均值和样本协方差矩阵：

data1 <- iris[1:50, 1:4] # 取第一类鸢尾花Setosa的四个属性 p1 <- colMeans(data1) # 计算样本均值 S1 <- cov(data1) # 计算样本协方差矩阵 Mo <- c(5.0, 3.4, 1.5, 0.2) # 均值假设值


    * 计算检验统计量：

    ```R
n1 <- nrow(data1)  # 样本容量
Tsq <- (n1 - 1) * t(p1 - Mo) %*% solve(S1) %*% (p1 - Mo)  # Hotelling's T-squared统计量

* 根据F分布计算拒绝域：

```R

p <- length(p1) # 变量个数 df1 <- p df2 <- n1 - p alpha <- 0.05 # 检验水平 F_crit <- qf(1 - alpha, df1, df2) # F分布的临界值


    * 判断是否拒绝原假设：

    ```R
if(Tsq > F_crit) {
  cat("拒绝原假设，样本均值与假设值不同。\n")
} else {
  cat("不拒绝原假设，样本均值与假设值相同。\n")
}

结果为：拒绝原假设，样本均值与假设值不同。

把第二类鸢尾花Versicolor的四个属性的均值记为p2，检验H0:P2=Po-H1:ufp，假设第一类和第二类鸢尾花的四个属性的协差阵不相等。
- 计算样本均值和样本协方差矩阵：

data2 <- iris[51:100, 1:4] # 取第二类鸢尾花Versicolor的四个属性 p2 <- colMeans(data2) # 计算样本均值 S2 <- cov(data2) # 计算样本协方差矩阵


    * 计算检验统计量：

    ```R
n2 <- nrow(data2)  # 样本容量
Sp <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)  # 求出混合协方差矩阵
Tsq <- (n1 * n2 / (n1 + n2)) * t(p1 - p2) %*% solve(Sp) %*% (p1 - p2)  # Hotelling's T-squared统计量

* 根据F分布计算拒绝域：

```R

df1 <- p df2 <- n1 + n2 - p - 1 alpha <- 0.05 # 检验水平 F_crit <- qf(1 - alpha, df1, df2) # F分布的临界值


    * 判断是否拒绝原假设：

    ```R
if(Tsq > F_crit) {
  cat("拒绝原假设，样本均值与假设值不同。\n")
} else {
  cat("不拒绝原假设，样本均值与假设值相同。\n")
}

结果为：拒绝原假设，样本均值与假设值不同。

把第三类鸢尾花virginica的四个属性的均值记为u，检验H0:P1=P2=P3-H1:P不全相等。
- 计算样本均值和样本协方差矩阵：

data3 <- iris[101:150, 1:4] # 取第三类鸢尾花virginica的四个属性 p3 <- colMeans(data3) # 计算样本均值 S3 <- cov(data3) # 计算样本协方差矩阵


    * 计算检验统计量：

    ```R
Sp <- ((n1 - 1) * S1 + (n2 - 1) * S2 + (n3 - 1) * S3) / (n1 + n2 + n3 - p - 1)  # 求出混合协方差矩阵
Tsq1 <- (n1 * n2 / (n1 + n2)) * t(p1 - p2) %*% solve(Sp) %*% (p1 - p2)  # 计算P1和P2的Hotelling's T-squared统计量
Tsq2 <- (n2 * n3 / (n2 + n3)) * t(p2 - p3) %*% solve(Sp) %*% (p2 - p3)  # 计算P2和P3的Hotelling's T-squared统计量
Tsq3 <- (n1 * n3 / (n1 + n3)) * t(p1 - p3) %*% solve(Sp) %*% (p1 - p3)  # 计算P1和P3的Hotelling's T-squared统计量
Tsq <- (Tsq1 + Tsq2 + Tsq3) / 3  # 计算P1、P2和P3的平均Hotelling's T-squared统计量

* 根据F分布计算拒绝域：

```R

df1 <- p * (n1 + n2 + n3 - p - 1) df2 <- n1 + n2 + n3 - p - 1 alpha <- 0.05 # 检验水平 F_crit <- qf(1 - alpha, df1, df2) # F分布的临界值


    * 判断是否拒绝原假设：

    ```R
if(Tsq > F_crit) {
  cat("拒绝原假设，样本均值不全相等。\n")
} else {
  cat("不拒绝原假设，样本均值全相等。\n")
}

结果为：拒绝原假设，样本均值不全相等。

结论

经过检验，我们发现三类鸢尾花四个属性的均值不全相等。

代码附后

# 载入iris数据集
data(iris)

# 1. 检验第一类鸢尾花
data1 <- iris[1:50, 1:4]
p1 <- colMeans(data1)
S1 <- cov(data1)
Mo <- c(5.0, 3.4, 1.5, 0.2)
n1 <- nrow(data1)
p <- length(p1)
df1 <- p
df2 <- n1 - p
alpha <- 0.05
Tsq <- (n1 - 1) * t(p1 - Mo) %*% solve(S1) %*% (p1 - Mo)
F_crit <- qf(1 - alpha, df1, df2)
if(Tsq > F_crit) {
  cat("拒绝原假设，样本均值与假设值不同。\n")
} else {
  cat("不拒绝原假设，样本均值与假设值相同。\n")
}

# 2. 检验第二类鸢尾花
data2 <- iris[51:100, 1:4]
p2 <- colMeans(data2)
S2 <- cov(data2)
n2 <- nrow(data2)
Sp <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)
Tsq <- (n1 * n2 / (n1 + n2)) * t(p1 - p2) %*% solve(Sp) %*% (p1 - p2)
df1 <- p
df2 <- n1 + n2 - p - 1
F_crit <- qf(1 - alpha, df1, df2)
if(Tsq > F_crit) {
  cat("拒绝原假设，样本均值与假设值不同。\n")
} else {
  cat("不拒绝原假设，样本均值与假设值相同。\n")
}

# 3. 检验第三类鸢尾花
data3 <- iris[101:150, 1:4]
p3 <- colMeans(data3)
S3 <- cov(data3)
n3 <- nrow(data3)
Sp <- ((n1 - 1) * S1 + (n2 - 1) * S2 + (n3 - 1) * S3) / (n1 + n2 + n3 - p - 1)
Tsq1 <- (n1 * n2 / (n1 + n2)) * t(p1 - p2) %*% solve(Sp) %*% (p1 - p2)
Tsq2 <- (n2 * n3 / (n2 + n3)) * t(p2 - p3) %*% solve(Sp) %*% (p2 - p3)
Tsq3 <- (n1 * n3 / (n1 + n3)) * t(p1 - p3) %*% solve(Sp) %*% (p1 - p3)
Tsq <- (Tsq1 + Tsq2 + Tsq3) / 3
df1 <- p * (n1 + n2 + n3 - p - 1)
df2 <- n1 + n2 + n3 - p - 1
F_crit <- qf(1 - alpha, df1, df2)
if(Tsq > F_crit) {
  cat("拒绝原假设，样本均值不全相等。\n")
} else {
  cat("不拒绝原假设，样本均值全相等。\n")
}