2  描述统计

2.1 引言

本章将介绍 描述统计 的相关内容，在介绍 描述统计 的过程中，我们将学习如何描述和总结数据。

小节 2.2 会讨论如何描述数据集。如果一个数据集中的数据差异相对较小，小节 2.2.1 和 小节 2.2.2 介绍了如何使用 频率表（frequency tables）来描述这样的数据集，小节 2.2.3 则介绍如何将数据集划分到不同的分组。

小节 2.3 会讨论如何使用统计量来总结数据集，统计量往往是由数据集确定的量化数值。小节 2.3.1 介绍了用于描述数据集中心（center）的三个统计量：样本均值（sample mean）、样本中位数（sample median）和样本众数（sample mode）。小节 2.3.2 会介绍样本方差（sample variance）及样本标准差（sample standard deviation）。样本均值、样本中位数、样本众数、样本方差和标准差这些统计数据用于描述数据集中数值的分布。小节 2.3.3 会介绍样本的百分位数（sample percentiles），样本百分位数这个统计量可以告诉我们数据集中哪个值大于所有数据的 95%。

小节 2.4 节会介绍数据样本集中的切比雪夫不等式（Chebyshev’s inequality），切比雪夫不等式给出了与样本均值相差超过 \(k\) 倍样本标准差的数据比例的上限。虽然切比雪夫不等式适用于所有数据集，但在某些情况下（例如 小节 2.5 中讨论的正态分布数据集），我们可以获得更精准的数据比例。

在 小节 2.5 中，我们注意到，当数据分布符合钟形（bell-shaped）时，我们认为该数据集是近似正态的，此时，位于样本均值 \(k\) 个样本标准差内的数据比例可以根据所谓的经验法则（empirical rule）给出更精确的估计。

小节 2.6 中关注的是成对数据（paired data），并介绍了如何使用散点图（scatter diagram）来描述成对数据以及样本相关系数（sample correlation coefficient）。样本相关系数是一个统计量，该值表示成对数据中第一个值与第二个值之间相关联的程度。

小节 2.7 会介绍国民收入分配，洛伦兹曲线（Lorenz curve）（ L(p)）和基尼系数（Gini index）。洛伦兹曲线给出了收入较低的百分之 \(p\) 的居民收入占总收入的比例，基尼系数则是用于衡量收入分配不平等程度的指标。

小节 2.8 会介绍如何使用 R 来分析数据集。

2.2 描述数据集

应当用数据清晰、简洁的显示研究结果，以便研究者可以迅速把握数据的关键特征。多年来，人们发现，对于展示数据而言，表格和图是特别有用的方式。图表常常可以揭示数据的重要特征，如数据的范围、集中程度和对称性。本节将介绍一些常见的数据图表展示方式。

2.2.1 频率表和图

如果一个数据集中的数据差异较小，我们便可以以频率表（frequency table）的形式来方便的展现该数据集。例如，表 2.1 就是一个数据集的频率表，该数据集由 42 名最近毕业的电子工程专业的大学生的起始年薪（以千美元为单位，四舍五入）构成。表 2.1 告诉我们，最低起薪为 57000 美元，并且有四名毕业生的起薪为 57000 美元；最高起薪为 70000 美元，并且有一名学生的起薪为 70000 美元；最为普遍的起薪为 62000 美元，并且有 10 名学生拿到了 62000 美元的起薪。

代码 
library(knitr)
ss <- c(57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 70)
f <- c(4, 1, 3, 5, 8, 10, 0, 5, 2, 3, 1)
df <- data.frame("Starting Salary" = ss, "Frequency" = f)
kable(df, align = "l")
表 2.1: 起薪分布
Starting.Salary Frequency
57 4
58 1
59 3
60 5
61 8
62 10
63 0
64 5
66 2
67 3
70 1

可以用线图（line graph）来直观地显示频率表中的数据。如 图 2.1，在线图中，水平坐标轴用于绘制不同的数据值，而垂直直线的高度用于表示对应数据值出现的频率。

代码 
library(ggplot2)  

df <- data.frame(  
  salary = c(57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 70),  
  value_end = c(4, 1, 3, 5, 8, 10, 0, 5, 2, 3, 1),  
  value_start = rep(0, 11) 
)  

# 将类别转换为因子，以确保正确的顺序  
df$salary <- factor(df$salary, levels = df$salary, ordered = TRUE)  
 
ggplot(df, aes(x = salary, y = value_start, yend = value_end, color = salary)) +  
  geom_segment(size = 2) +                    # 绘制水平线表示每个“条形”  
  geom_point(aes(y = value_end), size = 3) +  # 在每个“条形”的末端添加点  
  theme_minimal() +  
  labs(x = "Starting Salary", y = "Frequency") +  
  scale_y_continuous(breaks = seq(1, 10, by = 1)) +  
  theme(legend.position = "none")           # 调整图例位置
图 2.1: 起薪线图

当为 图 2.1 中的线条增加厚度时，图 2.1 就变成了 图 2.2 所示的条形图（bar graph）。

代码 
library(ggplot2)  

df <- data.frame(  
  salary = c(57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 70),  
  value = c(4, 1, 3, 5, 8, 10, 0, 5, 2, 3, 1)
)  

# 将类别转换为因子，以确保正确的顺序  
df$salary <- factor(df$salary, levels = df$salary, ordered = TRUE)  
 
ggplot(df, aes(x = salary, y = value)) +  
  geom_bar(stat = "identity", fill = "steelblue") +  
  theme_minimal() +  
  xlab("Starting Salary") +  
  ylab("Frequency") +
  scale_y_continuous(breaks = seq(1, 10, by = 2))
图 2.2: 起薪条形图

另一种用于表示频数表的图形是折线图（frequency polygon），在折线图中，对于在数据集中出现的不同的数据值，其出现的次数为纵坐标，以此画出对应的点，然后用直线将这些点依次连接起来。图 2.3 给出了 表 2.1 对应的折线图。

代码 
library(ggplot2)  

df <- data.frame(  
  salary = c(57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 70),  
  value = c(4, 1, 3, 5, 8, 10, 0, 5, 2, 3, 1)
)  

ggplot(df, aes(x=salary, y=value)) +   
  geom_line() +  
  labs(x="Starting Salary", y="Frequency") + 
  scale_y_continuous(breaks = seq(1, 10, by = 2)) + 
  theme_minimal()
图 2.3: 起薪折线图

2.2.2 相对频率表和图

一个数据集中有 \(n\) 个值，如果 \(f\) 是某个特定值的频率，那么\(\frac{f}{n}\) 为该特定值的相对频率（relative frequency）。也就是说，一个值的相对频率是数据集中具有该值的数据所占的比例。可以用线图、条形图、折线图来表示相对频率。事实上，在相对频率图中，除了纵坐标的值是绝对频率图中的值除以数据点总数外，相对频率图看起来和绝对频率图没什么不同。

例 2.1 表 2.2 是 表 2.1 的相对频率表。将 表 2.1 中相应的频率除以数据集的大小（42） 得到其对应的相对频率表。

代码 
library(knitr)
ss <- c(57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 70)
f <- c(4, 1, 3, 5, 8, 10, 0, 5, 2, 3, 1)
cnt = sum(f)
f <- f / cnt

df <- data.frame("Starting Salary" = ss, "Relative Frequency" = f)
kable(df, align = "l")
表 2.2: 起薪分布相对频率表
Starting.Salary Relative.Frequency
57 0.0952381
58 0.0238095
59 0.0714286
60 0.1190476
61 0.1904762
62 0.2380952
63 0.0000000
64 0.1190476
66 0.0476190
67 0.0714286
70 0.0238095

当数据集中的数据不是数值型数据时，通常使用饼图（pie chart）来表示相对频率。首先画一个圆，然后将其切分成不同的扇形区域，每一个扇区对应数据集中的一种数据类型。在饼图中，扇区的面积表示数据值的相对频率，其面积等于圆的总面积乘以数据值的相对频率。

例 2.2 根据一家肿瘤专科诊所登记的最近 200 名患者的数据，得到了一个关于不同类型肿瘤患者数的数据集。图 2.4 使用饼图的形式来显示如上的数据集。

代码 
library(knitr)

type <- c("Melanoma", "Bladder", "Colon", "Lung", "Breast", "Prostate")
num <- c(9, 12, 32, 42, 50, 55)  
value <- num / sum(num)

df <- data.frame("Type of Cancer" = type, "Number of New Cases" = num, "Relative Frequency" = value)

kable(df, align = "l")
Type.of.Cancer Number.of.New.Cases Relative.Frequency
Melanoma 9 0.045
Bladder 12 0.060
Colon 32 0.160
Lung 42 0.210
Breast 50 0.250
Prostate 55 0.275
代码 
library(ggplot2)
  
# 提供的数据  
type <- c("Melanoma", "Bladder", "Colon", "Lung", "Breast", "Prostate")
num <- c(9, 12, 32, 42, 50, 55)  
value <- num / sum(num)

df <- data.frame(type = type, num = num, value = value)
df$label <- paste(df$type, "\n", round(df$value * 100, 1), "%", sep="") 

# 绘制饼图  
ggplot(df, aes(x = "", y = value, fill = type)) +  
  geom_bar(width = 1, stat = "identity") +  
  coord_polar("y", start = 0) +  
  theme_void() +  
  theme(legend.title = element_blank()) +  
  scale_fill_brewer(palette = "Pastel1") +  
  labs(fill = "type") +  
  geom_text(aes(label = label),   
            position = position_stack(vjust = 0.5)) + # 添加百分比标签
  theme(legend.position = "none")
图 2.4: 诊所数据饼图
注记

虽然我们经常使用饼图，尤其是在商业分析中，饼图更是得到了广泛的使用，但是统计学家则一直很排斥使用饼图。因为人们对长度的判断比面积更为精确，因此统计学家更倾向于使用条形图。在 R 中，对饼图的支持要比其他的工具差很多。

2.2.3 分组数据，直方图，肩形图和茎叶图

如 小节 2.2.2 所述，使用线图或条形图绘制数据值的频率通常是描述数据集的有效方法。然而，对于那些不同数值数量太多的数据集而言，就无法使用线图或条形图的方法来对数据进行描述。在这种情况下，我们可以先将数值分组或按类别进行分类，然后绘制落在每个分组中的数据值的数量。分组数量的选取需要在以下两者之间进行取舍：

  • 较少的分组，但是将会丢失实际数据值中的很多信息；
  • 较多多的分组，但是将导致每个分组中的数值频率太小，以至于无法从图中辨别出数据模型。

尽管典型的分组数为 5 ~ 10 个，但是适当的分组数量来自主观选择。当然，我们可以尝试不同的分组数，并对比看看哪张图表的效果最好。尽管不是必须的，但是在选择分组数时，通常会让不同分组之间的间隔长度保持一致。

分组的端点（endpoint）称之为分组边界（class boundaries）。我们将会采用左端包含原则，该原则规定分组中的数据包含其左端点，但不包含其右端。因此，例如，20–30 这个分组包含所有 \(\ge\) 20 且 \(<\) 30 的值。

表 2.3 给出了 200 只白炽灯的使用寿命。表 2.4 为 表 2.3 数据的分组频率表，其分组的间隔长度为 100，第一个分组从 500 开始。

代码 
library(knitr)

data <- c(1067, 919, 855, 1092, 1157, 1195, 1022, 978, 923, 1333, 521, 933, 930, 807, 999, 932, 901, 1324, 996, 780, 1187, 1067, 824, 653, 844, 814, 1037, 1151, 1026, 1147, 1039, 1083, 1023, 984, 1134, 932, 998, 996, 610, 916, 1196, 785, 1162, 1170, 1195, 1340, 832, 1009, 811, 1217, 928, 1153, 954, 1063, 1035, 944, 818, 1250, 900, 1106, 1118, 1037, 980, 935, 1103, 1000, 863, 990, 883, 867, 1040, 1289, 856, 924, 938, 1078, 1133, 765, 1001, 895, 1126, 936, 929, 950, 1122, 938, 1157, 1151, 1085, 896, 946, 858, 1002, 909, 1049, 940, 1203, 1078, 704, 621, 958, 760, 878, 934, 788, 1143, 1035, 1112, 990, 1258, 699, 1083, 801, 1122, 1180, 1106, 775, 1105, 709, 860, 918, 1156, 920, 948, 905, 972, 1035, 1045, 970, 1237, 956, 1102, 1009, 765, 958, 902, 958, 1311, 1037, 702, 1071, 1069, 830, 1063, 1077, 1021, 1062, 1157, 1122, 1115, 833, 1320, 890, 1303, 1011, 1102, 854, 1178, 1138, 951, 1101, 949, 992, 966, 910, 1058, 730, 980, 935, 1069, 1170, 1067, 931, 970, 932, 904, 1192, 922, 1150, 1091, 880, 1029, 658, 912, 1292, 1116, 880, 1173, 1184, 954, 824, 529, 1081, 1171, 705, 1425, 1110, 1149, 972, 1002)

kable(matrix(data, nrow = 20, ncol = 10, byrow = TRUE))
表 2.3: 200 个 白炽灯的寿命（Hour）
1067 919 855 1092 1157 1195 1022 978 923 1333
521 933 930 807 999 932 901 1324 996 780
1187 1067 824 653 844 814 1037 1151 1026 1147
1039 1083 1023 984 1134 932 998 996 610 916
1196 785 1162 1170 1195 1340 832 1009 811 1217
928 1153 954 1063 1035 944 818 1250 900 1106
1118 1037 980 935 1103 1000 863 990 883 867
1040 1289 856 924 938 1078 1133 765 1001 895
1126 936 929 950 1122 938 1157 1151 1085 896
946 858 1002 909 1049 940 1203 1078 704 621
958 760 878 934 788 1143 1035 1112 990 1258
699 1083 801 1122 1180 1106 775 1105 709 860
918 1156 920 948 905 972 1035 1045 970 1237
956 1102 1009 765 958 902 958 1311 1037 702
1071 1069 830 1063 1077 1021 1062 1157 1122 1115
833 1320 890 1303 1011 1102 854 1178 1138 951
1101 949 992 966 910 1058 730 980 935 1069
1170 1067 931 970 932 904 1192 922 1150 1091
880 1029 658 912 1292 1116 880 1173 1184 954
824 529 1081 1171 705 1425 1110 1149 972 1002
代码 
library(knitr)

data <- c(1067, 919, 855, 1092, 1157, 1195, 1022, 978, 923, 1333, 521, 933, 930, 807, 999, 932, 901, 1324, 996, 780, 1187, 1067, 824, 653, 844, 814, 1037, 1151, 1026, 1147, 1039, 1083, 1023, 984, 1134, 932, 998, 996, 610, 916, 1196, 785, 1162, 1170, 1195, 1340, 832, 1009, 811, 1217, 928, 1153, 954, 1063, 1035, 944, 818, 1250, 900, 1106, 1118, 1037, 980, 935, 1103, 1000, 863, 990, 883, 867, 1040, 1289, 856, 924, 938, 1078, 1133, 765, 1001, 895, 1126, 936, 929, 950, 1122, 938, 1157, 1151, 1085, 896, 946, 858, 1002, 909, 1049, 940, 1203, 1078, 704, 621, 958, 760, 878, 934, 788, 1143, 1035, 1112, 990, 1258, 699, 1083, 801, 1122, 1180, 1106, 775, 1105, 709, 860, 918, 1156, 920, 948, 905, 972, 1035, 1045, 970, 1237, 956, 1102, 1009, 765, 958, 902, 958, 1311, 1037, 702, 1071, 1069, 830, 1063, 1077, 1021, 1062, 1157, 1122, 1115, 833, 1320, 890, 1303, 1011, 1102, 854, 1178, 1138, 951, 1101, 949, 992, 966, 910, 1058, 730, 980, 935, 1069, 1170, 1067, 931, 970, 932, 904, 1192, 922, 1150, 1091, 880, 1029, 658, 912, 1292, 1116, 880, 1173, 1184, 954, 824, 529, 1081, 1171, 705, 1425, 1110, 1149, 972, 1002)

intervals <- seq(500, 1500, by = 100)
group <- cut(data, intervals, right=FALSE)
kable(table("Class Interval" = group))
表 2.4: 200 个 白炽灯的寿命（Hour）分布表
Class.Interval Freq
[500,600) 2
[600,700) 5
[700,800) 12
[800,900) 25
[900,1e+03) 58
[1e+03,1.1e+03) 41
[1.1e+03,1.2e+03) 43
[1.2e+03,1.3e+03) 7
[1.3e+03,1.4e+03) 6
[1.4e+03,1.5e+03) 1

分组数据的条形图中，其不同的分组彼此相邻，我们称之为直方图（histogram）。直方图的纵坐标可以表示分组的频率或相对频率，当表示分组频率时该图为频率直方图，当表示分组相对频率时为相对频率直方图。图 2.5 给出了 表 2.4 中数据的频率直方图。

代码 
library(ggplot2)

data <- c(1067, 919, 855, 1092, 1157, 1195, 1022, 978, 923, 1333, 521, 933, 930, 807, 999, 932, 901, 1324, 996, 780, 1187, 1067, 824, 653, 844, 814, 1037, 1151, 1026, 1147, 1039, 1083, 1023, 984, 1134, 932, 998, 996, 610, 916, 1196, 785, 1162, 1170, 1195, 1340, 832, 1009, 811, 1217, 928, 1153, 954, 1063, 1035, 944, 818, 1250, 900, 1106, 1118, 1037, 980, 935, 1103, 1000, 863, 990, 883, 867, 1040, 1289, 856, 924, 938, 1078, 1133, 765, 1001, 895, 1126, 936, 929, 950, 1122, 938, 1157, 1151, 1085, 896, 946, 858, 1002, 909, 1049, 940, 1203, 1078, 704, 621, 958, 760, 878, 934, 788, 1143, 1035, 1112, 990, 1258, 699, 1083, 801, 1122, 1180, 1106, 775, 1105, 709, 860, 918, 1156, 920, 948, 905, 972, 1035, 1045, 970, 1237, 956, 1102, 1009, 765, 958, 902, 958, 1311, 1037, 702, 1071, 1069, 830, 1063, 1077, 1021, 1062, 1157, 1122, 1115, 833, 1320, 890, 1303, 1011, 1102, 854, 1178, 1138, 951, 1101, 949, 992, 966, 910, 1058, 730, 980, 935, 1069, 1170, 1067, 931, 970, 932, 904, 1192, 922, 1150, 1091, 880, 1029, 658, 912, 1292, 1116, 880, 1173, 1184, 954, 824, 529, 1081, 1171, 705, 1425, 1110, 1149, 972, 1002)

df <- data.frame(value = data)
breaks_vector <- seq(from = 500, to = 1500, by = 100)

ggplot(df, aes(x = value)) +  
  geom_histogram(breaks = breaks_vector, fill = "lightblue", color = "black") +  
  xlab("incandescent lamps lifetime(Hour)") +  
  ylab("Frequency")
图 2.5: 200 个 白炽灯的寿命（Hour）分布直方图

很多时候，我们需要绘制累积频率图（cumulative frequency graph）或者累积相对频率图（cumulative relative frequency graph）。在累积频率（相对频率）图的横坐标上，一个点代表一个可能的数据值，其对应的纵坐标为值小于或等于该数据值的数据数量（或比例）。图 2.6 给出了 表 2.3 中所示数据的累积相对频率图。从这个图中我们可以发现，使用寿命小于 1500（小时） 的白炽灯占比为 100%，使用寿命 \(\le\) 900 的占比大约为 40% ，使用寿命 \(\le\) 1100 的占比大约为 80%，……我们通常又把累积频率图称为肩形图（ogive）。

代码 
library(ggplot2)

data <- c(1067, 919, 855, 1092, 1157, 1195, 1022, 978, 923, 1333, 521, 933, 930, 807, 999, 932, 901, 1324, 996, 780, 1187, 1067, 824, 653, 844, 814, 1037, 1151, 1026, 1147, 1039, 1083, 1023, 984, 1134, 932, 998, 996, 610, 916, 1196, 785, 1162, 1170, 1195, 1340, 832, 1009, 811, 1217, 928, 1153, 954, 1063, 1035, 944, 818, 1250, 900, 1106, 1118, 1037, 980, 935, 1103, 1000, 863, 990, 883, 867, 1040, 1289, 856, 924, 938, 1078, 1133, 765, 1001, 895, 1126, 936, 929, 950, 1122, 938, 1157, 1151, 1085, 896, 946, 858, 1002, 909, 1049, 940, 1203, 1078, 704, 621, 958, 760, 878, 934, 788, 1143, 1035, 1112, 990, 1258, 699, 1083, 801, 1122, 1180, 1106, 775, 1105, 709, 860, 918, 1156, 920, 948, 905, 972, 1035, 1045, 970, 1237, 956, 1102, 1009, 765, 958, 902, 958, 1311, 1037, 702, 1071, 1069, 830, 1063, 1077, 1021, 1062, 1157, 1122, 1115, 833, 1320, 890, 1303, 1011, 1102, 854, 1178, 1138, 951, 1101, 949, 992, 966, 910, 1058, 730, 980, 935, 1069, 1170, 1067, 931, 970, 932, 904, 1192, 922, 1150, 1091, 880, 1029, 658, 912, 1292, 1116, 880, 1173, 1184, 954, 824, 529, 1081, 1171, 705, 1425, 1110, 1149, 972, 1002)

df <- data.frame(value = data)
df$value <- df$value[order(df$value)]
df$cumulative_freq <- cumsum(rep(1/nrow(df), nrow(df))) 

ggplot(df, aes(x = value, y = cumulative_freq)) +  
  geom_step(direction = "vh", color = "blue") +  
  xlab("incandescent lamps lifetime") +  
  ylab("Cumulative Relative Frequency") +   
  theme_minimal()
图 2.6: 200 个 白炽灯的寿命（Hour）累积频率图

对于中小型数据集而言，茎叶图（stem-leaf plot）是一种有效的数据组织方式。可以把每个数据值分为茎和叶两部分来获得茎叶图。例如，如果数据都是两位整数值，那么我们可以把十位数作为茎，个位数作为叶。例如，值 62 的茎叶图表示为：

62 的茎叶图
Stem Leaf
  6   2

62，67 两个数据的茎叶图表示为：

62，67 的茎叶图
Stem Leaf
  6  2,7

例 2.3 表 2.5 给出了美国 35 个城市的每日最低气温的月平均值和年平均值。表 2.5 中的每日最低温度的年平均值可以用下面的茎叶图表示。

代码 
df <- read.table("../misc/US_minimum_temperature.csv", header=TRUE, sep=",")
stem(df$Annual.avg.)

  The decimal point is 1 digit(s) to the right of the |

  2 | 9
  3 | 034
  3 | 56699
  4 | 0001244
  4 | 55567899
  5 | 112
  5 | 677899
  6 | 
  6 | 9
  7 | 0
注记

由于 R 中的 stem() 并不支持小数，因此对于 表 2.5 中所示的年平均值会先进行四舍五入转成整数，然后再画茎叶图。因此，这里给的茎叶图和原书中的不一致。

代码 
library(knitr)
df <- read.table("../misc/US_minimum_temperature.csv", header=TRUE, sep=",")
kable(df)
表 2.5: 1961年~1990年，美国指定城市的每日最低气温表。数据来源：Source: U.S. National Oceanic and Atmospheric Administration, Climatography of the United States, No. 81。
State Station Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec. Annual.avg.
AL Mobile 40.0 42.7 50.1 57.1 64.4 70.7 73.2 72.9 68.7 57.3 49.1 43.1 57.4
AK Juneau 19.0 22.7 26.7 32.1 38.9 45.0 48.1 47.3 42.9 37.2 27.2 22.6 34.1
AZ Phoenix 41.2 44.7 48.8 55.3 63.9 72.9 81.0 79.2 72.8 60.8 48.9 41.8 59.3
AR Little Rock 29.1 33.2 42.2 50.7 59.0 67.4 71.5 69.8 63.5 50.9 41.5 33.1 51.0
CA Los Angeles 47.8 49.3 50.5 52.8 56.3 59.5 62.8 64.2 63.2 59.2 52.8 47.9 55.5
CA Sacramento 37.7 41.4 43.2 45.5 50.3 55.3 58.1 58.0 55.7 50.4 43.4 37.8 48.1
CA San Diego 48.9 50.7 52.8 55.6 59.1 61.9 65.7 67.3 65.6 60.9 53.9 48.8 57.6
CA San Francisco 41.8 45.0 45.8 47.2 49.7 52.6 53.9 55.0 55.2 51.8 47.1 42.7 49.0
CO Denver 16.1 20.2 25.8 34.5 43.6 52.4 58.6 56.9 47.6 36.4 25.4 17.4 36.2
CT Hartford 15.8 18.6 28.1 37.5 47.6 56.9 62.2 60.4 51.8 40.7 32.8 21.3 39.5
DE Wilmington 22.4 24.8 33.1 41.8 52.2 61.6 67.1 65.9 58.2 45.7 37.0 27.6 44.8
DC Washington 26.8 29.1 37.7 46.4 56.6 66.5 71.4 70.0 62.5 50.3 41.1 31.7 49.2
FL Jacksonville 40.5 43.3 49.2 54.9 62.1 69.1 71.9 71.8 69.0 59.3 50.2 43.4 57.1
FL Miami 59.2 60.4 64.2 67.8 72.1 75.1 76.2 76.7 75.9 72.1 66.7 61.5 69.0
GA Atlanta 31.5 34.5 42.5 50.2 58.7 66.2 69.5 69.0 63.5 51.9 42.8 35.0 51.3
HI Honolulu 65.6 65.4 67.2 68.7 70.3 72.2 73.5 74.2 73.5 72.3 70.3 67.0 70.0
ID Boise 21.6 27.5 31.9 36.7 43.9 52.1 57.7 56.8 48.2 39.0 31.1 22.5 39.1
IL Chicago 12.9 17.2 28.5 38.6 47.7 57.5 62.6 61.6 53.9 42.2 31.6 19.1 39.5
IL Peoria 13.2 17.7 29.8 40.8 50.9 60.7 65.4 63.1 55.2 43.1 32.5 19.3 41.0
IN Indianapolis 17.2 20.9 31.9 41.5 51.7 61.0 65.2 62.8 55.6 43.5 34.1 23.2 42.4
IA Des Moines 10.7 15.6 27.6 40.0 51.5 61.2 66.5 63.6 54.5 42.7 29.9 16.1 40.0
KS Wichita 19.2 23.7 33.6 44.5 54.3 64.6 69.9 67.9 59.2 46.6 33.9 23.0 45.0
KY Louisville 23.2 26.5 36.2 45.4 54.7 62.9 67.3 65.8 58.7 45.8 37.3 28.6 46.0
LA New Orleans 41.8 44.4 51.6 58.4 65.2 70.8 73.1 72.8 69.5 58.7 51.0 44.8 58.5
ME Portland 11.4 13.5 24.5 34.1 43.4 52.1 58.3 57.1 48.9 38.3 30.4 17.8 35.8
MD Baltimore 23.4 25.9 34.1 42.5 52.6 61.8 66.8 65.7 58.4 45.9 37.1 28.2 45.2
MA Boston 21.6 23.0 31.3 40.2 49.8 59.1 65.1 64.0 56.8 46.9 38.3 26.7 43.6
MI Detroit 15.6 17.6 27.0 36.8 47.1 56.3 61.3 59.6 52.5 40.9 32.2 21.4 39.0
MI Sault Ste. Marie 4.6 4.8 15.3 28.4 38.4 45.5 51.3 51.3 44.3 36.2 25.9 11.8 29.8
MN Duluth -2.2 2.8 15.7 28.9 39.6 48.5 55.1 53.3 44.5 35.1 21.5 4.9 29.0
MN Minneapolis-St. Paul 2.8 9.2 22.7 36.2 47.6 57.6 63.1 60.3 50.3 38.8 25.2 10.2 35.3
MS Jackson 32.7 35.7 44.1 51.9 60.0 67.1 70.5 69.7 63.7 50.3 42.3 36.1 52.0
MO Kansas City 16.7 21.8 32.6 43.8 53.9 63.1 68.2 65.7 56.9 45.7 33.6 21.9 43.7
MO St. Louis 20.8 25.1 35.5 46.4 56.0 65.7 70.4 67.9 60.5 48.3 37.7 26.0 46.7
MT Great Falls 11.6 17.2 22.8 31.9 40.9 48.6 53.2 52.2 43.5 35.8 24.3 14.6 33.1

2.3 对数据集进行汇总

如今的实验涉及到的数据规模通常比较大。例如，为了了解某些行为习惯对健康的影响，1951 年，医学统计学家 R. Doll 和 A. B. Hill 向英国的所有医生发起了一项调查问卷，并且回收到了大约 40,000 份问卷。问卷的问题涉及年龄、饮食习惯和吸烟习惯。随后，Doll 和 Hill 对问卷回复者进行了为期 10 年的跟踪调查，并对其中的死亡人员的死因进行了监测。为了感知如此大量的数据，需要选择适当的方法对数据进行总结。在本节中，我们将介绍一些汇总数据的统计量，其中统计量是一个数值，数据集决定了其统计量的具体数值。

The British Doctors’ Study (1951–2001)

1951 年，R. Doll 和 A. B. Hill 向英国的 6 万名医生发出调查问卷，并跟踪他们的吸烟情况和健康情况。根据跟踪情况，Doll 和 Hill 提交了两份报告（统称为英国医生研究）：1954 年提交的 The Mortality of Doctors in Relation to Their Smoking Habits 和 1956 年提交的 Lung Cancer and Other Causes of Death in Relation to Smoking。报告中公布了吸烟与疾病发病率的确凿科学证据，报告一经发布就震惊了世界，就连 Doll 本人也吓得戒掉了自己的重度烟瘾。

针对英国医生的跟踪调查一直持续到 2001 年，当时参与问卷的大多数参与者已经故去，并且还在世的最年轻的参与者也已经七十多岁了。在近 50 年的研究中，研究人员公布了多项报告，这些报告强化了公众对吸烟会导致呼吸道疾病和心血管疾病的认知。

关于这项研究的内容可以参阅：The British Doctors’ Study (1951–2001)。

2.3.1 样本均值，样本中位数和样本众数

在这一节中，我们将介绍一些用于描述数据集中心（center）的统计量。首先，假设我们有一个由 \(n\) 个数据构成的数据集：\(x_1，x_2，...，x_n\) 。这 \(n\) 个数的算术平均数（arithmetic average）就是样本均值（sample mean）。

定义 2.1 样本均值（sample mean） \(\overline{x}\) 的定义如下：

\[ \overline{x} = \frac{\sum_{i=1}^{n}{x_i}}{n} \tag{2.1}\]

如果 \(a\)，\(b\) 为常数，并且

\(y_i = ax_i+b, i\in[1,n]\)

那么，数据集 \(y_1,...,y_n\) 的样本均值为：

\[ \begin{align} \overline{y}&=\frac{\sum_{i=1}^{n}{ax_i+b}}{n} \\ &= \frac{\sum_{i=1}^{n}{ax_i}}{n} + \frac{\sum_{i=1}^{n}{b}}{n} \\ &= a\overline{x}+b \end{align} \tag{2.2}\]

习题 2.1 2004 年至 2013 年，美国高尔夫球大师赛的获胜分数如下：

\(280, 278, 272, 276, 281, 279, 276, 281, 289, 280\)

请计算如上分数的 样本均值。

解 2.1. 比起直接将这些分数相加，对于每个分数而言，先减去 280 以获得新值 \(y_i = x_i - 280\) 会更容易计算其样本均值：

\(0, −2, −8, −4, 1, −1, −4, 1, 9, 0\)

\(y_i\) 的算术平均数为 \(-\frac{8}{10}\)，因此 \(\overline{y}=-\frac{8}{10}\)，所以 \(\overline{x}=\overline{y}+280=279.2\)。

有时，我们想要确定一个用频率表来描述的数据集的样本均值，在频率表中存在 \(k\) 个不同的值 \(v_1,...,v_k\)，以及它们对应的频率 \(f_1,...,f_k\)。该数据集由 \(n=\sum_{i=1}^{k}{f_i}\) 个值组成，其中值 \(v_i\) 出现 \(f_i\) 次，因此这 \(n\) 个数据的样本均值为：

\[ \overline{x}=\sum_{i=1}^{k}{\frac{v_if_i}{n}} \tag{2.3}\]

对上述计算进行展开，得到：

\(\overline{x}=\frac{f_1}{n}{v_1} + \frac{f_2}{n}{v_2} + ... + \frac{f_k}{n}{v_k}\)

因此，样本均值是数据集中不同值的加权平均数，其中 \(v_i\) 的权重等于 \(n\) 个数据中值等于 \(v_i\) 的数据的占比。

习题 2.2 如下的频率表给出了青少年交响乐团的成员年龄分布。

代码 
library(knitr)

df <-data.frame(
  Age = c(15, 16, 17, 18, 19, 20),
  Frequency = c(2, 5, 11, 9, 14, 13)
)

kable(df, align = "l")
Age Frequency
15 2
16 5
17 11
18 9
19 14
20 13

计算该交响乐团中 54 位成员的年龄的样本均值。

解 2.2. \(\overline{x} = (15·2 + 16·5 + 17·11 + 18·9 + 19·14 + 20·13) / 54 ≈ 18.24\)

另一个用于表示数据集中心的统计量是样本中位数（sample median），简单来说，中位数就是当数据集中的数据按递增顺序排列后，位于数据集中间的数。

定义 2.2 把大小为 \(n\) 的数据集中的数据从小到大排序。如果 \(n\) 是奇数，则 样本中位数 就是第 \(\frac{(n + 1)}{2}\) 个数；如果 \(n\) 是偶数，则 样本中位数 就是第 \(\frac{n}{2}\) 个数和第 \(\frac{n}{2} + 1\) 个数的平均值。

因此，3 个数的样本中位数是其第 2 小的数；4 个数的样本中位数是第 2 小和第 3 小的数的平均值。

习题 2.3 找出 习题 2.2 中的样本中位数。

解 2.3. 由于有 54 个数，因此当数据按递增顺序排列后，样本中位数是第 27 个数和第 28 个数的平均值。因此，样本中位数是 18.5。

样本均值和样本中位数都是描述数据集中心趋势的、非常有用的统计量。样本均值的计算会用到数据集中的所有数据，并且容易受到数据集中极端值（远远大于或小于数据集中大多数数据的数据）的影响。样本中位数的计算只会到数据集中的一个或两个数据，因此样本中位数不受数据集中的极端值的影响。样本均值更好还是样本中位数更好，取决于我们试图从数据中学习什么。

  • 如果政府有统一的个人所得税税率，并试图估计税收总额，那么可以使用居民收入的样本均值作为统计量指标。
  • 如果政府正在考虑建造保障房，并希望确定能够负担得起这种住房的居民人口比例，那么更好的方式则是使用样本中位数。

习题 2.4 1972 年，Hoel, D. G. 发表了一篇论文 A representation of mor- tality data by competing risks。在这篇论文中，Hoel 首先让一组 5 周大的小白鼠接受 300rad 的辐射剂量，然后把这些小白鼠分成两组：第一组置于无菌环境中，第二组则置于常规的实验室环境。随后，Hoel 观察并记录这些小白鼠的生存天数。以下的茎叶图（茎以百天为单位）记录了因胸腺淋巴瘤死亡的小白鼠的数据：第一张图是生活在无菌条件下的小白鼠，第二张图是生活在普通实验室条件下的小白鼠。

表 2.6: 不同环境下的小白鼠实验
无菌环境下的小白鼠
Stem Leaf
1 58, 92, 93, 94, 95
2 02, 12, 15, 29, 30, 37, 40, 44, 47, 59
3 01, 01, 21, 37
4 15, 34, 44, 85, 96
5 29,37
6 24
7 07
8 00
正常实验室环境下的小白鼠
Stem Leaf
1 59, 89, 91, 98
2 35, 45, 50, 56, 61, 65, 66, 80
3 43, 56, 83
4 03, 14, 28, 32

计算两组小白鼠的生存时间的样本均值和样本中位数。

解 2.4. 从茎叶图中可以清晰地看出，无菌环境下小白鼠的样本均值（344.07）大于常规实验室环境下小白鼠的样本均值（292.32）。另一方面，由于无菌环境下的小白鼠有 29 个数据值，所以其样本中位数是第 15-最大 的数据值，即 259。另一组小白鼠的样本中位数是第 10-最大 的数据值，即 265。因此，尽管第一组数据的样本均值大于第二组，但这两组数据的样本中位数却大致基本一致。其中的原因在于，第一组数据中大于 500 的 5 个数据对其样本均值的影响较大，但对其样本中位数的影响要小得多。事实上，如果把大于 500 的 5 个数据替换为任何其他五个大于或等于 259 的数据，其样本中位数将保持不变。从茎叶图来看，无菌环境可能延长了五只寿命最长的小白鼠的寿命，但无法确定无菌环境是否对其他小白鼠的寿命产生了什么别的影响。

另一个用于表示数据集中心趋势的统计量是 样本众数（sample mode）。在数据集中，出现频率最高的值为 样本众数。如果出现频率最高的值有多个，那么所有出现频率最高的这些数值都是 样本众数。

习题 2.5 以下的频率表为抛 40 次骰子获得的结果。

代码 
library(knitr)

df <- data.frame(
    Value = c(1, 2, 3, 4, 5, 6),
    Frequency = c(9, 8, 5, 5, 6, 7)
)
kable(df, align = "l")
Value Frequency
1 9
2 8
3 5
4 5
5 6
6 7
  • 计算其样本均值
  • 计算其样本中位数
  • 计算器样本众数

解 2.5.

  • 样本均值为 \(\overline{x}=(9 + 16 + 15 + 20 + 30 + 42) / 40 = 3.05\)
  • 样本中位数为 第 20-最小 和 第 21-最小 的平均数 3
  • 样本众数为 1

2.3.2 样本方差和样本标准差

虽然我们已经介绍了用于描述数据集中心趋势的统计量，但我们也对描述数据波动的统计量感兴趣。可以通过计算数据值与样本均值之间距离的平方的平均值来描述数据集中数据的波动，这就是样本方差（sample variance）。出于技术原因，在计算样本方差时，需要除以 \(n-1\) 而不是 \(n\)（\(n\) 是数据集的大小）。

定义 2.3 对于数据集 \({x_1,...,x_n}\)，其 样本方差 \(s^2\) 的定义如下：

\[ s^2 = \frac{\sum_{i=1}^{n}{(x_i-\overline{x})^2}}{n - 1} \tag{2.4}\]

习题 2.6 计算如下两个数据集的样本方差：

  • A: 3, 4, 6, 7, 10
  • B: -20, 5, 15, 24

解 2.6. 对于 A 而言，其样本均值 \(\overline{x} = (3 + 4 + 6 + 7 + 10) / 5 = 6\)，因此，其样本方差为： \(s^2 = [(−3)^2 + (−2)^2 + 0^2 + 1^2 + 4^2] / 4 = 7.5\)

对于 B 而言，其样本均值也是 6，因此，其样本方差为： \(s^2 = [(−26)^2 + (−1)^2 + 9^2 + (18)^2] / 3 ≈ 360.67\)

通过 Solution 2.6 可发现，虽然 A 和 B 的样本均值是一样的，但是 B 的样本方差远大于 A。

在计算样本方差时，会经常用到如下的数学等式：

\[ \sum_{i=1}^{n}{(x_i-\overline{x})^2} = \sum_{i=1}^{n}{x_i^2} - n\overline{x}^2 \tag{2.5}\]

式 2.5 的证明过程如下：

\[ \begin{align} \sum_{i=1}^{n}{(x_i-\overline{x})^2} &= \sum_{i=1}^{n}{(x_i^2 - 2x_i\overline{x} + \overline{x}^2)} \\ &=\sum_{i=1}^{n}{x_i^2} - 2\overline{x}\sum_{i=1}^{n}{x_i} + \sum_{i=1}^{n}{\overline{x}^2} \\ &=\sum_{i=1}^{n}{x_i^2} - 2n\overline{x}^2 + n\overline{x}^2 \\ &=\sum_{i=1}^{n}{x_i^2} - n\overline{x}^2 \end{align} \]

如果 \(y_i = ax_i + b\)，\(i = 1,...,n\)，则 \(\overline{y} = a\overline{x} + b\)，因此：

\(\sum_{i=1}^{n}{(y_i - \overline{y})^2} = a^2\sum_{i=1}^{n}{(x_i-\overline{x})^2}\)

所以：

\[ s_y^2 = a^2s_x^2 \tag{2.6}\]

换句话说，为数据集中的每个数据增加一个常数不会改变其样本方差；而让数据集中的每个数据乘以一个常数则会得到一个新的样本方差，新的样本方差等于原来的样本方差乘以该常数的平方。

习题 2.7 下表给出了 1997 年至 2005 年间全球商业航空运输中导致人员死亡的事故数量。

代码 
library(knitr)

df <- data.frame(
    Year=seq(1997, 2005, by=1),
    Accidents = c(25, 20, 21, 18, 13, 13, 7, 9, 18))

kable(df, align="l")
Year Accidents
1997 25
1998 20
1999 21
2000 18
2001 13
2002 13
2003 7
2004 9
2005 18

根据如上的数据，计算航工事故的样本方差。

解 2.7. 将数据集中的数据都减去 18，我们得到一个新的数据集：

\(y = \{7,2,3,0,−5,−5,−11,−9,0\}\)

\(\overline{y} = \sum_{i=1}^{9}{y_i} / 9 = -2\)

\(\sum_{i=1}^{9}{y_i^2} = 49 + 4 + 9 + 25 + 25 + 121 + 81 = 314\)

因为新的数据集的方差和原数据集是一样的，因此根据 式 2.5：

\(s^2 = \frac{\sum_{i=1}^{9}{y_i^2} - n\overline{x}^2}{n - 1} = \frac{314 - 9·4}{8} = 34.75\)

注记

因为 \(y_i = x_i - 20\)，因此，根据 式 2.6， \(y\) 和 \(x\) 的方差是一样的。

样本方差的平方根称之为 样本标准差（sample standard deviation）。

定义 2.4 样本标准差 \(s\)的定义如下所示：

\[ s = \sqrt{\frac{\sum_{i=1}^{n}{(x_i-\overline{x})^2}}{n - 1}} \tag{2.7}\]

2.3.3 样本百分位数和箱线图

简单来说，数据集的 样本\(100p\)-百分位数 是指小于或等于数据集中 \(100p\%\) 的数据所对应的数据值，其中 \(0 \le p \le 1\)。

定义 2.5 样本\(100p\)-百分位数 是数据集中的一个数据值，这个值满足至少 \(100p\%\) 的数据小于或等于它，并且至少 \(100(1-p)\%\) 的数据大于或等于它。如果存在两个数据值满足这个条件，那么该数据集的 样本\(100p\)-百分位数 是这两个值的算术平均数。

为了确定大小为 \(n\) 的数据集的样本\(100p\)-百分位数，我们需要确定这样的数据值，以使得该值满足如下的条件：

  1. 至少有 \(np\) 个数据小于或等于该值
  2. 至少有 \(n(1-p)\) 个数据大于或等于该值

为了计算样本\(100p\)-百分位数，首先要对数据集进行升序排列。

  • 如果 \(np\) 不是整数，那么唯一满足前述条件的数据值是数据从小到大排序后，位置大于 \(np\) 的最小整数所对应的值。例如，如果 \(n=22\)，\(p=0.8\)，那么我们需要一个数据值，使得至少有 17.6 个数据小于或等于它，并且至少有 4.4 个数据大于或等于它。显然，只有第 18 小的值同时满足这两个条件，这就是 样本80-百分位数。
  • 如果 \(np\) 是整数，那么很容易发现 \(np\) 和 \(np+1\) 这两个位置的数据都满足前述条件，因此样本 \(100p\)-百分位数就是这两个数据值的算术平均数。例如，如果我们想计算大小为 20 的数据集的 90-百分位数，那么第 18-小 和第 19-小 的值都会使得至少有 90% 的数据小于或等于它们，并且至少有 10% 的数据大于或等于它们。因此，90-百分位数是这两个值的平均数。

习题 2.8 表 2.7 列出了 2006 年美国人口最多的 25 个城市的人口数。对于这组数据，找出：

  • 样本10-百分位数
  • 样本80-百分位数。
代码 
library(knitr)
df <- read.table("../misc/US_population.csv", header=TRUE, sep=",")
kable(df, align="l")
表 2.7: 2006 年 7 月，美国 25 个最大城市的人口数据
Rank City Population
1 New York, NY 8250567
2 Los Angeles, CA 3849378
3 Chicago, IL 2833321
4 Houston, TX 2144491
5 Phoenix, AR 1512986
6 Philadelphia, PA 1448394
7 San Antonio, TX 1296682
8 San Diego, CA 1256951
9 Dallas, TX 1232940
10 San Jose, CA 929936
11 Detroit, MI 918849
12 Jacksonville, FL 794555
13 Indianapolis, IN 785597
14 San Francisco, CA 744041
15 Columbus, OH 733203
16 Austin, TX 709893
17 Memphis, TN 670902
18 Fort Worth, TX 653320
19 Baltimore, MD 640961
20 Charlotte, NC 630478
21 El Paso, TX 609415
22 Milwaukee, WI 602782
23 Boston, MA 590763
24 Seattle, WA 582454
25 Washington, DC 581530

解 2.8.

  • 因为样本的大小是 25，25·0.1 = 2.5，因此，样本10-百分位数就是样本第 3-小的数据，也就是 590763。
  • 因为 25·0.8 = 20，因此，样本80-百分位数是第 20-小和第 21-小的算术平均数，\(\frac{1512986 + 1448394}{2}\)，也就是 1480690。

样本50-百分位数，就是样本中位数。样本25-百分位数、样本50-百分位数、样本75-百分位数一起构成了样本四分位数（sample quartiles）。

定义 2.6 样本25-百分位数称之为 第一四分位数（the first quartile），样本50-百分位数称之为 样本中位数 或 第二四分位数（the second quartile），样本75-百分位数称之为 第三四分位数（the third quartile）

注记

有时，也会把样本四分位数称为下四分位数和上四分位数：下四分位数（25-百分位数），中位数（50-百分位数），上四分位数（75-百分位数）。

四分位数把数据集分为四个部分，大约 25% 的数据小于第一四分位数，25% 的数据位于第一四分位数和第二四分位数之间，25% 的数据位于第二四分位数和第三四分位数之间，剩余的 25% 的数据大于第三四分位数。

习题 2.9 我们使用分贝（dB）来测量噪音。听力正常的人在安静环境下能听到的最弱的声音大约是 1dB，窃窃私语时的声音大约是 30dB，人们正常交谈时的声音大约是 70dB，收音机的音量大约是 100dB。当声音超过 120 dB 时，耳朵就开始感知到不适了。

曼哈顿中央车站外测量到的 36 个不同时间点的噪音水平如下所示：

曼哈顿中央车站噪音水平
82, 89, 94, 110, 74, 122, 112, 95, 100, 78, 65, 60, 90, 83, 87, 75, 114, 85 69, 94, 124, 115, 107, 88, 97, 74, 72, 68, 83, 91, 90, 102, 77, 125, 108, 65

计算如上数据的四分位数。

解 2.9. 因为 \(n = 36\)，所以四分位数对应的 \(np\) 均为整数，因此其对应的四分位数值为第 \(np\)-小和第 \((np+1)\)-小的平均数。

  • 第一四分位数（\(p=0.25\)）：第 9-小和第 10-小的平均值，也就是 74.5。
  • 第二四分位数（\(p=0.5\)）：第 18-小和第 19-小的平均值，也就是 89.5。
  • 第三四分位数（\(p=0.75\)）：第 27-小和第 28-小的平均值，也就是 104.5。

我们经常用箱线图（box plot）绘制数据集的汇总统计量。

  1. 在横坐标轴上绘制一条从最小值到最大值的直线段。
  2. 在这条直线上叠加一个 “箱子”，该箱子从第一四分位数开始一直延伸到第三四分位数，并用垂直线表示第二四分位数。

例如，表 2.1 中给出的 42 个起薪数据覆盖了从最低值 57000$ 到最高值 70000$ 之间的数。第一四分位数（第 11-小的值）是 60000\(，第二四分位数（第 21-小和第 22-小的平均值）是 61500\)，第三四分位数（第 32-小的值）是 64000$。该数据集的箱线图如 图 2.7 所示。

代码 
library(ggplot2)
ss <- c(57, 58, 59, 60, 61, 62, 63, 64, 66, 67, 70)
f <- c(4, 1, 3, 5, 8, 10, 0, 5, 2, 3, 1)
df <- data.frame(value=rep(ss, f))
ggplot(df, aes(x=value, y="")) +
  geom_boxplot() + 
  labs(x="Starting Salary", y="")
图 2.7: 起薪线图
注记

图 2.7 中的箱线图和文中对各四分位数的计算并不一致，这主要是因为该图使用 R 的 ggplot2 包绘制，对于该包而言，在绘制箱线图时，会剔除数据中可能的离群点（\(\pm 1.5IQR\)），离群点在图中显示为一个点。

箱线图上直线的长度等于数据集中的最大值减去数据集中的最小值，称为数据范围。箱线图中箱子的长度等于第三四分位数减去第一四分位数，称为四分位距（IQR: inter quaritle range）。

2.4 切比雪夫不等式

设 \(\overline {x}\) 和 \(s\) 分别是数据集的样本均值和样本标准差。假设 \(s > 0\)，切比雪夫不等式表明，对于 \(\forall k \ge 1\) ，至少有 \((100 \cdot (1−\frac {1}{k^2})) \%\) 的数据位于 \([ \overline {x} − ks, \overline {x} + ks]\) 的区间内。

  • 令 \(k = \frac {3}{2}\)，根据切比雪夫不等式，至少有 \(55.56\%\) 的数据会位于 \(\overline {x} \pm 1.5s\) 的范围内。
  • 令 \(k=2\)，则至少有 \(75\%\) 的数据位于 \(\overline {x} \pm 2s\) 的范围内。
  • 令 \(k=3\)，则至少有 \(88.9\%\) 的数据位于 \(\overline {x} \pm 3s\) 的范围内。

2.4.1 切比雪夫不等式

对于数据集 \(\{x_1,...,x_n\}\) 而言，令 \(\overline{x}\) 和 \(s\) 分别是其样本均值和样本标准差，其中 \(s \gt 0\)。

令 \(S_k=\{|x_i - \overline{x}| \lt ks, i \in [1, n]\}\)，即 \(S_k\) 为距离 \(\overline{x}\) \(k\) 个标准差以内的数据构成的集合。

令 \(|S_k|\) 为集合 \(S_k\) 中的元素个数。

则有：

\[ \frac{|S_k|}{n} \ge {1 - \frac{n-1}{nk^2}} \gt {1 - \frac{1}{k^2}}, \ \ \forall k \ge 1 \tag{2.8}\]

式 2.8 就是切比雪夫不等式，其证明过程如下：

\[ \begin{align} (n-1)s^2 &= \sum_{i=1}^{n}{(x_i-\overline{x})^2} \\ &= \sum_{i \in S_k}{{(x_i-\overline{x})^2}} + \sum_{i \notin S_k}{{(x_i-\overline{x})^2}} \\ & \ge \sum_{i \notin S_k}{{(x_i-\overline{x})^2}} \end{align} \]

\[ \begin{align} &\because \forall i \in S_k, \ \ |x_i - \overline{x}| \lt ks, \\ &\therefore \forall i \notin S_k, \ \ |x_i - \overline{x}| \ge ks, \\ &\therefore \sum_{i \notin S_k}{{(x_i-\overline{x})^2}} \ge \sum_{i \notin S_k}{k^2s^2} \\ &=k^2s^2(n-|S_k|) \end{align} \]

\[ \therefore (n-1)s^2 \ge k^2s^2(n-|S_k|) \tag{2.9}\]

当 \(s^2 \ne 0\) 时，对于 式 2.9 的等式两边都除以 \(nk^2s^2\)，得到：

\[ \begin{align} &\frac{n-1}{nk^2} \ge \frac{n-|S_k|}{n}, \\ &\therefore \frac{n-1}{nk^2} \ge 1-\frac{|S_k|}{n} \\ &\therefore \frac{|S_k|}{n} \ge 1-\frac{n}{nk^2}+\frac{1}{nk^2} \\ &\therefore \frac{|S_k|}{n} \ge 1-\frac{1}{k^2} \end{align} \]

由于切比雪夫不等式的通用性，对于给定的数据集，实际位于 \([\overline{x}-ks, \overline{x}+ks]\) 区间内的数据的百分比可能比切比雪夫不等式给出的下限要大一点。

例 2.4 表 2.8 给出了 2013 年 6 月美国最畅销的 10 款汽车。根据数据可知：\(\overline{x}=35.33\)，\(s=11.86\)，当 \(k=\frac{3}{2}\) 时，至少有 55.55%（\(100(1-\frac{1}{k^2})\%\)）的数据位于 [17.54, 53.12]（\([\overline{x}-ks, \overline{x}+ks]\)）区间内。而实际上，对于表中的数据而言，有 \(90\%\) 的数据落在了这个区间。

代码 
library(knitr)
df <- read.table("../misc/US_cars.csv", header=TRUE, sep=",")
kable(df, align="l")
表 2.8: 2013 年 6 月美国最畅销的 10 款汽车
Rank Model Sales.Volume.in.thousands.of.vehicles.
1 Ford F Series 68.0
2 Chevrolet Silverado 43.3
3 Toyota Camry 35.9
4 Chevrolet Cruze 32.9
5 Honda Accord 31.7
6 Honda Civic 29.7
7 Dodge Ram 29.6
8 Ford Escape 28.7
9 Nissan Altima 26.9
10 Honda CR-V 26.6

对于数据集 \(\{x_1,...,x_n\}\) 而言，令 \(\overline{x}\) 和 \(s\) 分别是其样本均值和样本标准差，其中 \(s \gt 0\)。那么，\(\{x_i - \overline{x} \ge ks, i \in [1, n]\}\) 的数据占比是多少呢？

假设令 \(N_k=\{x_i - \overline{x} \ge ks, i \in [1, n]\}\)，我们是否可以计算出 \(\frac{|N_k|}{n}\)？

当然！！！

根据 式 2.8，我们可以得到：

\[ \frac{|N_k|}{n} \le \frac{1}{k^2} \tag{2.10}\]

不过，利用接下来要介绍的单边切比雪夫不等式，我们可以得到一个更为精准的数据。

2.4.2 单边切比雪夫不等式

对于数据集 \(\{x_1,...,x_n\}\) 而言，令 \(\overline{x}\) 和 \(s\) 分别是其样本均值和样本标准差，其中 \(s \gt 0\)。

令 \(N_k=\{x_i - \overline{x} \ge ks, i \in [1, n]\}\)， \(|N_k|\) 为集合 \(N_k\) 中的元素个数，则有：

\[ \frac{|N_k|}{n} \le \frac{1}{1+k^2} \tag{2.11}\]

我们称 式 2.11 为单边切比雪夫不等式（one-sided Chebyshev inequality），其证明过程如下：

令 \(y_i=x_i-\overline{x}\)，\(i \in [1,n]\)，对于 \(\forall b \gt 0\)，有：

\[ \begin{align} \sum_{i=1}^{n}{(y_i+b)^2} &\ge \sum_{i:y_i \ge ks}{(y_i+b)^2} \\ & \ge \sum_{i:y_i \ge ks}{(ks+b)^2} \\ & = |N_k| \cdot (ks+b)^2 \end{align} \tag{2.12}\]

因为 \(|N_k|\) 为 \(y_i \ge ks\) 的元素个数，\(ks\) 和 \(b\) 都是正数，因此有： \[ \begin{align} \sum_{i=1}^{n}{(y_i+b)^2} &= \sum_{i=1}^{n}{(y_i^2 + 2by_i + b^2)} \\ &= \sum_{i=1}^{n}{y_i^2} + 2b\sum_{i=1}^{n}{y_i} + nb^2 \end{align} \]

\[ \begin{align} &\because y_i = x_i - \overline{x} \\ &\therefore \sum_{i=1}^{n}{y_i} = \sum_{i=1}^{n}{(x_i - \overline{x})} = \sum_{i=1}^{n}{x_i} - n\overline{x} = 0 \end{align} \]

所以有：

\[ \begin{align} \sum_{i=1}^{n}{(y_i+b)^2} &= \sum_{i=1}^{n}{y_i^2} + nb^2 \\ &=(n-1)s^2 + nb^2 \end{align} \tag{2.13}\]

根据 式 2.12 和 式 2.13，有：

\[ (n-1)s^2 + nb^2 \ge |N_k| \cdot (ks+b)^2 \]

进而得到：

\[ \frac{|N_k|}{n} \le \frac{s^2 + b^2}{(ks + b)^2} \tag{2.14}\]

因为 \(b\) 是任意大于 0 的数，因此我们可以令：

\[ b=\frac{s}{k} \tag{2.15}\]

把 式 2.15 代入 式 2.14 得到 式 2.11：

\[ \frac{|N_k|}{n} \le \frac{1}{k^2 + 1} \]

根据 式 2.10，我们可以发现，在一个数据集中，超过样本均值的 2 倍标准差的数据最多占比 25%。但是利用 式 2.11 所示的单边切比雪夫不等式，我们可以将这个比例精确到最多占比 20%。

2.5 正态分布数据集

我们在实践中观察到的许多大数据集都具有形状相似的直方图。通常，这些直方图在样本中位数处达到峰值，然后在中位数的两侧以钟形对称的形式下降。我们称这样的数据集为正态数据集，称其直方图为正态直方图。图 2.8 (a) 展示了一个正态数据集的直方图。

如果一个数据集的直方图接近正态直方图，我们说该数据集近似正态分布。例如，我们会说图 图 2.8 (b) 给出的直方图来自一个近似正态的数据集，而图 图 2.8 (c) 和图 图 2.8 (d) 给出的直方图则不是近似正态的数据集（因为每个直方图都太不对称）。与样本中位数不是近似对称的任何数据集都是偏态（skewed）分布。如果数据集的长尾在中位数右侧，则称为 “右偏”；如果数据集的长尾在中位数左侧，则称为 “左偏”。因此，图 图 2.8 (c) 中的数据集为左偏分布，而 图 2.8 (d) 中的数据集为右偏分布。

代码 
library(ggplot2)

x <- seq(1, 19)
y <- 100 - 10 * abs(x - mean(x))
df <- data.frame(value=rep(x, y))

ggplot(df, aes(x=value)) +
  geom_histogram(bins=19)


x <- seq(1, 19)
y <- 100 - 10 * abs(x - mean(x))
value <- rep(x, y)
value <- value[value != 17]
df <- data.frame(value=value)

ggplot(df, aes(x=value)) +
  geom_histogram(bins=19)

df <- data.frame(value=rbeta(10000, 5, 2))
ggplot(df, aes(x=value)) +
  geom_histogram(bins=19)

df <- data.frame(value=rbeta(10000, 1, 5))
ggplot(df, aes(x=value)) +
  geom_histogram(bins=19)
(a) 正态分布
(b) 近似正态分布
(c) 左偏分布
(d) 右偏分布
图 2.8: 不同数据集的直方图

从正态直方图的对称性可以看出，一个近似正态分布的数据集，其样本均值和样本中位数大致相等。

假设 \(\overline{x}\) 和 \(s\) 分别是近似正态分布数据集的样本均值和样本标准差。经验法则（empirical rule）指定了观察到的数据位于 \(s\)、\(2s\)、\(3s\) 个样本均值 \(x\overline{x}\) 内的大致比例。

2.5.1 经验法则

如果一个数据集是一个近似正态分布的数据集，其样本均值和样本标准差分别是 \(\overline{x}\) 和 \(s\)，那么有如下的正确结论：

  1. 大约有 68% 的观察数据会位于 \(\overline{x} \pm s\) 区间内
  2. 大约有 95% 的观察数据会位于 \(\overline{x} \pm 2s\) 区间内
  3. 大约有 99% 的观察数据会位于 \(\overline{x} \pm 3s\) 区间内

如下的数据是工业工程专业的学生在统计学考试中的分数：

43, 46, 52, 55, 55, 56, 58, 60, 62, 63, 64, 66, 66, 72, 74, 74, 75, 77, 77, 78, 83, 85, 85, 87, 88, 90, 91, 94

通过其茎叶图，我们会发现如上数据的直方图近似正态分布。接下来我们用这份数据来评估一下经验法则。

代码 
scores <- c(43, 46, 52, 55, 55, 56, 58, 60, 62, 63, 64, 66, 66, 72, 74, 74, 75, 77, 77, 78, 83, 85, 85, 87, 88, 90, 91, 94)
stem(scores)

  The decimal point is 1 digit(s) to the right of the |

  4 | 36
  5 | 25568
  6 | 023466
  7 | 2445778
  8 | 35578
  9 | 014

解 2.10. 计算得到数据集的样本均值和样本方差：\(\overline{x} \thickapprox 70.571\)，\(s \thickapprox 14.354\)。因此，根据经验法则：

  • 有 68% 的数据会位于 56.2~84.9 之间，但是实际上，数据集中只有 \(\frac{15}{28} \%\) 也就是 53.6% 的数据位于56.2~84.9 之间。
  • 有 95% 的数据会位于 41.68~99.28 之间，但是实际上，数据集中所有的数据都位于 41.68~99.28 之间。

如果一个总体是由多个不同类型的子总体而构成，那么从这个总体中抽样得到的数据集通常不是正态分布。这种数据集的直方图通常看起来像是多个正态直方图的组合或叠加，因此在这种直方图上通常会有多个局部峰值。由于这些局部峰值处的直方图比其邻近值更高，因此这些峰值类似于众数。如 图 2.9 所示，在直方图中，具有两个局部峰值的数据集称之为双峰分布。

代码 
samples1 <- rnorm(10000, mean=1, sd=1)  
samples2 <- rnorm(10000, mean=5, sd=1)
df <- data.frame(value=c(samples1, samples2))

ggplot(df, aes(x=value)) +
    geom_histogram()
图 2.9: 双峰分布直方图

2.6 成对数据集和样本相关系数

我们经常关注由成对数据值构成的数据集，这些数据集中的成对数据之间存在着某种关系。在这样的数据集中，如果每个成对数据都有一个 \(x\) 值和一个 \(y\) 值，那么我们就用 \((x_i, y_i)\) 表示第 \(i\) 个数据。例如，为了确定每天中午的温度（以摄氏度为单位）与当天生产的不合格零件数量之间的关系，一家公司记录了 表 2.9 中的数据。对于这个数据集，\(x_i\) 代表第 \(i\) 天中午的摄氏温度，\(y_i\) 代表第 \(i\) 天生产的不合格零件数量。

代码 
library(knitr)
df <- read.table("../misc/temp_defect.csv", header=TRUE, sep=",")
kable(df, align="l")
表 2.9: 天气温度和产品缺陷量数据
Day Temperature Number.of.Defects
1 24.2 25
2 22.7 31
3 30.5 36
4 28.6 33
5 25.5 19
6 32.0 24
7 28.6 27
8 26.5 25
9 25.3 16
10 26.0 14
11 24.4 22
12 24.8 23
13 20.6 20
14 25.1 25
15 21.4 25
16 23.7 23
17 23.9 27
18 25.2 30
19 27.4 33
20 28.3 32
21 28.8 35
22 26.6 24

在二维图上绘制成对数据集中的数据是一种表示成对数据集的有效方法，其中 \(x\) 轴代表成对数据的 \(x\) 值，\(y\) 轴代表成对数据的 \(y\) 值，这样的图称之为散点图（scatter diagram）。图 2.10 为 表 2.9 中数据的散点图。

代码 
library(ggplot2)
df <- read.table("../misc/temp_defect.csv", header=TRUE, sep=",")
ggplot(df, aes(x=Temperature, y=Number.of.Defects)) + 
    geom_point()
图 2.10: 天气温度和产品缺陷量散点图

对于成对数据集而言，一个有趣问题是，比较大的 \(x\) 值是否倾向于与比较大的 \(y\) 值配对，小的 \(x\) 值是否与小的 \(y\) 值配对？如果情况并非如此，那么我们可能会质疑其中一个变量的较大的值是否倾向于与另一个变量的较小的值配对。通常，可以通过散点图来大致回答这些问题。例如，图 2.10 表明，高温和缺陷产品数量高之间似乎存在某种联系。为了对成对数据之间的这种关系进行定量度量，我们需要开发一个新的统计量，以衡量 \(x\) 值与 \(y\) 值之间的配对程度。

假设数据集由 \(\{(x_i, y_i)\}\) 组成，其中 \(i = 1,...,n\)，\(\overline{x}\) 和 \(\overline{y}\) 分别是 \(x\) 和 \(y\) 的样本均值。对于第 \(i\) 对数据，使用 \(x_i - \overline{x}\) 作为 \(x_i\) 与其样本均值的偏差，使用 \(y_i - \overline{y}\) 作为 \(y_i\) 与其样本均值的偏差。如果 \(x_i\) 的值比较大，那么 \(x_i \gt \overline{x}\)，所以 \(x_i - \overline{x} \gt 0\)。类似地， 如果 \(x_i\) 的值比较小，那么 \(x_i \lt \overline{x}\)，所以 \(x_i - \overline{x} \lt 0\)。对于 \(y_i\) 而言，同样如此。于是，我们可以得出以下结论：

重要

当 \(x\) 的较大值趋向于和 \(y\) 的较大值相关联，\(x\) 的较小值趋向于和 \(y\) 的较小值相关联时，则 \(x_i - \overline{x}\) 和 \(y_i- \overline{y}\) 的符号（无论正负）都将趋于相同。

如果 \(x_i - \overline{x}\) 和 \(y_i - \overline{y}\) 的符号相同（无论正负），那么它们的乘积 \((x_i - \overline{x})(\)y_i - )$ 是正数。因此，当 \(x\) 的较大值趋向于和 \(y\) 的较大值相关联，\(x\) 的较小值趋向于和 \(y\) 的较小值相关联时，\(\sum_{i=1}^{n}{(x_i - \overline{x})(y_i - \overline{y})}\) 将趋向于一个较大的正数。事实上，当较大（小）的 \(x\) 与较大（小）的 \(y\) 配对时，不但所有的 \((x_i - \overline{x})(y_i - \overline{y})\) 的符号都是正的，而且当 \((x_i - \overline{x})\) 的最大值与 \((y_i - \overline{y})\) 的最大值配对、其对应的次大值配对、依次类推对 \((x_i - \overline{x})\) 和 \((y_i - \overline{y})\) 配对时，\(\sum_{i=1}^{n}{(x_i - \overline{x})(y_i - \overline{y})}\) 将获得最大值。同理，当 \(x_i\) 的较大值倾向于和 \(y_i\) 的较小值配对时，\(x_i - \overline{x}\) 和 \(y_i - \overline{y}\) 的符号相反，所以 \(\sum_{i=1}^{n}{(x_i - \overline{x})(y_i - \overline{y})}\) 将是一个较大的负数。

定义 2.7 对于成对数据集 \(\{(x_i, y_i), i = 1,...,n\}\)，\(s_x\) 和 \(s_y\) 分别表示 \(x\) 和 \(y\) 的样本标准差，则 \((x_i, y_i)\) 的样本相关系数（sample correlation coefficient）\(r\) 的定义如下：

\[ \begin{align} r&=\frac{\sum_{i=1}^{n}{(x_i-\overline{x})(y_i-\overline{y})}}{(n-1)s_xs_y} \\ &=\frac{\sum_{i=1}^{n}{(x_i-\overline{x})(y_i-\overline{y})}}{\sqrt{\sum_{i=1}^{n}{(x_i-\overline{x})^2} \cdot \sum_{i=1}^{n}{(y_i-\overline{y})^2}}} \end{align} \tag{2.16}\]

当 \(r \gt 0\) 时，成对数据集中的样本成正相关（positively correlated），当 \(r \lt 0\) 时，样本成负相关（negatively correlated）。

成对数据集的样本相关系数 \(r\) 具备如下的特性：

  1. \(-1 \le r \le 1\)
  2. \(\forall a>0\)， 如果 \(y_i = ax_i + b\)，则 \(r=1\)
  3. \(\forall a<0\)， 如果 \(y_i = ax_i + b\)，则 \(r=-1\)
  4. \(\forall a\) 和 \(\forall c\)，且 \(a \cdot c \gt 0\)，如果 \(r\) 是 \(\{(x_i, y_i)\}\) 的相关系数，则 \(r\) 也是 \(\{(ax_i + b, cy_i + d)\}\) 的相关系数
  • 特性 1 说明样本相关系数 \(r\) 总是介于 -1 和 +1 之间。
  • 特性 2 说明当成对数据之间存在线性关系并且 \(x\) 的较大（小）值趋向于和 \(y\) 的较大（小）值相关联时，\(r\) 等于 +1。
  • 特性 3 说明当成对数据之间存在线性关系并且 \(x\) 的较大（小）值趋向于和 \(y\) 的较小（大）值相关联时，\(r\) 等于 -1。
  • 特性 4 指出，当给每个 \(x\)（或 \(y\)） 变量加上一个常数，或将每个 \(x\)（或 \(y\)） 变量都乘以一个正数时，\(r\) 的值保持不变。特性 4 意味着 \(r\) 不依赖于测量数据的单位。例如，无论身高数据的测量单位是英尺还是英寸，也无论体重数据的测量单位是磅还是千克，一个人的身高和体重之间的样本相关系数都是一致的。此外，如果成对数据中的一个值是温度，那么无论该值是华氏温度还是摄氏温度，样本相关系数都是相同的。

样本相关系数 \(r\) 的绝对值 \(|r|\) 用于表示 \(x\) 和 \(y\) 之间线性关系的强度。\(|r| = 1\) 意味着完美的线性关系，也就是说，可以用一条直线穿过 \(\{(xi,yi), i = 1, ..., n\}\) 中的所有数据点 。\(|r| \thickapprox 0.8\) 意味着相对较强的线性关系，此时虽然没有任何一条直线可以穿过所有的数据点，但存在一条直线 “接近” 穿过所有的数据点。\(|r| \thickapprox 0.3\) 意味着成对数据之间的线性关系相对较弱。

\(r\) 的符号给定了相关性的方向。当较小的 \(y\) 倾向于和较小的 \(x\) 配对，而较大的 \(y\) 倾向于和较大的 \(x\) 配对时，对于这样的线性关系，\(r \gt 0\)。而当较小的 \(y\) 倾向于和较大的 \(x\) 配对，而较大的 \(y\) 倾向于和较小的 \(x\) 配对时，对于这样的线性关系，\(r \lt 0\)。图 2.11 显示了具有不同 \(r\) 值的数据集的散点图。

代码 
library(MASS)  
library(ggplot2)

# 设置你想要的相关系数  
rho <- -0.5  # 例如，我们想要的相关系数是0.5  

# 创建一个协方差矩阵。假设两个变量的方差都是1（对于标准化数据）  
# 协方差矩阵的对角线元素是方差，非对角线元素是协方差（与相关系数相关）  
Sigma <- matrix(c(1, rho, rho, 1), nrow = 2)  
data <- mvrnorm(n = 100, mu = c(0, 0), Sigma = Sigma)
df1 <- as.data.frame(data)
colnames(df1) <- c("X", "Y")

rho <- 0   
Sigma <- matrix(c(1, rho, rho, 1), nrow = 2)  
data <- mvrnorm(n = 100, mu = c(0, 0), Sigma = Sigma)
df2 <- as.data.frame(data)
colnames(df2) <- c("X", "Y")

rho <- 0.9   
Sigma <- matrix(c(1, rho, rho, 1), nrow = 2)  
data <- mvrnorm(n = 100, mu = c(0, 0), Sigma = Sigma)
df3 <- as.data.frame(data)
colnames(df3) <- c("X", "Y")

rho <- 1
Sigma <- matrix(c(1, rho, rho, 1), nrow = 2)  
data <- mvrnorm(n = 100, mu = c(0, 0), Sigma = Sigma)
df4 <- as.data.frame(data)
colnames(df4) <- c("X", "Y")

ggplot(df1, aes(x=X, y=Y)) +
    geom_point()

ggplot(df2, aes(x=X, y=Y)) +
    geom_point()

ggplot(df3, aes(x=X, y=Y)) +
    geom_point()

ggplot(df4, aes(x=X, y=Y)) +
    geom_point()
(a) r = -0.5
(b) r = 0
(c) r = 0.9
(d) r = 1
图 2.11: 不同样本相关系数的数据散点图

习题 2.10 计算 表 2.9 所示的样本相关系数。

解 2.11.

代码 
df <- read.table("../misc/temp_defect.csv", header=TRUE, sep=",")
r <- cor(df$Temperature, df$Number.of.Defects)
print(r)
[1] 0.418944

这意味着气温和当天的缺陷产品数量之间存在着弱正相关性。

例 2.5 以下数据表给出了 10 个学生的静息心率（每分钟心跳次数）和其对应的受教育年限。图 2.12 展现了这些数据的散点图。这些数据的样本相关系数 \(r=−0.7638\)。负相关系数表明，对于这个数据集来说，高静息心率与低受教育年限紧密相关，而低静息心率则与高受教育年限紧密相关。

代码 
library(knitr)
person <- seq(1:10)
Years <- c(12, 16, 13, 18, 19, 12, 18, 19, 12, 14)
BPM <- c(73, 67, 74, 63, 73, 84, 60, 62, 76, 71)
df <- data.frame(Person=person, Years.of.School=Years, BPM=BPM)

kable(df, align="l")
Person Years.of.School BPM
1 12 73
2 16 67
3 13 74
4 18 63
5 19 73
6 12 84
7 18 60
8 19 62
9 12 76
10 14 71
代码 
library(ggplot2)
person <- seq(1:10)
Years <- c(12, 16, 13, 18, 19, 12, 18, 19, 12, 14)
BPM <- c(73, 67, 74, 63, 73, 84, 60, 62, 76, 71)
df <- data.frame(Person=person, Years.of.School=Years, BPM=BPM)

ggplot(df, aes(x=Years.of.School, y=BPM)) +
    geom_point()
图 2.12: 不同数据集的直方图
相关性衡量的是变量之间的关联关系，而不是因果关系

例 2.5 中的数据集仅考虑了 10 名学生，因此不足以让人对受教育年限和心率之间的关系得出任何确凿的结论。此外，即使有更大规模的数据集，并且受教育年限高低与其静息心率之间同样存在着较强的负相关性，我们也没有理由得出 多接受几年教育会直接降低一个人的心率 的结论。也就是说，尽管高受教育年限往往与较低的静息心率有关联，但这并不意味着多受几年教育是导致较低心率的直接原因。通常，对这种关联性的解释在于，存在一个与所考虑的这两个变量都有相关性的隐藏因素。

在 例 2.5 中，可能是受教育年限高的人更了解健康领域的最新发现，因此可能更了解锻炼和良好营养的重要性。亦或者，可能不是知识在起作用，而是受过更多教育的人所从事的工作会让他们有更多时间进行锻炼，同时也可以获取更多薪水以补充良好的营养。受教育年限和静息心率之间的强烈负相关性可能是由受教育年限以及其他潜在因素的综合结果。

接下来，我们将证明样本相关系数 \(r\) 的第一个特性：\(|r| \le 1\)，当且仅当所有的数据点都在一条直线上时，\(=\) 成立。

\[ \begin{align} \because &\sum{(\frac{x_i - \overline{x}}{s_x} - \frac{y_i - \overline{y}}{s_y})^2} \ge 0 \\ \therefore & \sum{\frac{(x_i - \overline{x})^2}{s_x}} + \sum{\frac{(x_i - \overline{x})^2}{s_x}} - 2 \sum{\frac{(x_i - \overline{x})(y_i - \overline{y})}{s_xs_y}} \\ \therefore & (n - 1) + (n - 1) -2(n - 1)r \ge 0 \\ \therefore & 2(n - 1)(1 - r) \ge 0 \\ \therefore & r \le 1 \end{align} \tag{2.17}\]

假设所有的数据点 \({(x_i, y_i), i=1,...,n}\) 都位于直线 \(y_i = ax_i + b;\ i=1,...n \ \& \ a>0\) 上，则：

\[ \begin{align} & s_y^2=a^2s_x^2 \\ & \overline{y}=a\overline{x} + b \\ \therefore & a = \frac{s_y}{s_x}, \ b = \overline{y} - \frac{s_y}{s_x}\overline{x} \end{align} \tag{2.18}\]

根据 式 2.17， 当且仅当 \(r = 1\) 时，\(\sum{(\frac{x_i - \overline{x}}{s_x} - \frac{y_i - \overline{y}}{s_y})^2} = 0\)。

也即当且仅当 \(r = 1\) 时，\(\frac{y_i - \overline{y}}{s_y} = \frac{x_i - \overline{x}}{s_x}\)。

所以，当且仅当 \(r = 1\) 时，\(y_i = \overline{y} - \frac{s_y}{s_x}\overline{x} + \frac{s_y}{s_x}x_i\)，把 式 2.18 代入得到 \(y_i=ax_i + b\)。

因此，当且仅当 \(r = 1\) 时，\((x_i, y_i)\) 的所有的点都位于直线 \(y_i = ax_i + b, a>0\)。

同理，我们可以证明，当且仅当 \(r = -1\) 时，\((x_i, y_i)\) 的所有的点都位于直线 \(y_i = ax_i + b, a<0\)。

2.7 洛伦兹曲线和基尼系数

洛伦兹曲线（Lorenz Curve）\(L(p), 0 \le p \le 1\)，是与群体中成员的收入相关的图表。\(L(p)\) 表示群体中，收入最低的 \(100p\%\) 的人的总收入在群体总收入中的收入占比。例如，\(L(0.3)\) 是收入最低的 30% 的人的总收入在群体中的收入占比。一般来说，假设群体中有 \(n\) 个人，其个体收入按升序排列为 \(x_1 \le x_2 \le x_3 ... \le x_n\)。因为 \(x_1 + ... + x_j\) 是收入最低的 \(j\) 个人的总收入，而 \(x_1 + ... + x_n\) 是该群体所有人的总收入。因此，收入最低的 \(100\frac{j}{n}\%\) 的人的收入占比就是：

\[ L(\frac{j}{n}) = \frac{x_1 + ... + x_j}{x_1 + ... + x_n},\ j = 1, ..., n \]

一般而言，令 \(L(0) = 0\)，并通过直线连接 \(\frac{j}{n}\) 和 \(\frac{j + 1}{n}\) 来绘制洛伦兹曲线。

例 2.6 假设我们想绘制 \(n=5\) 的洛伦兹曲线，其中该群体的收入分别是：9，7，22，5，17。该群体的收入按照升序排序后是：5，7，9，17，22。因为群体的总收入是 60，因此：

  • L(0.2) = 5 / 60
  • L(0.4) = 12 / 60
  • L(0.6) = 21 / 60
  • L(0.8) = 38 / 60
  • L(1) = 60 /60

将如上的点用直线连接起来就形成了洛伦兹曲线，如 图 2.13。

代码 
library(ggplot2)
x <- seq(0, 1, 0.2)
y <- c(0, 5, 7, 9, 17, 22)
y <- cumsum(y)
df <- data.frame(x=x, y=y, yf=(y / 60))
ggplot(df, aes(x=x, y=yf)) +
  geom_line() + 
  scale_x_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  scale_y_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  geom_point(color = "black", size = 1) + 
  labs(x="people proportion", y="incomes proportion")
图 2.13: 例 2.6 的洛伦兹曲线

也可以用洛伦茨曲线来展现群体中的个体财富，其中 \(L(p)\) 代表最贫穷的 \(100p%\) 的人所拥有的全部财富的比例。例如，如果全部人口由 1000 人组成，那么 \(L(0.22)\) 就是 220 个最贫穷的人所拥有的财富比例。

现在，当向集合中添加一个大于之前所有值的新值时，集合的平均值总是会增加。因此，如果群体成员的收入递增排序后是 \(x_1 \le x_2 \le x_3... \le x_n\)，那么 \(\forall j = 1,...,n\)，有：

\[ \begin{align} & \frac{x_1 + ... + x_j}{j} \le \frac{x_1 + ... + x_{j+1}}{j + 1} \le \frac{x_1 + ... + x_n}{n} \\ \therefore & \frac{x_1 + ... + x_j}{x_1 + ... + x_n} \le \frac{j}{n} \\ \therefore & L{(\frac{j}{n}}) \le \frac{j}{n} \end{align} \]

此外，可以验证，当且仅当所有人的收入都相等时，\(L(\frac{j}{n}) = \frac{j}{n}, j=1,...,n\)。因此，除非所有人的收入都相等，否则洛伦茨曲线总是位于 \((0,0)\) 和 \((1,1)\) 构成的直线之下。对于 例 2.6，图 2.14 正说明了这一点。

代码 
library(ggplot2)

x <- seq(0, 1, 0.2)
y <- c(0, 5, 7, 9, 17, 22)
y <- cumsum(y)
df <- data.frame(x=x, y=y, yf=(y / 60))
ggplot(df, aes(x=x, y=yf)) +
  geom_line() + 
  scale_x_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  scale_y_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  geom_point(color = "black", size = 1) + 
  geom_segment(x=df$x[1], y=df$yf[1], xend=df$x[length(df$x)], yend=df$yf[length(df$yf)], color="red") +
  labs(x="people proportion", y="incomes proportion")
图 2.14: 例 2.6 的包含直线的洛伦兹曲线

因此，当所有人的收入都相等时，洛伦茨曲线与从 \((0, 0)\) 到 \((1, 1)\) 的直线重合；而当群体中不同个体收入不相等时，洛伦茨曲线则位于该直线的下方。当群体成员的收入越不平等，\((0, 0)\) 到 \((1, 1)\) 的直线和洛伦茨曲线之间的区域则越大（如 图 2.15 中的阴影区域所示）。

代码 
library(ggplot2)
x <- seq(0, 1, 0.2)
y <- c(0, 5, 7, 9, 17, 22)
y <- cumsum(y)
df <- data.frame(x=x, y=y, yf=(y / 60))
ggplot(df, aes(x=x, y=yf)) +
  geom_line() + 
  scale_x_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  scale_y_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  geom_point(color = "black", size = 1) + 
  geom_polygon(aes(x=x, y=yf), fill="grey", alpha=0.7) +
  geom_segment(x=df$x[1], y=df$yf[1], xend=df$x[length(df$x)], yend=df$yf[length(df$yf)], color="red") + 
  annotate("text", x = 0.75, y = 0.25, label = "B", hjust = 1, vjust = 0, size = 6) +
  labs(x="people proportion", y="incomes proportion")
图 2.15: 直线[(0, 0), (1, 1)]和洛伦兹曲线之间的面积

可以用基尼系数（Gini Index）衡量 \((0, 0)\) 到 \((1, 1)\) 的直线和洛伦茨曲线之间的区域大小，进而衡量群体中成员的收入不平等性。基尼系数等于 图 2.15 中阴影区域面积与直线下方区域的面积之比。由于三角形的面积是底乘以高的一半，所以从 \((0, 0)\) 到 \((1, 1)\) 的直线下方的面积等于 \(\frac{1}{2}\)。

因此，基尼系数 \(G\) 的定义如下：

\[ G = \frac{L(p) 曲线和直线围成的区域面积}{\frac{1}{2}} \]

令洛伦兹曲线下方的区域为 \(B\)，则洛伦兹曲线和直线围成的面积就是直线下方的面积减去洛伦兹曲线下方的面积。因此，基尼系数

\[ G=\frac{\frac{1}{2} - B}{\frac{1}{2}} = 1 - 2B \tag{2.19}\]

习题 2.11 计算 例 2.6 中的数据的基尼系数。

代码 
x <- seq(0, 1, 0.2)
y <- c(0, 5, 7, 9, 17, 22)
y <- cumsum(y)
df <- data.frame(x=x, y=y, yf=(y / 60))
ggplot(df, aes(x=x, y=yf)) +
  geom_line() + 
  scale_x_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  scale_y_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  geom_point(color = "black", size = 1) + 
  geom_polygon(aes(x=x, y=yf), fill="grey", alpha=0.7) +
  geom_segment(x=df$x[1], y=df$yf[1], xend=df$x[length(df$x)], yend=df$yf[length(df$yf)], color="red") + 
  geom_segment(x=df$x[2], y=df$yf[2], xend=df$x[2], yend=0) +
  geom_segment(x=df$x[3], y=df$yf[3], xend=df$x[3], yend=0) +
  geom_segment(x=df$x[4], y=df$yf[4], xend=df$x[4], yend=0) +
  geom_segment(x=df$x[5], y=df$yf[5], xend=df$x[5], yend=0) +
  geom_segment(x=df$x[6], y=df$yf[6], xend=df$x[6], yend=0) +
  annotate("text", x = 0.15, y = 0.02, label = "B1", hjust = 1, vjust = 0, size = 3) +
  annotate("text", x = 0.3, y = 0.1, label = "B2", hjust = 1, vjust = 0, size = 3) +
  annotate("text", x = 0.5, y = 0.2, label = "B3", hjust = 1, vjust = 0, size = 3) +
  annotate("text", x = 0.7, y = 0.4, label = "B4", hjust = 1, vjust = 0, size = 3) +
  annotate("text", x = 0.9, y = 0.6, label = "B5", hjust = 1, vjust = 0, size = 3) +
  labs(x="people proportion", y="incomes proportion")
图 2.16: 求解 B 的面积

解 2.12. 为了计算 \(B\) 的面积（即洛伦兹曲线下的面积），如 图 2.16 所示，我们令 \(B = B1 + B2 + B3 + B4 + B5\)，其中 \(B1\) 是0~0.2 的洛伦兹曲线下的区域面积，\(B2\) 是 0.2~0.4 之间的区域面积，\(B3\) 是 0.4~0.6 之间的区域面积，\(B4\) 是 0.6~0.8 之间的区域面积；\(B5\) 是 0.8~1 之间的区域面积。现在，\(B1\) 是一个三角形的面积，其底为 0.2，高为 5/60，因此：

\[ B1 = \frac{1}{2} \cdot 0.2 \cdot \frac{5}{60} = \frac{5}{600} \]

对于 \(B2, B3, B4, B5\) 而言，这些区域的面积由两部区域构成：顶部的三角形的面积和底部的矩形面积。如下图所示，以 \(B2\) 为例，其面积由 \(B2-1\) 和 \(B2-2\) 两部分构成：

代码 
x <- seq(0, 1, 0.2)
y <- c(0, 5, 7, 9, 17, 22)
y <- cumsum(y)
df <- data.frame(x=x, y=y, yf=(y / 60))
ggplot(df, aes(x=x, y=yf)) +
  scale_x_continuous(limits=c(0, 0.6), breaks=seq(0, 0.6, 0.2)) +
  scale_y_continuous(limits=c(0, 0.25), breaks=seq(0, 0.25, 0.05)) +
  geom_segment(x=df$x[2], y=df$yf[2], xend=df$x[3], yend=df$yf[3]) + 
  geom_segment(x=df$x[2], y=df$yf[2], xend=df$x[2], yend=0) +
  geom_segment(x=df$x[3], y=df$yf[3], xend=df$x[3], yend=0) +
  geom_segment(x=df$x[2], y=df$yf[2], xend=df$x[3], yend=df$yf[2], linetype="dashed") +
  annotate("text", x = 0.3, y = 0.1, label = "B2-1", hjust = 1, vjust = 0, size = 3) +
  annotate("text", x = 0.3, y = 0.03, label = "B2-2", hjust = 1, vjust = 0, size = 3) +
  labs(x="people proportion", y="incomes proportion")

因此，所有的三角形的底的长度都是 0.2，其对应的高分别为 \(\frac{5}{60}, \frac{7}{60}, \frac{9}{60}, \frac{17}{60}, \frac{22}{60}\)。因此，所有的三角形的面积为：

\[ S_\triangle = \frac{1}{2} \cdot 0.2 \cdot \frac{5+7+9+17+22}{60} = 0.1 \]

四个矩形的长均是 0.2，其对应的高分别是 \(\frac{5}{60}, \frac{12}{60}, \frac{21}{60}, \frac{38}{60}\)。因此，所有的矩形的面积为：

\[ S_\square = 0.2 \cdot \frac{5 + 12 + 21 + 38}{60} \thickapprox 0.25333 \]

所以 \(B = S_\triangle + S_\square =\) \(0.1 + 0.25333 = 0.35333\)，故而 \(G = 1 - 2B =\) \(1 - 2 \cdot 0.35333 =\) \(0.29334\)。

一般而言，收入数据以升序进行排列，\(x_1 \le x_2 \le x_3 \le ... \le x_n\)，令 \(s_j = x_1 + ... + x_j, j=1,...,n\)，则洛伦兹曲线下方的所有三角形的底均是 \(\frac{1}{n}\)，且其对应的高分别为 \(\frac{x_1}{s_n}\)，\(\frac{x_2}{s_n}\)，……，\(\frac{x_n}{s_n}\)。

因此，所有三角形的面积为：

\[ S_\triangle=\frac{1}{2n}\frac{(x_1 + ... + x_n)}{s_n}=\frac{1}{2n} \]

所有的矩形的长均为 \(\frac{1}{n}\)，且其对应的高分别为 \(\frac{s_1}{s_n}\)，\(\frac{s_2}{s_n}\)，……，\(\frac{s_{n-1}}{s_n}\)。

因此，所有矩形的面积为：

\[ S_\square = \frac{s_1 + s_2 + ... + s_{n-1}}{ns_n} \]

所以：

\[ B = \frac{1}{2n} + \frac{s_1 + s_2 + ... + s_{n-1}}{ns_n} \tag{2.20}\]

当所有人的收入都相等时，基尼系数为 0，因此直线与洛伦兹曲线之间的面积为 0。另一个极端的情况就是，当群体的 \(n\) 个人中，只有一个人有收入，那么此时基尼系数达到最大值。对于第二种场景，洛伦兹曲线下的面积是一个三角形的面积，其底边长度为 \(\frac{1}{n}\)，高度为 1（如 图 2.17）。此时 \(B=\frac{1}{2n}\)，因此 \(G=1-\frac{1}{n}\)。

代码 
x <- seq(0, 1, 0.2)
y <- c(0, 0, 0, 0, 0, 60)
y <- cumsum(y)
df <- data.frame(x=x, y=y, yf=(y / 60))
ggplot(df, aes(x=x, y=yf)) +
  geom_line() + 
  scale_x_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  scale_y_continuous(limits=c(0, 1), breaks=seq(0, 1, 0.2)) +
  geom_point(color = "black", size = 1) + 
  geom_polygon(aes(x=x, y=yf), fill="grey", alpha=0.7) +
  geom_segment(x=df$x[1], y=df$yf[1], xend=df$x[length(df$x)], yend=df$yf[length(df$yf)], color="red") + 
  geom_segment(x=df$x[2], y=df$yf[2], xend=df$x[2], yend=0) +
  geom_segment(x=df$x[3], y=df$yf[3], xend=df$x[3], yend=0) +
  geom_segment(x=df$x[4], y=df$yf[4], xend=df$x[4], yend=0) +
  geom_segment(x=df$x[5], y=df$yf[5], xend=df$x[5], yend=0) +
  geom_segment(x=df$x[6], y=df$yf[6], xend=df$x[6], yend=0) +
  labs(x="people proportion", y="incomes proportion") +
  theme(  
    axis.text.x = element_blank(),   
    axis.text.y = element_blank(), 
  ) + 
  annotate("text", x = 0.9, y = 0.25, label = "B", hjust = 1, vjust = 0, size = 6) +
  annotate("text", x = 0.8, y = 0, label = "1-(1/n)", hjust = 1, vjust = 0, size = 3) + 
  annotate("text", x = 1, y = 0, label = "1", hjust = 1, vjust = 0, size = 3) + 
  annotate("text", x = 0, y = 1, label = "1", hjust = 1, vjust = 0, size = 3)
图 2.17: 最大基尼系数

2.8 R

如果我们要计算 \(x_1,...,x_n\) 的样本均值和样本方差，我们可以在 R 中输入：

x <- c(x_1,...,x_n)

然后可以按下 回车 键，并输入 mean(x)，然后当再次按下 回车 时，就会打印出样本均值。当输入 var(x)，然后按下 回车 时，就会打印出样本方差。

例如，我们要计算：4，6，12，9，21，14 的样本均值和样本方差，我们可以用 R 按照如下方式获取：

> x <- c(4,6,12,9,21,14)
> mean(x)
[1] 11
> var(x)
[1] 37.6
>

在 R 中，我们不需要输入行首的 >，R 会自动为每一行补充提示符 >，以提示我们输入代码。R 中的命令还有：

  • sum(x)：返回向量 x 中所有数据的和。
  • median(x)：返回向量 x 中的所有数据的中位数。
  • sd(x)：返回向量 x 中所有数据的标准差。

如果 \(x=c(x_1,...,x_n)\)，\(y=c(y_1,...,y_n)\)，则：

  • \(x+y=c(x_1+y_1,...,x_n+y_n)\)
  • \(x-y=c(x_1-y_1,...,x_n-y_n)\)
  • \(x^2=c(x_1^2,...,x_n^2)\)
  • \(\frac{x}{y}=c(\frac{x_1}{y_1},...,\frac{x_n}{y_n}), \ \forall y_i \ne 0\)

如果 \(a\) 是一个实数，则：

  • \(a + x = c(a+x_1,...,a+x_n)\)
  • \(ax = c(ax_1,...,ax_n)\)
  • \(\frac{x}{a}=c(\frac{x_1}{a},...,\frac{x_n}{a}), \ \forall a \ne 0\)
  • \(\frac{a}{x}=c(\frac{a}{x_1},...,\frac{a}{x_n}), \ \forall x_i \ne 0\)

对于成对数据 \(\{(x_i,y_i)\}, \ i=1,...,n\)，我们可以使用如下命令获取其相关系数：

> x <- c(x1,...,xn) 
> y <- c(y1,...,yn) 
> cor(x, y)

可以用 plot(x, y) 获取 \(\{(x_i,y_i)\}\) 的散点图：

> plot(x, y)

假设我们有如下的数据：(4,10), (6,13), (12,22), (9,15), (21,30), (14,15)，我们可以使用如下的 R 命令获取其样本相关系数和散点图：

代码 
x <- c(4, 6, 12, 9, 21, 14)
y <- c(10, 13, 22, 15, 30, 15) 
cor(x, y)
[1] 0.9041494
代码 
plot(x, y)

R 还可以用于数学计算，例如我们可以使用如下命令计算 \(\frac{18\sqrt{177}}{677}\)：

代码 
18*sqrt(177)/677 
[1] 0.3537288
  • 如果我们从向量 x <- c(x1, ..., xn) 中选择第 \(i\) 个位置的元素 \(x_i\)，我们可以使用命令 x[i]。
  • 如果希望从向量 x 中选择 \(i\) 到 \(j\) 的元素并构成一个新的向量，我们可以使用命令 x[i:j]。
代码 
x <- c(3, 18, 9, 7, 22, 5, 17)
x[3]
[1] 9
代码 
x[2:5]
[1] 18  9  7 22

我们可以使用 R 计算基尼系数。对于向量 x <- c(x1, ..., xn)，先利用 sort(x) 得到 关于向量 x 中数据的递增排序。

代码 
x <- c(9, 7, 22, 5, 17)
sort(x)
[1]  5  7  9 17 22

如果我们想使用 式 2.19 和 式 2.20 来计算 习题 2.11 中的基尼系数，我们可以使用如下的代码：

代码 
y <- c(9, 7, 22, 5, 17)
x <- sort(y)
s <- cumsum(x)[1:4]
B <- 1/10 + sum(s) / (5 * sum(x))
G <- 1 - 2 * B
G 
[1] 0.2933333

我们也可以在 R 中定义 函数。例如，我们想在 R 中定义函数 \(f\) 以实现 \(f(x) = x^2\)，可以使用如下的代码：

代码 
f = function(x){x * x}
f(4)
[1] 16

如果我们想实现 \(f(x) = e^x\)，可以使用如下的代码：

代码 
f = function(x){exp(x)}
f(1)
[1] 2.718282

使用 函数 的概念，我们可以使用如下的代码以获取 9, 7, 22, 5, 17 的基尼系数：

代码 
s = function(x, j){sum(x[1:j])}

y <- c(9, 7, 22, 5, 17)
x <- sort(y)
B <- 1/10 + (s(x, 1) + s(x, 2) + s(x, 3) + s(x, 4)) / (5 * s(x, 5))
G <- 1 - 2 * B
G 
[1] 0.2933333

习题

  1. 以下是 1997 年 6 月，旧金山湾区每加仑标准无铅汽油的价格的样本数据。

    3.88, 3.90, 3.93, 3.90, 3.93, 3.96, 3.88, 3.94, 3.96, 3.88, 3.94, 3.99, 3.98

    使用如下的方式对数据进行描述：

    • 频率表
    • 相对频率的线图

  2. 解释一下如何构建一个饼图。如果数据集中的某个数据值的相对频率为 r，那么在饼图中，该扇区将会位于什么角度？

  3. 以下是西半球四个地区的石油储量估计数据（以十亿桶为单位）：

    • United States 38.7
    • South America 22.6
    • Canada 8.8
    • Mexico 60.0

    请使用饼图来描述如上的数据。

  4. 选择一本书或一篇文章，计算前 100 个句子中每个句子的单词数，并使用茎叶图来展现这些数据。现在选择另一本由不同作者撰写的书或文章，并做同样的操作。这两个茎叶图看起来是否相似？你认为这种方法能否有效地判断不同的文章是否由不同的作者所写？

  5. 以下是每日上班通勤时间（以分钟为单位）的频率表：

    Travel time Frequency
    15 6
    18 5
    22 4
    23 3
    24 4
    25 2
    26 4
    32 3
    36 1
    48 1
    • 频率表中的数据包含了几天的数据？
    • 频率表中的总的通勤时间是多少？

  6. 表 2.10 列出了 1985 年至 2006 年间，美国商业航空事故的次数以及由此造成的死亡总人数。

    • 绘制事故数的频率表
    • 绘制事故数的折线图
    • 绘制事故数的累积相对频率图
    • 计算事故数的样本均值
    • 计算事故数的样本中位数
    • 找出事故数的样本众数
    • 计算事故数的样本标准差
代码 
library(knitr)
df <- read.table("../misc/US_airline.csv", header=TRUE, sep=",")
kable(df, align="l")
表 2.10: 美国商业航空公司 1985~2006 安全数据，数据来源： National Transportation Safety Board。
Year Departures..millions. Accidents Fatalities
1985 6.1 4 197
1986 6.4 2 5
1987 6.6 4 231
1988 6.7 3 285
1989 6.6 11 278
1990 7.8 6 39
1991 7.5 4 62
1992 7.5 4 33
1993 7.7 1 1
1994 7.8 4 239
1995 8.1 2 166
1996 7.9 3 342
1997 9.9 3 3
1998 10.5 1 1
1999 10.9 2 12
2000 11.1 2 89
2001 10.6 6 531
2002 10.3 0 0
2003 10.2 2 22
2004 10.8 1 13
2005 10.9 3 22
2006 11.2 2 50

  1. 根据 表 2.10 的数据，

    • 绘制事故导致的死亡人数的直方图
    • 绘制事故导致的死亡人数的茎叶图
    • 计算事故导致的死亡人数的的样本均值
    • 计算事故导致的死亡人数的的样本中位数
    • 计算事故导致的死亡人数的的样本标准差

  2. A 镇成年女性的体重样本均值大于 B 镇成年女性的体重样本均值。此外，A 镇成年男性的体重样本均值也大于 B 镇成年男性的体重样本均值。我们是否可以得出结论：A 镇成年人的体重样本均值大于 B 镇成年人的体重样本均值？请给出你的解释。

  3. 首位数定律（Benford’s law for first digits）指出：在许多现实生活中的数值数据集中，首位数字并不是以相等的比例出现，而是偏向于较小的数字。更具体地说，首位数定律指出，首位非零数字为 \(i, i=1,...,9\) 的数据比例大约为 \(log_{10}{(i+1)}\)。例如，\(log_{10}{(2)}=0.301\)，这表明大约 30.1% 的数据的首位数字是 1。表 2.11 给出了首位数定律中以 1~9 作为首位数字的数据的比例。

    有趣的是，首位数定律已被证明可以适用于各种现实生活数据集，包括电费账单、街道地址、股票价格、人口数量、死亡率、河流长度、物理和数学常数，并且当数据值广泛分布时似乎最为准确。首位数定律最初是由美国天文学家西蒙・纽科姆（Simon Newcomb）在 1881 年发表的。1938 年，物理学家弗兰克・本福德（Frank Benford）在 20 个不同领域的数据集上测试了首位数定律，并表明在大多数情况下它都是一个很好的拟合。弗兰克・本福德测试的数据包括：河流的表面积、美国城市的人口规模、物理常数和分子量（molecular weights）等。

    物理考试中的一个多选题为：以下哪个是 20℃ 下，100% 过氧化氢溶液的密度（单位为 \(g/{cm}^3\)），(a) 7.3316、(b) 6.2421、(c) 1.4512、(d) 8.1818。如果你对过氧化氢一无所知，你会猜哪个答案是正确的？

    表 2.11: 首位数定律
    First digit Proportion of data having it as first digit
    1 0.301
    2 0.176
    3 0.125
    4 0.097
    5 0.079
    6 0.067
    7 0.058
    8 0.051
    9 0.046

  4. A 公司总共有 100 名员工，而 B 公司总共有 110 名员工。假设 A 公司的所有员工的薪水总和比 B 公司高。

    • 对于 A 公司工资的中位数与 B 公司工资的中位数来说意味着什么？
    • 对于 A 公司工资的平均数与 B 公司工资的平均数来说意味着什么？

  5. 一个包含 198 个数据的数据集中，前 99 个数据的样本均值等于 120，而后 99 个数据的样本均值等于 100。关于整个数据集的样本均值，你能得出什么结论？

    • 关于整个数据集的样本中位数，你能得出什么结论？
    • 关于整个数据集的样本众数，你能得出什么结论？

  6. 下表给出了 1922 年英格兰的重大道路交通事故中，按照年龄、性别汇总的行人死亡人数数据。

    • 估算男性年龄的样本均值
    • 估算女性年龄的样本均值
    • 估算死亡男性的四分位数
    • 估算死亡女性的四分位数
    Age Range Number of Males Number of Females
    0-5 120 67
    5-10 184 120
    10-15 44 22
    15-20 24 15
    20-30 23 25
    30-40 50 22
    40-50 60 40
    50-60 102 76
    60-70 167 104
    70-80 150 90
    80-100 49 27

  7. 以下是 12 个相邻位置发现的煤炭样本中的含灰量占比数据：

    9.2, 14.1, 9.8, 12.4, 16.0, 12.6, 22.7, 18.9, 21.0, 14.5, 20.4, 16.9
    • 计算如上数据的样本均值
    • 计算如上数据的样本标准差

  8. 5 个数据的样本均值和样本方差分别是 \(\overline{x} = 104\)，\(s^2 = 16\)，如果其中的 3 个数是 102， 100， 105，另外两个数是什么？

  9. 假设你得到了美国 50 个州中每个州所有工人的平均工资。

    • 你认为这 50 个州的平均工资的样本均值会等于整个美国的工人平均工资吗？
    • 如果对（a）的回答是否定的，请解释除了这 50 个平均值之外，还需要什么其他信息来确定整个国家的样本平均薪资。同时，解释你将如何使用这些额外信息来计算这个整个美国的工人平均工资。

  10. 如下是 40 个晶体管的使用寿命（单位为小时）：

    112, 121, 126, 108, 141, 104, 136, 134, 
121, 118, 143, 116, 108, 122, 127, 140,
113, 117, 126, 130, 134, 120, 131, 133,
118, 125, 151, 147, 137, 140, 132, 119,
110, 124, 132, 152, 135, 130, 136, 128
    • 计算如上数据的样本均值、样本中位数、样本众数
    • 给出如上数据的累积相对频率图

  11. 一个实验测量了 50 个粘土样本干燥后的收缩比率，并记录了以下数据：

    18.2 21.2 23.1 18.5 15.6 20.8 19.4 15.4 21.2 13.4 
16.4 18.7 18.2 19.6 14.3 16.6 24.0 17.6 17.8 20.2 
17.4 23.6 17.5 20.3 16.6 19.3 18.5 19.3 21.2 13.9 
20.5 19.0 17.6 22.3 18.4 21.2 20.4 21.4 20.3 20.1 
19.6 20.6 14.8 19.7 20.5 18.0 20.8 15.8 23.1 17.0
      1. 画出如上数据的茎叶图
      1. 计算样本均值，样本中位数，样本众数
      1. 计算样本方差
      1. 从 13% 开始，以 1% 为间隔将如上的数据进行分组，并绘制直方图
      1. 对于分组数据而言，假设每个数据点实际上位于其所在区间的中点，计算样本均值和样本方差，并与 (b) 和 (c) 的结果进行比较。为什么它们会不同？

  12. 如下是计算 \(\{x_i;\ i=1,...,n\}\) 样本均值和样本方差的一种快速算法。首先计算前 \(j(j \ge 2)\) 个数据的样本均值和方差：\(\overline{x}_j = \frac{\sum_{i=1}^{j}{x_i}}{j}\)，\(s_j^2=\frac{\sum_{i=1}^{j}{(x_i - \overline{x})^2}}{j - 1}\)。其中，\(\overline{x}_1 = x_1\)，\(s_1^2 = 0\)。则：

    \[ \begin{align} & \overline{x}_{j+1} = \overline{x}_j + \frac{x_{j + 1} - \overline{x}}{j + 1} \\ & s_{j+1}^2 = (1 - \frac{1}{j})s_j^2 + (j + 1)(\overline{x}_{j+1} - \overline{x}_j)^2 \end{align} \]

      1. 使用如上的算法计算 3, 4, 7, 2, 9, 6 的样本均值和样本方差
      1. 使用普通的计算方式来校验 (a) 的结果
      1. 验证如上算法中的 \(\overline{x}_{j+1}\) 和 \(\overline{x}_{j}\) 的关系

  13. 对于 表 2.5 的数据，

    • 计算 1 月的平均气温的 90-分位值
    • 计算 7 月的平均气温的 75-分位值

  14. 根据 《纽约时报》在 2013 年 8 月 1 日两周前发布的讣告，找出如下的死亡年龄的四分位数。

    92, 90, 92, 74, 69, 80, 94, 98, 65, 96, 
84, 69, 86, 91, 88, 74, 97, 85, 88, 68, 
77, 94, 88, 65, 76, 75, 60, 69, 97, 92, 
85, 70, 80, 93, 91, 68, 82, 78, 89

  15. 我们按照各个大学在谷歌上的月搜索量对大学进行月度排名，在截止到 2013 年 6 月的 114 个月中，获得过月榜单 TOP 10 的大学的上榜次数如下表所示。

    University Number of Months in Top 10
    Harvard University 114
    University of Texas, Austin 114
    University of Michigan 114
    Stanford University 113
    University of California Los Angeles (UCLA) 111
    University of California Berkeley 97
    Penn State University 94
    Massachusetts Institute of Technology (MIT) 66
    University of Southern California (USC) 63
    Ohio State University 52
    Yale University 48
    University of Washington 33
    • 计算样本均值
    • 计算样本方差
    • 计算样本的四分位数

  16. 填写缺失的单词或短语以完成以下句子：“如果向一组数字中增加一个新的数字，如果新的数字____，则该数组的样本均值将增加。”

  17. 用箱线图表示第 20 题中的数据。

  18. 在一个石油化工企业中，在 36 个随机选择的时间点测量了平均颗粒物浓度（单位为 \(mg/m^3\)），得到了如下的浓度数据：

    5, 18,  15, 7,  23, 220, 130, 85, 103, 25, 
80, 7,  24, 6,  13, 65,  37,  25, 24,  65, 
82, 95, 77, 15, 70, 110, 44,  28, 33,  81, 
29, 14, 45, 92, 17, 53
    • 使用直方图描述如上数据
    • 如上的直方图是近似正态直方图吗？

  19. 一位化学工程师想要研究盐水蒸发池中水的蒸发率，他获得了 4 年中盐水蒸发池 7 月份每天蒸发英寸数的 55 个数据。这些数据以下面的茎叶图给出，其最小数据为 0.02 英寸，最大数据为 0.56 英寸。

    .0    2,6
.1    1,4
.2    1,1,1,3,3,4,5,5,5,6,9
.3    0,0,2,2,2,3,3,3,3,4,4,5,5,5,6,6,7,8,9 
.4    0,1,2,2,2,3,4,4,4,5,5,5,7,8,8,8,9,9 
.5    2,5,6

    计算如上数据的：

    • 样本均值
    • 样本中位数
    • 样本标准差
    • 如上的数据看起来符合近似正态分布吗？
    • 位于样本均值 1 个标准差以内的数据占比是多少？

  20. 以下是加利福尼亚大学伯克利分校工业工程与运筹学系录取的最近 30 名学生的 GPA（grade point averages）。

    3.46, 3.72, 3.95, 3.55, 3.62, 3.80, 3.86, 3.71, 3.56, 3.49, 
3.96, 3.90, 3.70, 3.61, 3.72, 3.65, 3.48, 3.87, 3.82, 3.91, 
3.69, 3.67, 3.72, 3.66, 3.79, 3.75, 3.93, 3.74, 3.50, 3.83
      1. 使用茎叶图来绘制如上的数据
      1. 计算样本均值 \(\overline{x}\)
      1. 计算样本标准差 \(s\)
      1. 计算位于 \([\overline{x} - 1.5s, \overline{x} + 1.5s]\) 的数据比例，并和切比雪夫不等式给出的下限比例进行对比
      1. 计算位于 \([\overline{x} - 2s, \overline{x} + 2s]\) 的数据比例，并和切比雪夫不等式给出的下限比例进行对比

  21. 第 26 题中的数据是否近似于正态分布？对于第 26 题的 (d) 和 (e) ，请使用经验法则给出的近似占比与实际占比进行比较。

  22. 你是否期望一个健身俱乐部所有会员的体重直方图会近似于正态分布？

  23. 对于第 16 题中的数据：

      1. 计算样本均值和样本中位数
      1. 这些数据符合近似正态分布吗？
      1. 计算样本标准差 \(s\)
      1. 计算位于 \([\overline{x} - 1.5s, \overline{x} + 1.5]\) 区间的数据的占比
      1. 对比 (d) 的结果和经验法则给出的结果
      1. 对比 (d) 的结果和切比雪夫不等式给出的结果

  24. 以下是 12 名考试成绩大致相同的法学院学生的身高和起薪数据：

    Height (inches) Salary
    64 91
    65 94
    66 88
    67 103
    69 77
    70 96
    72 105
    72 88
    74 122
    74 102
    75 90
    76 114
    • 绘制数据的散点图
    • 计算样本相关系数

  25. 根据人们的站立姿势数据形成一个随机样本。针对样本中的人，还额外记录了每个人在过去一年中经历背痛的天数。令研究人员惊讶的是，这些数据表明良好的站姿与背痛天数之间存在正相关关系。这是否意味着良好的站姿会导致背痛？

  26. 如果我们把美国 50 个州中每个州的居民平均收入和居住在该州的外国出生的移民数量绘制为成对的数据图，那么这些成对数据将呈现正相关关系。我们能否得出结论，移民居民往往比土生土长的美国人的收入更高？如果不是，还能如何解释这个现象？

  27. 随机抽取 12 名高中三年级学生，并要求他们估算自己每周平均学习的小时数。以下的数据给出了估算结果和学生的 GPA。

    Hours GPA
    6 2.8
    14 3.2
    3 3.1
    22 3.6
    9 3.0
    11 3.3
    12 3.4
    5 2.7
    18 3.1
    24 3.8
    15 3.0
    17 3.9
    • 计算学习小时数和 GPA 之间的样本相关系数。

  28. 验证样本相关系数的特性 3。

  29. 验证样本相关系数的特性 4。

  30. 在一项针对二至四年级的儿童研究中，研究人员针对每个学生都做了阅读测试。在研究结果数据时，研究人员注意到学生的阅读测试分数和身高之间存在正相关关系。研究人员得出结论，身高较高的孩子阅读能力更好，因为他们可以更容易地看到黑板。你认为呢？

  31. 最近的一项研究发现，母乳喂养的婴儿与 6 岁时进行的词汇测试分数之间存在正相关关系。讨论在解释这项研究结果时可能遇到的困难。

  32. 一个群体的收入分别为 25、32、60、40、38、50，绘制其洛伦茨曲线并计算其基尼系数。

  33. 根据如下的年收入（以千美元为单位）频率表，绘制该群体的洛伦茨曲线并计算其基尼系数：

    Value Frequency
    30 2
    50 4
    60 5
    90 4
    100 3
    120 2

  34. 如果样本中所有的数据都乘以一个大于 0 的常数 \(c\)，基尼系数将会发生什么变化？如果样本中所有的数据都加上一个大于 0 的常数 \(c\)，基尼系数又会发生什么变化？下降，保持不变，还是不好说？