\documentclass[a4paper, 10.5pt, twoside, openany]{book}
\usepackage{amsfonts}
\usepackage{array}
\usepackage{boxedminipage, fancybox}
\usepackage{caption}
\usepackage{color}
\usepackage[colorlinks,linkcolor=blue]{hyperref}
\usepackage{ctex}
\usepackage{datetime}
\usepackage[dvipsnames]{xcolor}
\usepackage{enumerate}
\usepackage{epsfig,graphicx,subfigure}
\usepackage{extarrows}
\usepackage{fancyheadings}
\usepackage{float}
\usepackage{geometry}
\usepackage{listings}
\usepackage{longtable}
\usepackage{makeidx}
\usepackage{mathrsfs}
\usepackage{multirow}
\usepackage{natbib}
\usepackage{pifont}
\usepackage{rotating}
\usepackage{setspace}
\usepackage{shadow}
\usepackage{stmaryrd, amssymb, amsmath}
\usepackage{tabularx}
\usepackage{url}
\usepackage{varioref}
\usepackage{verbatim}
\usepackage{wrapfig}
\usepackage{xcolor}
\geometry{left=2.0cm, right=2.0cm, top=2.5cm, bottom=2.5cm}
\linespread{1.5}
\definecolor{mygray}{rgb}{0.85, 0.85, 0.85}
\newcommand{\codeinline}[1]{\colorbox{mygray}{\lstinline|#1|}}
%% ----------------------------------------------------------------------------------------------------------------------------------------------------------------------
\title{\Huge \bf 《多元统计分析》课后作业}
\author{\kaishu 姓名：\underline{\quad 李倩倩 \quad} \\[5mm]
\kaishu 学号：\underline{\quad 2024017349 \quad} \\[5mm]
\kaishu 班级：\underline{\quad 统计 24-1 班 \quad} \\[50mm]
\kaishu 中国石油大学（北京）克拉玛依校区文理学院数学与统计系
}
\date{\today}
\begin{document}
% -------------------------------------------- 封面页 --------------------------------------------
\frontmatter
\maketitle
% -------------------------------------------- 作业要求 --------------------------------------------
\chapter{作业要求}
\begin{enumerate}
\item 可以和其他同学讨论作业当中的问题，但应当自己独立完成作业
\item 计算、证明等要有过程，要有主要步骤的说明
\item 请将计算、绘图所用的 R 代码以及生成的结果和图像一并添加在作业文件当中
\item 请使用 \LaTeX 编辑并生成 PDF 格式的文件，第 X 周作业文件命名方式：学号-姓名-X.pdf
\item 评分标准：每一问得分 $\in \left\{ 2 ,\, 1 ,\, 0 \right\}$
\begin{itemize}
\item 2:~ 按时完成并上交作业，且答案基本正确
\item 1:~ 按时完成并上交作业，且答案部分正确
\item 0:~ 答案完全错误，或者迟交作业(规定时间72小时之后)
\end{itemize}
\item 请将完成的 PDF 格式的作业文件发送至邮箱：xiaolei@cup.edu.cn
\item 每位同学可以有一次迟交作业的机会，但不得晚于规定时间三日之后
\item 第 8 周作业截止时间：2026年5月8日24:00
\end{enumerate}
\tableofcontents
% -------------------------------------------- 正文部分 --------------------------------------------
\mainmatter
% -------------------------------------------- 第 8 周作业 --------------------------------------------
\chapter{第 8 周作业}
{\kaishu \color{blue} 第 8 周作业截止时间：} 2026年5月8日24:00
{\kaishu \color{blue} 第 8 周作业完成时间：} \today \space \currenttime % 请勿编辑、删除本行！
\lstset{
language=R,
basicstyle=\small\ttfamily,
backgroundcolor=\color{mygray},
breaklines=true,
showstringspaces=false,
frame=single,
rulecolor=\color{black!30},
keywordstyle=\color{blue!80!black},
commentstyle=\color{green!55!black},
stringstyle=\color{red!70!black}
}
\begin{enumerate}
\item 利用数据矩阵的因子分解方法，简要分析瑞士银行钞票数据集 (mclust 包中的 banknote 数据集).
\begin{enumerate}
\item {\color{TealBlue} [2 分]} 利用 R 中的 scale() 函数对数据进行标准化，将标准化之后的数据集记为 $\mathcal{X}$.
{\color{red} \heiti 【解】}
加载 mclust 包并读取 banknote 数据集，提取 6 个数值变量（Length, Left, Right, Bottom, Top, Diagonal），利用 \codeinline{scale()} 对每列进行标准化（各列均值减为 0，样本标准差化为 1），所得矩阵记为 $\mathcal{X}$，维度为 $200 \times 6$，即 $n = 200$，$p = 6$.
\begin{lstlisting}[language=R]
rm(list = ls(all = TRUE))
library(mclust)
data(banknote)
n <- nrow(banknote) # n = 200
p <- ncol(banknote) - 1 # p = 6
# Extract 6 numeric variables (drop the first column: Status)
X_data <- banknote[, -1]
# Standardize: subtract column mean, divide by sample std (ddof = n-1)
X <- scale(X_data)
# Verify: column means = 0, column sds = 1
round(colMeans(X), 6)
round(apply(X, 2, sd), 6)
\end{lstlisting}
标准化后，$\mathcal{X}$ 各列均值均为 $0$，各列样本标准差均为 $1$，即
\[
\mathcal{X}_{ij} = \frac{x_{ij} - \bar{x}_j}{s_j}, \quad i=1,\ldots,200,\; j=1,\ldots,6.
\]
其中 $s_j = \sqrt{\dfrac{1}{n-1}\sum_{i=1}^{n}(x_{ij}-\bar{x}_j)^2}$ 为第 $j$ 个变量的样本标准差.
\item {\color{TealBlue} [2 分]} 为处理更为规范，我们将 $\mathcal{X}$ 的所有元素除以 $\sqrt{n-1}$，其中 $n$ 为样本容量，得到的数据矩阵记为 $\mathcal{Y}$，
现在 $\mathcal{Y}$ 即是我们要分析的数据矩阵.
{\color{red} \heiti 【解】}
将标准化数据矩阵 $\mathcal{X}$ 的每个元素除以 $\sqrt{n-1} = \sqrt{199} \approx 14.1067$，得到分析数据矩阵 $\mathcal{Y}$：
\begin{lstlisting}[language=R]
Y <- X / sqrt(n - 1) # sqrt(199) = 14.10674 (the analysis data matrix)
\end{lstlisting}
由此，$\mathcal{Y}_{ij} = \dfrac{x_{ij}-\bar{x}_j}{s_j\,\sqrt{n-1}}$，维度仍为 $200\times6$.
经此变换后，$\mathcal{R} = \mathcal{Y}^{\rm T}\mathcal{Y}$ 恰好等于原始数据的样本相关矩阵，详见下一小题.
\item {\color{TealBlue} [2 分]} 求矩阵 $\mathcal{R} = \mathcal{Y}^{\rm T} \mathcal{Y}$ 的特征值及其对应的单位特征向量. 矩阵 $\mathcal{R} = \mathcal{Y}^{\rm T} \mathcal{Y}$
是原始数据的相关矩阵.
{\color{red} \heiti 【解】}
计算 $\mathcal{R} = \mathcal{Y}^{\rm T}\mathcal{Y}$ 并对其进行谱分解，求得特征值与对应单位特征向量.
\begin{lstlisting}[language=R]
# Compute correlation matrix R = Y^T Y
R <- t(Y) %*% Y
# Verify: R equals the sample correlation matrix of the raw data
round(R, 4)
round(cor(X_data), 4) # should be identical to R above
# Spectral decomposition: eigen() returns eigenvalues in descending order
eig <- eigen(R)
lambda <- eig$values # eigenvalue vector
U <- eig$vectors # unit eigenvector matrix (columns are eigenvectors)
round(lambda, 4)
round(U, 4)
\end{lstlisting}
\textbf{相关矩阵 $\mathcal{R}$：}
\[
\mathcal{R} = \begin{pmatrix}
1.0000 & 0.2313 & 0.1518 & -0.1898 & -0.0613 & 0.1943 \\
0.2313 & 1.0000 & 0.7433 & 0.4138 & 0.3623 & -0.5032 \\
0.1518 & 0.7433 & 1.0000 & 0.4868 & 0.4007 & -0.5165 \\
-0.1898 & 0.4138 & 0.4868 & 1.0000 & 0.1419 & -0.6230 \\
-0.0613 & 0.3623 & 0.4007 & 0.1419 & 1.0000 & -0.5940 \\
0.1943 & -0.5032 & -0.5165 & -0.6230 & -0.5940 & 1.0000
\end{pmatrix}
\]
（行列对应顺序：Length, Left, Right, Bottom, Top, Diagonal）
\textbf{特征值（降序排列）：}
\[
\lambda_1 = 2.9456,\quad \lambda_2 = 1.2781,\quad \lambda_3 = 0.8690,\quad
\lambda_4 = 0.4498,\quad \lambda_5 = 0.2687,\quad \lambda_6 = 0.1889
\]
验证：$\sum_{k=1}^{6}\lambda_k = 6.0000 = p$，与相关矩阵迹 $\mathrm{tr}(\mathcal{R})=6$ 一致. \medskip
\textbf{对应的单位特征向量（列向量 $\boldsymbol{u}_k$）：}
\begin{center}
\small
\begin{tabular}{lcccccc}
\hline
变量 & $\boldsymbol{u}_1$ & $\boldsymbol{u}_2$ & $\boldsymbol{u}_3$ & $\boldsymbol{u}_4$ & $\boldsymbol{u}_5$ & $\boldsymbol{u}_6$ \\
\hline
Length & $-0.0070$ & $-0.8155$ & $-0.0177$ & $-0.5746$ & $-0.0588$ & $-0.0311$ \\
Left & $\phantom{-}0.4678$ & $-0.3420$ & $\phantom{-}0.1034$ & $\phantom{-}0.3949$ & $\phantom{-}0.6395$ & $\phantom{-}0.2977$ \\
Right & $\phantom{-}0.4867$ & $-0.2525$ & $\phantom{-}0.1235$ & $\phantom{-}0.4303$ & $-0.6141$ & $-0.3492$ \\
Bottom & $\phantom{-}0.4068$ & $\phantom{-}0.2662$ & $\phantom{-}0.5835$ & $-0.4037$ & $-0.2155$ & $\phantom{-}0.4624$ \\
Top & $\phantom{-}0.3679$ & $\phantom{-}0.0915$ & $-0.7876$ & $-0.1102$ & $-0.2198$ & $\phantom{-}0.4190$ \\
Diagonal & $-0.4935$ & $-0.2739$ & $\phantom{-}0.1139$ & $\phantom{-}0.3919$ & $-0.3402$ & $\phantom{-}0.6318$ \\
\hline
\end{tabular}
\end{center}
各 $\boldsymbol{u}_k$ 均满足 $\|\boldsymbol{u}_k\|=1$ 且两两正交.
\item {\color{TealBlue} [2 分]} 因子变量是数据矩阵 $\mathcal{Y}$ 的六个变量 $Y_1 \sim Y_6$ 的线性组合，写出前两个因子变量 $Z_1$ 和 $Z_2$ 的表达式.
{\color{red} \heiti 【解】}
由讲义第 10 章，第 $k$ 个因子变量 $Z_k = \mathcal{Y}\boldsymbol{u}_k$，其元素为初始变量 $Y_1,\ldots,Y_6$ 的线性组合，线性组合系数由 $\boldsymbol{u}_k$ 的各分量给出（$Y_1=$Length，$Y_2=$Left，$Y_3=$Right，$Y_4=$Bottom，$Y_5=$Top，$Y_6=$Diagonal）.
\begin{lstlisting}[language=R]
# Direction vectors of the first two factor axes
u1 <- U[, 1] # corresponds to the largest eigenvalue lambda_1
u2 <- U[, 2] # corresponds to the second largest eigenvalue lambda_2
\end{lstlisting}
将 $\boldsymbol{u}_1$ 和 $\boldsymbol{u}_2$ 的分量代入，得前两个因子变量的表达式：
\begin{align}
Z_1 &= -0.0070\,Y_1 + 0.4678\,Y_2 + 0.4867\,Y_3 + 0.4068\,Y_4 + 0.3679\,Y_5 - 0.4935\,Y_6 \notag\\
Z_2 &= -0.8155\,Y_1 - 0.3420\,Y_2 - 0.2525\,Y_3 + 0.2662\,Y_4 + 0.0915\,Y_5 - 0.2739\,Y_6 \notag
\end{align}
\textbf{说明：} $Z_1$ 的系数显示 Left、Right、Bottom、Top 对其有正向贡献，Diagonal 有较强负向贡献，Length 贡献几乎为零；$Z_2$ 主要由 Length（系数最大）主导.
\item {\color{TealBlue} [2 分]} 计算前两个因子变量 $Z_1$ 和 $Z_2$ 对应的特征值之和占所有特征值之和的比例.
{\color{red} \heiti 【解】}
由讲义定义，该比例 $\tau_2$ 反映了二维因子表示所能解释的总惯量（方差）的比例：
\begin{lstlisting}[language=R]
tau2 <- sum(lambda[1:2]) / sum(lambda)
cat("tau_2 =", round(tau2, 4), "\n")
# output: tau_2 = 0.7039
\end{lstlisting}
\[
\tau_2 = \frac{\lambda_1 + \lambda_2}{\displaystyle\sum_{k=1}^{6}\lambda_k}
= \frac{2.9456 + 1.2781}{6.0000}
= \frac{4.2237}{6.0000}
\approx \boxed{0.7039}
\]
前两个因子变量共解释了总方差的约 $\mathbf{70.39\%}$，二维图形可以较好地展现数据的主要结构.
\item {\color{TealBlue} [2 分]} 计算观测数据 (行) 在前两个因子变量 $Z_1$ 和 $Z_2$ 上的坐标值.
{\color{red} \heiti 【解】}
$n$ 个观测点在第 $k$ 个因子轴 $\boldsymbol{u}_k$ 上的坐标为 $\boldsymbol{z}_k = \mathcal{Y}\boldsymbol{u}_k$，全部观测的坐标矩阵为 $Z = \mathcal{Y} U_{[:,1:2]}$，维度为 $200\times 2$.
\begin{lstlisting}[language=R]
# Compute coordinates of 200 observations along Z1 and Z2
Z <- Y %*% U[, 1:2]
colnames(Z) <- c("Z1", "Z2")
# Preview first 6 rows (genuine) and rows 101-106 (counterfeit)
round(Z[1:6, ], 5)
round(Z[101:106, ], 5)
\end{lstlisting}
前 6 个观测（真钞）与第 101--106 个观测（伪钞）的坐标如下：
\begin{center}
\small
\begin{tabular}{cccc}
\hline
观测序号 & $Z_1$ & $Z_2$ & 类别 \\
\hline
1 & $\phantom{-}0.12356$ & $-0.11673$ & genuine \\
2 & $-0.16082$ & $\phantom{-}0.03810$ & genuine \\
3 & $-0.16104$ & $\phantom{-}0.00761$ & genuine \\
4 & $-0.16147$ & $\phantom{-}0.00620$ & genuine \\
5 & $-0.18612$ & $-0.00277$ & genuine \\
6 & $\phantom{-}0.05363$ & $-0.21841$ & genuine \\
\hline
101 & $\phantom{-}0.07817$ & $-0.05174$ & counterfeit \\
102 & $\phantom{-}0.10553$ & $-0.01264$ & counterfeit \\
103 & $\phantom{-}0.09100$ & $\phantom{-}0.02020$ & counterfeit \\
104 & $\phantom{-}0.11840$ & $-0.03290$ & counterfeit \\
105 & $\phantom{-}0.13066$ & $\phantom{-}0.01028$ & counterfeit \\
106 & $\phantom{-}0.08716$ & $-0.07267$ & counterfeit \\
\hline
\end{tabular}
\end{center}
真钞的 $Z_1$ 均值约为 $-0.1059$，伪钞的 $Z_1$ 均值约为 $+0.1059$，两类在 $Z_1$ 方向上有明显分离.
\item {\color{TealBlue} [2 分]} 作观测数据在前两个因子变量 $Z_1$ 和 $Z_2$ 上的散点图，将真钞与伪钞的数据点分别用不同的颜色表示，你能看到什么现象.
{\color{red} \heiti 【解】}
\begin{lstlisting}[language=R]
# Get banknote class labels
status <- banknote$Status
# Scatter plot on Z1-Z2: genuine = blue circle, counterfeit = red triangle
plot(Z[, 1], Z[, 2],
col = ifelse(status == "genuine", "steelblue", "firebrick"),
pch = ifelse(status == "genuine", 16, 17),
xlab = expression(Z[1]),
ylab = expression(Z[2]),
main = "Banknote: Observation Scatter Plot on Z1-Z2")
abline(h = 0, v = 0, lty = 2, col = "gray50")
legend("topright",
legend = c("Genuine (genuine)", "Counterfeit (counterfeit)"),
col = c("steelblue", "firebrick"),
pch = c(16, 17),
bty = "n")
\end{lstlisting}
\begin{figure}[H]
\centering
\includegraphics[width=0.75\textwidth]{C:/Users/35297/Documents/Rplot.png}
\caption{200 个观测在因子变量 $Z_1$-$Z_2$ 上的散点图（蓝色圆点为真钞，红色三角为伪钞）}
\label{fig:kde}
\end{figure}
\textbf{现象分析：}
\begin{itemize}
\item 沿 $Z_1$ 轴（第一因子变量），真钞（蓝色）集中于负半轴，伪钞（红色）集中于正半轴，两组有明显的分离.
\item 沿 $Z_2$ 轴（第二因子变量），两组的分离程度远不如 $Z_1$ 轴明显，两组在 $Z_2$ 方向上有较大重叠.
\item 这说明第一因子变量 $Z_1$ 是区分真伪钞票的主要判别方向，该因子主要反映 Left、Right、Bottom、Top（对角线两侧尺寸）与 Diagonal（对角线长度）之间的相对大小关系.
\end{itemize}
\item {\color{TealBlue} [2 分]} 计算变量数据 (列) 在前两个因子变量 $W_1$ 和 $W_2$ 上的坐标值.
{\color{red} \heiti 【解】}
由讲义定理 10.4，$p$ 个变量点在第 $k$ 个因子轴（$G_k$ 方向，对应 $\mathcal{Y}\mathcal{Y}^{\rm T}$ 的特征向量 $\boldsymbol{v}_k$）上的坐标为
\[
\boldsymbol{w}_k = \mathcal{Y}^{\rm T}\boldsymbol{v}_k = \sqrt{\lambda_k}\,\boldsymbol{u}_k,
\]
即 $W_k$ 是对应特征值的平方根与单位特征向量的乘积，不需要显式计算 $\boldsymbol{v}_k$.
\begin{lstlisting}[language=R]
# Variable coordinates: W_k = sqrt(lambda_k) * u_k
# Matrix form: W = U[,1:2] %*% diag(sqrt(lambda[1:2]))
W <- U[, 1:2] %*% diag(sqrt(lambda[1:2]))
rownames(W) <- colnames(X_data)
colnames(W) <- c("W1", "W2")
round(W, 4)
\end{lstlisting}
6 个变量在 $W_1$-$W_2$ 平面上的坐标为：
\begin{center}
\begin{tabular}{lcc}
\hline
变量 & $W_1 = \sqrt{\lambda_1}\,u_{j1}$ & $W_2 = \sqrt{\lambda_2}\,u_{j2}$ \\
\hline
Length & $-0.0120$ & $-0.9219$ \\
Left & $\phantom{-}0.8028$ & $-0.3866$ \\
Right & $\phantom{-}0.8353$ & $-0.2854$ \\
Bottom & $\phantom{-}0.6981$ & $\phantom{-}0.3010$ \\
Top & $\phantom{-}0.6314$ & $\phantom{-}0.1034$ \\
Diagonal & $-0.8469$ & $-0.3097$ \\
\hline
\end{tabular}
\end{center}
注意 $\boldsymbol{w}_k$ 即是变量 $Y_j$ 与第 $k$ 个因子变量 $Z_k$ 的协方差（相关系数），因此 $|W_{jk}|$ 越大表明该变量与第 $k$ 个因子的关联越密切.
\item {\color{TealBlue} [2 分]} 作变量点在前两个因子变量 $W_1$ 和 $W_2$ 上的散点图，你能看到什么现象.
{\color{red} \heiti 【解】}
\begin{lstlisting}[language=R]
# Scatter plot of variable points on W1-W2
plot(W[, 1], W[, 2],
type = "n",
xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1),
xlab = expression(W[1]),
ylab = expression(W[2]),
main = "Banknote: Variable Scatter Plot on W1-W2",
asp = 1)
text(W[, 1], W[, 2], rownames(W), cex = 1.0, col = "darkgreen")
abline(h = 0, v = 0, lty = 2, col = "gray50")
# Draw unit circle as reference
theta <- seq(0, 2 * pi, length.out = 200)
lines(cos(theta), sin(theta), col = "gray70", lty = 3)
\end{lstlisting}
\begin{figure}[H]
\centering
\includegraphics[width=0.65\textwidth]{C:/Users/35297/Documents/Rplot03.png}
\caption{6 个变量在因子变量 $W_1$-$W_2$ 上的散点图}
\label{fig:kde}
\end{figure}
\textbf{现象分析：}
\begin{itemize}
\item 在 $W_1$ 方向（第一因子轴）上，Left、Right、Bottom、Top 的坐标均为正值（约在 $0.63 \sim 0.84$），而 Diagonal 的坐标为负值（$-0.85$），Length 接近零. 这说明第一因子 $Z_1$ 主要刻画了钞票宽度/高度尺寸（Left、Right、Bottom、Top）与对角线长度（Diagonal）之间的对立关系.
\item 在 $W_2$ 方向（第二因子轴）上，Length 的坐标绝对值最大（$-0.9219$），独立于其他变量，说明第二因子 $Z_2$ 主要反映钞票长度（Length）的信息.
\item Left 与 Right 靠近，Bottom 与 Top 靠近，说明这两对变量高度相关（与相关矩阵中 Left-Right 相关系数 $0.74$，Bottom-Top 为 $0.14$ 相呼应）.
\item 从变量散点图可知：真钞与伪钞在钞票尺寸上的系统差异主要体现在 Diagonal 与 Left/Right/Bottom/Top 的比例关系上，这与第一因子变量 $Z_1$ 能够有效区分真伪钞票的结论一致.
\end{itemize}
\end{enumerate}
\end{enumerate}
\end{document}