本文与数理统计-代码合集的整理思路类似，以R代码为主，python对应代码辅助。

矩阵运算

    A <- matrix(c(3, 1, -1,
        5, 2, 4), nrow = 2, byrow = TRUE)
t(A) # 转置
B <- matrix(c(1, 2, -2,
            5, -3, 1), nrow = 2, byrow = TRUE)
4 * A # 矩阵标量乘法
A + B # 矩阵加法
A %*% t(B) # 矩阵乘法 

C <- c(-2,7,9)

A %*% C # 矩阵向量乘法（自动确定行/列向量）

D <- matrix(c(3, 4,
            2, 1), nrow = 2, byrow = TRUE)
solve(D) # 矩阵求逆

F <- matrix(c(3/2, -1/2, 0,
            -1/2, 3/2, 0,
            0, 0, 1), nrow = 3, byrow = TRUE)

F_eigen <- eigen(F)
F_value <- F_eigen$values # 矩阵特征值（向量形式）
F_vector <- F_eigen$vectors # 矩阵特征向量（矩阵形式）

library(psych)
tr(F) # 矩阵的迹
det(F) # 矩阵行列式
    import numpy as np
A = np.array([[3, 1, -1],
            [5, 2, 4]])
print(A.T)

B = np.array([[1, 2, -2],
            [5, -3, 1]])
print(4 * A)
print(A + B)
print(A @ B.T)

C = np.array([-2, 7, 9])
print(A @ C)

D = np.array([[3, 4],
            [2, 1]])
print(np.linalg.inv(D))

F = np.array([[3/2, -1/2, 0],
            [-1/2, 3/2, 0],
            [0, 0, 1]])
F_value, F_vector = np.linalg.eig(F)
print(F_value)
print(F_vector)
print(np.linalg.det(F))
print(np.trace(F))

    绘制椭圆

    • 旋转椭圆的二次型形式为

      [xy]A1[xy]=[xy][a1a32a32a2][xy]=a1x2+a2y2+a3xy=1\begin{bmatrix} x \\ y \end{bmatrix} A_1 \begin{bmatrix} x & y \end{bmatrix}=\begin{bmatrix} x \\ y \end{bmatrix} \begin{bmatrix} a_1 & \frac{a_3}{2} \\ \frac{a_3}{2} & a_2 \end{bmatrix} \begin{bmatrix} x & y \end{bmatrix} = a_1x^2 + a_2y^2 + a_3xy = 1 
        library(plotrix)
par(pty = "s")
plot(c(-2, 2), c(-2, 2), type = "n", xlab = "", ylab = "")
grid()
# 原始椭圆：(x,y)确定中心坐标，a和b确定轴长
draw.ellipse(x = 0, y = 0, a = 2, b = 1, angle = 0, deg = FALSE)

# 旋转椭圆的画法
A1 <- matrix(c(3, -sqrt(2),
            -sqrt(2), 2), nrow=2)
par(pty = "s")
plot(c(-1, 1), c(-1, 1),
    type = "n", xlab = "", ylab = "")
grid()
draw.ellipse(x = 0, y = 0,
            a = 1/sqrt(eigen(A1)$values[1]), # 对应特征值
            b = 1/sqrt(eigen(A1)$values[2]),
            angle = atan(eigen(A1)$vectors[2,1] / eigen(A1)$vectors[1,1]), # 对应最大特征值对应特征向量与x正半轴夹角
            deg = FALSE) 
# 原始坐标轴
arrows(0, 0, 1, 0)
arrows(0, 0, 0, 1)
text(1, -0.1, label = "x1")
text(-0.1, 1, label = "y1")
# 旋转后坐标轴
arrows(0, 0, -eigen(A1)$vectors[1,1], -eigen(A1)$vector[2,1], lty = 2)
arrows(0, 0, -eigen(A1)$vectors[1,2], -eigen(A1)$vector[2,2], lty = 2)
text(-eigen(A1)$vectors[1,1], -eigen(A1)$vector[2,1], label = "x2")
text(-eigen(A1)$vectors[1,2], -eigen(A1)$vector[2,2], label = "y2")

# 使用contour函数绘图
d <- function(x, y){
    3 * x^2 + 2 * y^2 - 2 * sqrt(2) * x * y
}
x1 <- y1 <- seq(-1.5, 1.5, length = 100)
z <- outer(x1, y1, d)
contour(x = x1, y = y1, z = z, levels = c(1, 2),
        labels=c('c^2=1', 'c^2=4'),
        main = expression(c^2 == 3*x^2 + 2*y^2 - 2*sqrt(2)*x*y)) 
grid()
legend('bottomright', lty=c(1, 2),
    legend=c(expression(c^2==1),
                expression(c^2==4)))
        import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

# 绘制椭圆
def draw_ellipse(ax, x, y, a, b, angle_rad, **kwargs):
    ellipse = Ellipse(xy=(x, y), width=2*a, height=2*b, 
                    angle=np.degrees(angle_rad), **kwargs)
    ax.add_patch(ellipse)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

ax1.set_aspect('equal') 
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.set_xticks(np.arange(-2, 2.5, 0.5))
ax1.set_yticks(np.arange(-2, 2.5, 0.5))
ax1.grid(True)
ax1.set_xticklabels([])
ax1.set_yticklabels([])

# 画原始椭圆
draw_ellipse(ax1, x=0, y=0, a=2, b=1, angle_rad=0, edgecolor='blue', facecolor='none')
ax1.set_title("Original Ellipse")
ax2.set_aspect('equal')
ax2.set_xlim(-1, 1)
ax2.set_ylim(-1, 1)
ax2.grid(True)
ax2.set_xticklabels([])
ax2.set_yticklabels([])

A1 = np.array([[3, -np.sqrt(2)],
            [-np.sqrt(2), 2]])

eigenvalues, eigenvectors = np.linalg.eig(A1)
idx = np.argsort(eigenvalues)[::-1] 
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

a = 1 / np.sqrt(eigenvalues[0]) 
b = 1 / np.sqrt(eigenvalues[1])

angle_rad = np.arctan(eigenvectors[1, 0] / eigenvectors[0, 0])

# 画旋转椭圆
draw_ellipse(ax2, x=0, y=0, a=a, b=b, angle_rad=angle_rad, edgecolor='blue', facecolor='none')

ax2.arrow(0, 0, 1, 0, head_width=0.03, head_length=0.05, fc='k', ec='k')
ax2.arrow(0, 0, 0, 1, head_width=0.03, head_length=0.05, fc='k', ec='k')
ax2.text(1.05, -0.05, "x1", fontsize=12)
ax2.text(-0.1, 1.05, "y1", fontsize=12)

v1 = -eigenvectors[:, 0]
v2 = -eigenvectors[:, 1]

ax2.arrow(0, 0, v1[0], v1[1], head_width=0.03, head_length=0.05, fc='r', ec='r', linestyle='--')
ax2.arrow(0, 0, v2[0], v2[1], head_width=0.03, head_length=0.05, fc='g', ec='g', linestyle='--')

ax2.text(v1[0] * 1.1, v1[1] * 1.1, "x2", fontsize=12, color='r')
ax2.text(v2[0] * 1.1, v2[1] * 1.1, "y2", fontsize=12, color='g')

ax2.set_title("Rotated Ellipse (Quadratic Form)")

plt.tight_layout()
plt.show()


from matplotlib.lines import Line2D

def d(x, y):
    return 3 * x**2 + 2 * y**2 - 2 * np.sqrt(2) * x * y

x1 = np.linspace(-1.5, 1.5, 100)
y1 = np.linspace(-1.5, 1.5, 100)
X, Y = np.meshgrid(x1, y1)
Z = d(X, Y)

fig, ax = plt.subplots(figsize=(8, 6))

contour = ax.contour(X, Y, Z, levels=[1, 2], colors='tab:blue')

fmt = {1: 'c^2=1', 2: 'c^2=4'}
ax.clabel(contour, inline=True, fontsize=10, fmt=fmt)

ax.set_title(r'$c^2 = 3x^2 + 2y^2 - 2\sqrt{2}xy$', fontsize=14)

ax.grid(True, linestyle='-', alpha=0.7)

custom_lines = [
    Line2D([0], [0], color='tab:blue', lw=2, linestyle='-'),
    Line2D([0], [0], color='tab:blue', lw=2, linestyle='--')
]
ax.legend(custom_lines, [r'$c^2=1$', r'$c^2=4$'], loc='lower right')

ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)

plt.show()

      正态分布

        # 设置随机数种子
set.seed(seed = 0)
# 产生来自 N(1, 2^2) 的 40 个样本
x1 <- rnorm(n = 40, mean = 1, sd = 2)
# 使用直方图查看样本分布
hist(x1)
        import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
x1 = np.random.normal(loc=1, scale=2, size=40)
plt.hist(x1)
plt.show()

        • 考察不同协方差矩阵对正态分布的影响：

            library(MASS)
# 均值向量
mu <- c(2,1)
# 协方差矩阵
Sigma1 <- matrix(c(5, 4,
                4, 5), ncol = 2, byrow = TRUE)
Sigma2 <- diag(c(3,3)) # diag里第一个参数是对角化向量，后面为矩阵大小（默认与向量长度一致）
Sigma3 <- matrix(c(5, -4,
                -4, 5), ncol = 2, byrow = TRUE)
# 设置随机数种子
set.seed(seed = 0)
# 生成 3 组样本（均值向量相同，协方差矩阵各异）
X1 <- mvrnorm(n = 40, mu = mu, Sigma = Sigma1)
X2 <- mvrnorm(n = 40, mu = mu, Sigma = Sigma2)
X3 <- mvrnorm(n = 40, mu = mu, Sigma = Sigma3)


plot(X1, pch = 19)
plot(X2, pch = 19)
plot(X3, pch = 19)
            import numpy as np
import matplotlib.pyplot as plt

mu = np.array([2, 1])

Sigma1 = np.array([[5,  4],
                [4,  5]])
Sigma2 = np.diag([3, 3])  
Sigma3 = np.array([[5, -4],
                [-4,  5]])

np.random.seed(0)

X1 = np.random.multivariate_normal(mean=mu, cov=Sigma1, size=40)
X2 = np.random.multivariate_normal(mean=mu, cov=Sigma2, size=40)
X3 = np.random.multivariate_normal(mean=mu, cov=Sigma3, size=40)
for X, name in [(X1, "X1"), (X2, "X2"), (X3, "X3")]:
    plt.figure()
    plt.scatter(X[:, 0], X[:, 1], marker='o')  # pch = 19 ≈ 实心圆
    plt.xlabel("x")
    plt.ylabel("y")
    plt.title(name)
    plt.show()

          • 三维正态分布密度图：

              library(mvtnorm)
# 多元正态分布的均值
mu1 <- 0.5
mu2 <- 0.5
# 对应书上 Figure 4.2(a) 的\rho_{12}，这里意味着 X1 和 X2 独立
rho1 <- 0
# 对应书上 Figure 4.2(b) 的\rho_{12}
rho2 <- 0.75
# 对应书上\sigma_{11} = \sigma_{22}
sigma <- 0.2
# 协方差矩阵的维数
p <- 2
# 构造协方差矩阵，用 rho 来表示
Sigma1 <- (diag(c(1-rho1, 1-rho1)) + matrix(rep(rho1, p^2), ncol = p)) * sigma
Sigma2 <- (diag(c(1-rho2, 1-rho2)) + matrix(rep(rho2, p^2), ncol = p)) * sigma
# 张成一个 40x40 的网格
x1 <- seq(from = mu1 - 7*sigma, to = mu1 + 7*sigma, length.out = 40)
x2 <- seq(from = mu2 - 7*sigma, to = mu2 + 7*sigma, length.out = 40)
## 将 x1 和 x2" 编织" 成网格 (可类比 JMP 中的表> 堆叠操作）
## 对于 x2 的每一个取值都有 x1 的所有 40 个取值与之对应
## 相当于把网格里的每一个点（二维）都枚举了出来，每一行就是一个点
X <- expand.grid(x1 = x1, x2 = x2)

dX1 <- dmvnorm(X, mean = c(mu1, mu2), sigma = Sigma1)
dX2 <- dmvnorm(X, mean = c(mu1, mu2), sigma = Sigma2)
## 将一维的向量转换成矩阵
dX1_ <- matrix(dX1, nrow = 40, ncol = 40)
dX2_ <- matrix(dX2, nrow = 40, ncol = 40)
# 画出三维透视图
persp(x = x1, y = x2, z = dX1_, theta = 30, phi = 30)
persp(x = x1, y = x2, z = dX2_, theta = 30, phi = 30)

# 等高线图
contour(x = x1, y = x2, z = dX1_, levels = c(0.1, 0.5)) 
contour(x = x1, y = x2, z = dX2_, levels = c(0.1, 0.5))

# 绘制样本椭圆（python中无对应）
library(car)
dataEllipse(X1, levels = c(0.95, 0.5))
              import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

mu1 = 0.5
mu2 = 0.5
rho1 = 0
rho2 = 0.75
sigma = 0.2
p = 2

Sigma1 = (np.diag([1 - rho1, 1 - rho1]) + np.full((p, p), rho1)) * sigma
Sigma2 = (np.diag([1 - rho2, 1 - rho2]) + np.full((p, p), rho2)) * sigma

x1 = np.linspace(mu1 - 7 * sigma, mu1 + 7 * sigma, 40)
x2 = np.linspace(mu2 - 7 * sigma, mu2 + 7 * sigma, 40)

X1, X2 = np.meshgrid(x1, x2)
pos = np.dstack((X1, X2))

dX1 = multivariate_normal.pdf(pos, mean=[mu1, mu2], cov=Sigma1)
dX2 = multivariate_normal.pdf(pos, mean=[mu1, mu2], cov=Sigma2)

fig1 = plt.figure()
ax1 = fig1.add_subplot(111, projection='3d')
ax1.plot_surface(X1, X2, dX1, cmap='viridis')
ax1.view_init(elev=30, azim=30)

fig2 = plt.figure()
ax2 = fig2.add_subplot(111, projection='3d')
ax2.plot_surface(X1, X2, dX2, cmap='viridis')
ax2.view_init(elev=30, azim=30)

plt.show()

plt.contour(x1, x2, dX1, levels=[0.1, 0.5])
plt.show()

plt.contour(x1, x2, dX2, levels=[0.1, 0.5])
plt.show()

            正态性检验

            1. Q-Q plot绘图检验

              • 将样本x1,x2,…,xnx_1,x_2,\dots,x_n从小到大排序得到x(1),x(2),…,x(n)x_{(1)},x_{(2)},\dots,x_{(n)}，它们各自对应的概率值为1−12n,2−12n,…,n−12n\displaystyle\frac{1-\frac{1}{2}}{n},\frac{2-\frac{1}{2}}{n},\dots,\frac{n-\frac{1}{2}}{n}。
              • 计算标准正态分布的下分位数q(1),q(2),…,q(n)q_{(1)},q_{(2)},\dots,q_{(n)}。
              • 画出(q(1),x(1)),(q(2),x(2)),…,(q(n),x(n))(q_{(1)},x_{(1)}),(q_{(2)},x_{(2)}),\dots,(q_{(n)},x_{(n)})的图，并检查这些点是否都近似地在一条直线上。

              当然也可以直接用qqnorm()函数，此时概率值变为1−38n−14,2−38n−14,…,n−38n−14\displaystyle \frac{1-\frac{3}{8}}{n-\frac{1}{4}},\frac{2-\frac{3}{8}}{n-\frac{1}{4}},\dots,\frac{n-\frac{3}{8}}{n-\frac{1}{4}}。

              • 同时也可以用qqline()绘制参考线（基于样本四分位数与正态分布四分位数点连线得到）

              • 更严谨的方法是计算Q-Q plot的相关系数，可参考下述检验表格：

                n .01 .05 .10
                5 .8299 .8788 .9032
                10 .8801 .9198 .9351
                15 .9126 .9389 .9503

                当指定样本数与显著性水平下相关系数小于临界值，则拒绝正态性假设。

                # 输入原始观测值
x <- c(-1, -0.1, 0.16, 0.41, 0.62, 0.8, 1.26, 1.54, 1.71, 2.3)
## 使用 sort() 函数进行排序
x_sorted <- sort(x)
## 代入公式，计算样本各自对应的概率值
prob_values <- (rank(x) - 1/2) / length(x)
q <- qnorm(prob_values)

# 利用 qqnorm() 函数内置的计算方式绘制出的点
qqnorm(x, pch = 19)
# 利用我们自己写的代码绘制出的点（蓝色）
points(q, x, pch = 19, col = "blue")

qqline(x)

cor(x, q) # 相关系数检验正态性
                import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

x = np.array([-1, -0.1, 0.16, 0.41, 0.62, 0.8, 1.26, 1.54, 1.71, 2.3])
n = len(x)

prob_qqnorm = (np.arange(1, n + 1) - 3/8) / (n + 1/4)
q_qqnorm = stats.norm.ppf(prob_qqnorm)
x_sorted = np.sort(x)

rank = stats.rankdata(x)
prob_values = (rank - 0.5) / n
q = stats.norm.ppf(prob_values)

fig, ax = plt.subplots(figsize=(6, 6))

ax.scatter(q_qqnorm, x_sorted, color='black', label='qqnorm()', zorder=3)
ax.scatter(q, x, color='blue', label='custom (rank-1/2)/n', zorder=3)

probs = np.array([0.25, 0.75])
y_q = np.quantile(x, probs, method='linear')
x_q = stats.norm.ppf(probs)
slope = (y_q[1] - y_q[0]) / (x_q[1] - x_q[0])
intercept = y_q[0] - slope * x_q[0]

xx = np.linspace(q_qqnorm.min(), q_qqnorm.max(), 100)
ax.plot(xx, intercept + slope * xx, color='red', label='qqline')

ax.set_xlabel('Theoretical Quantiles')
ax.set_ylabel('Sample Quantiles')
ax.set_title('Normal Q-Q Plot')
ax.legend()
plt.show()

r = np.corrcoef(x_sorted, q_qqnorm)[0, 1]
print(r)
                • R语言中也可以用shapiro.test(x)进行正态性检验。（具体略）

              • 多变量正态性检验
                使用Chi-square plot绘图检验。具体方法：计算马氏距离平方

                dj2=(xj−x‾)⊤S−1(xj−x‾),j=1,2,…,nd_j^2=(\mathbf{x}_j-\overline{\mathbf{x}})^\top\mathbf{S}^{-1}(\mathbf{x}_j-\overline{\mathbf{x}}),\quad j=1,2,\dots,n

                并把它们从小到大排序得到d(1)2,d(2)2,…,d(n)2d_{(1)}^2,d_{(2)}^2,\dots,d_{(n)}^2。然后画出

                (qc,p(j−12n),d(j)2)\left(q_{c,p}\left(\frac{j-\frac{1}{2}}{n}\right),d_{(j)}^2\right)

                这里qc,p(j−12n)\displaystyle q_{c,p}\left(\frac{j-\frac{1}{2}}{n}\right)是χp2\chi_p^2的100⋅j−12n100\cdot\dfrac{j-\frac{1}{2}}{n}下分位数。

                  # 读入数据
x1 <- c(108.28, 152.36, 95.04, 65.45, 62.97, 263.99, 265.19, 285.06, 92.01, 165.68)
x2 <- c(17.05, 16.59, 10.91, 14.14, 9.52, 25.33, 18.54, 15.73, 8.10, 11.13)
X <- cbind(x1, x2)
# 计算 d_j^2
d2 <- mahalanobis(X,
                center = colMeans(X), # 样本的均值向量
                cov = cov(X)) # 样本的协方差矩阵
# 1. 排序（从小到大）
d2.sorted <- sort(d2)
# 2. 绘图
d2.q <- qchisq((rank(d2.sorted) - 1/2) / length(d2.sorted),
            ncol(X))
plot(d2.q, d2.sorted, pch = 19, col = "blue")

library(MVN)
library(ggplot2)

p <- multivariate_diagnostic_plot(X, type = "qq")

# 在返回的 ggplot 对象上叠加散点
p + geom_point(aes(x = d2.sorted, y = d2.q), 
            color = "blue", 
            shape = 19)
                  import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

x1 = np.array([108.28, 152.36, 95.04, 65.45, 62.97, 263.99, 265.19, 285.06, 92.01, 165.68])
x2 = np.array([17.05, 16.59, 10.91, 14.14, 9.52, 25.33, 18.54, 15.73, 8.10, 11.13])
X = np.column_stack([x1, x2])

mean = X.mean(axis=0)
cov = np.cov(X, rowvar=False)
inv_cov = np.linalg.inv(cov)
diff = X - mean
d2 = np.einsum('ij,jk,ik->i', diff, inv_cov, diff)

d2_sorted = np.sort(d2)
n = len(d2_sorted)

probs = (np.arange(1, n + 1) - 0.5) / n
d2_q = stats.chi2.ppf(probs, df=X.shape[1])

plt.figure(figsize=(6, 5))
plt.scatter(d2_q, d2_sorted, color='blue', marker='o')
plt.plot([0, max(d2_q.max(), d2_sorted.max())],
        [0, max(d2_q.max(), d2_sorted.max())], 'r--', lw=1)
plt.xlabel('Chi-square quantiles')
plt.ylabel('Squared Mahalanobis distance')
plt.title('Chi-square Q-Q plot')
plt.tight_layout()
plt.show()

                异常值检测

                检测异常点的步骤如下：

                1. 为每个变量画 dot plot，关注离群点。
                2. 为每对变量画 scatter plot，关注离群点。
                3. 计算zjk=xjk−x‾kskk,j=1,2,…,n,k=1,2,…,pz_{jk}=\dfrac{x_{jk}-\overline{x}_k}{\sqrt{s_{kk}}},j=1,2,\dots,n,k=1,2,\dots,p，关注绝对值大的zjkz_{jk}。
                4. 计算(xj−x‾)⊤S−1(xj−x‾),j=1,2,…,n(\mathbf{x}_j-\overline{\mathbf{x}})^\top\mathbf{S}^{-1}(\mathbf{x}_j-\overline{\mathbf{x}}),j=1,2,\dots,n，关注离原点远的点。
                  data4.3 <- read.table(text = "
1889  1651  1561  1778  .60
2403  2048  2087  2197  5.48
2119  1700  1815  2222  7.62
1645  1627  1110  1533  5.21
1976  1916  1614  1883  1.40
1712  1712  1439  1546  2.22
1943  1685  1271  1671  4.99
2104  1820  1717  1874  1.49
2983  2794  2412  2581  12.26
1745  1600  1384  1508  .77
1710  1591  1518  1667  1.93
2046  1907  1627  1898  .46
1840  1841  1595  1741  2.70
1867  1685  1493  1678  .13
1859  1649  1389  1714  1.08
1954  2149  1180  1281  16.85
1325  1170  1002  1176  3.50
1419  1371  1252  1308  3.99
1828  1634  1602  1755  1.36
1725  1594  1313  1646  1.46
2276  2189  1547  2111  9.90
1899  1614  1422  1477  5.06
1633  1513  1290  1516  .80
2061  1867  1646  2037  2.54
1856  1493  1356  1533  4.58
1727  1412  1238  1469  3.40
2168  1896  1701  1834  2.38
1655  1675  1414  1597  3.00
2326  2301  2065  2234  6.28
1490  1382  1214  1284  2.58
", header = FALSE)

X <- as.matrix(data4.3[, -5])

# 为每个变量画 dot plot，关注离群点
stripchart(data4.3[, -5], 
        method = "stack", 
        pch = 19, 
        col = c("brown", "blue", "green4", "purple"))

# 为每对变量画 scatter plot，关注离群点
pairs(data4.3[, -5])

# 计算标准化变量 z
Z <- scale(X, center = TRUE, scale = TRUE)
# 计算马氏距离
d2 <- mahalanobis(X, center = colMeans(X), cov = cov(X))

cbind(Z, d2)
                  import pandas as pd
import numpy as np
import io
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import zscore
from scipy.spatial.distance import mahalanobis

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

data_str = """
    1889  1651  1561  1778  .60
    2403  2048  2087  2197  5.48
    2119  1700  1815  2222  7.62
    1645  1627  1110  1533  5.21
    1976  1916  1614  1883  1.40
    1712  1712  1439  1546  2.22
    1943  1685  1271  1671  4.99
    2104  1820  1717  1874  1.49
    2983  2794  2412  2581  12.26
    1745  1600  1384  1508  .77
    1710  1591  1518  1667  1.93
    2046  1907  1627  1898  .46
    1840  1841  1595  1741  2.70
    1867  1685  1493  1678  .13
    1859  1649  1389  1714  1.08
    1954  2149  1180  1281  16.85
    1325  1170  1002  1176  3.50
    1419  1371  1252  1308  3.99
    1828  1634  1602  1755  1.36
    1725  1594  1313  1646  1.46
    2276  2189  1547  2111  9.90
    1899  1614  1422  1477  5.06
    1633  1513  1290  1516  .80
    2061  1867  1646  2037  2.54
    1856  1493  1356  1533  4.58
    1727  1412  1238  1469  3.40
    2168  1896  1701  1834  2.38
    1655  1675  1414  1597  3.00
    2326  2301  2065  2234  6.28
    1490  1382  1214  1284  2.58
"""

data4_3 = pd.read_csv(io.StringIO(data_str), sep=r'\s+', header=None)
data4_3.columns = ['V1', 'V2', 'V3', 'V4', 'V5']
X = data4_3.iloc[:, :-1].values

df_plot = data4_3.iloc[:, :-1]
df_melted = df_plot.melt(var_name='Variable', value_name='Value')

plt.figure(figsize=(8, 5))
sns.stripplot(data=df_melted, x='Value', y='Variable', 
            palette=["brown", "blue", "green4", "purple"], 
            size=6, jitter=0, alpha=0.8)
plt.gca().invert_yaxis()
plt.title("Stripchart (Dot plot)")
plt.show()

sns.pairplot(df_plot, diag_kind='kde', plot_kws={'alpha': 0.6})
plt.suptitle("Scatter Plot Matrix (Pairs)", y=1.02)
plt.show()

Z = zscore(X, axis=0, ddof=1)

mu = np.mean(X, axis=0)
cov_matrix = np.cov(X, rowvar=False)
inv_cov_matrix = np.linalg.inv(cov_matrix)

d2 = np.array([mahalanobis(row, mu, inv_cov_matrix) ** 2 for row in X])

result_df = pd.DataFrame(Z, columns=['Z_V1', 'Z_V2', 'Z_V3', 'Z_V4'])
result_df['d2'] = d2

pd.set_option('display.max_rows', None)
print(result_df)

                  近似正态变换

                  1. 单变量情形

                      x <- c(
0.15, 0.09, 0.18, 0.10, 0.05, 0.12, 0.08, 0.05,
0.08, 0.10, 0.07, 0.02, 0.01, 0.10, 0.10, 0.10,
0.02, 0.10, 0.01, 0.40, 0.10, 0.05, 0.03, 0.05,
0.15, 0.10, 0.15, 0.09, 0.08, 0.18, 0.10, 0.20,
0.11, 0.30, 0.02, 0.20, 0.20, 0.30, 0.30, 0.40,
0.30, 0.05
)

qqnorm(x, pch = 19)
qqline(x)

## 计算 lambda
library(car)
best_lambda <- powerTransform(x, family = "bcPower")
best_lambda

library(forecast)
x_bc <- BoxCox(x, lambda = 0.276) 
qqnorm(x_bc, pch = 19)
qqline(x_bc)
                      import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

x = np.array([
    0.15, 0.09, 0.18, 0.10, 0.05, 0.12, 0.08, 0.05,
    0.08, 0.10, 0.07, 0.02, 0.01, 0.10, 0.10, 0.10,
    0.02, 0.10, 0.01, 0.40, 0.10, 0.05, 0.03, 0.05,
    0.15, 0.10, 0.15, 0.09, 0.08, 0.18, 0.10, 0.20,
    0.11, 0.30, 0.02, 0.20, 0.20, 0.30, 0.30, 0.40,
    0.30, 0.05
])

fig, ax = plt.subplots(figsize=(6, 4))
stats.probplot(x, dist="norm", plot=ax)
ax.set_title("Normal Q-Q Plot")
ax.set_xlabel("Theoretical Quantiles")
ax.set_ylabel("Sample Quantiles")
plt.show()

best_lambda = stats.boxcox_normmax(x, method="mle")
print("Estimated lambda:", best_lambda)

def boxcox_forecast(x, lam):
    if np.isclose(lam, 0):
        return np.log(x)
    return (x**lam - 1) / lam

x_bc = boxcox_forecast(x, best_lambda)

fig, ax = plt.subplots(figsize=(6, 4))
stats.probplot(x_bc, dist="norm", plot=ax)
ax.set_title("Normal Q-Q Plot (Box-Cox transformed)")
ax.set_xlabel("Theoretical Quantiles")
ax.set_ylabel("Sample Quantiles")
plt.show()

                    • 多元情形

                        x1 <- x
x2 <- c(
0.30, 0.09, 0.30, 0.10, 0.10, 0.12, 0.09, 0.10, 0.09, 0.10,
0.07, 0.05, 0.01, 0.45, 0.12, 0.20, 0.04, 0.10, 0.01, 0.60,
0.12, 0.10, 0.05, 0.05, 0.15, 0.30, 0.15, 0.09, 0.09, 0.28,
0.10, 0.10, 0.10, 0.30, 0.12, 0.25, 0.20, 0.40, 0.33, 0.32,
0.12, 0.12
)
X <- cbind(x1,x2)
# 从联合分布的角度，选 lambda
powerTransform(X, family = "bcPower")

# 从边际分布的角度，选 lambda
powerTransform(x1, family = "bcPower")
powerTransform(x2, family = "bcPower")
                        import numpy as np
from scipy.optimize import minimize
from scipy import stats

x = np.array([
    0.15, 0.09, 0.18, 0.10, 0.05, 0.12, 0.08, 0.05,
    0.08, 0.10, 0.07, 0.02, 0.01, 0.10, 0.10, 0.10,
    0.02, 0.10, 0.01, 0.40, 0.10, 0.05, 0.03, 0.05,
    0.15, 0.10, 0.15, 0.09, 0.08, 0.18, 0.10, 0.20,
    0.11, 0.30, 0.02, 0.20, 0.20, 0.30, 0.30, 0.40,
    0.30, 0.05
])

x1 = x.copy()
x2 = np.array([
    0.30, 0.09, 0.30, 0.10, 0.10, 0.12, 0.09, 0.10, 0.09, 0.10,
    0.07, 0.05, 0.01, 0.45, 0.12, 0.20, 0.04, 0.10, 0.01, 0.60,
    0.12, 0.10, 0.05, 0.05, 0.15, 0.30, 0.15, 0.09, 0.09, 0.28,
    0.10, 0.10, 0.10, 0.30, 0.12, 0.25, 0.20, 0.40, 0.33, 0.32,
    0.12, 0.12
])

X = np.column_stack((x1, x2))

def boxcox_transform(x, lam):
    if np.isclose(lam, 0):
        return np.log(x)
    return (x**lam - 1) / lam

def neg_log_likelihood_multivariate(lambdas, X):
    n, p = X.shape
    Y = np.empty_like(X, dtype=float)
    log_jac = 0.0
    for j in range(p):
        lam = lambdas[j]
        Y[:, j] = boxcox_transform(X[:, j], lam)
        log_jac += (lam - 1) * np.sum(np.log(X[:, j]))
    Y_centered = Y - Y.mean(axis=0)
    Sigma = (Y_centered.T @ Y_centered) / n
    sign, logdet = np.linalg.slogdet(Sigma)
    nll = (n / 2) * logdet - log_jac
    return nll

lam1_init = stats.boxcox_normmax(x1, method="mle")
lam2_init = stats.boxcox_normmax(x2, method="mle")
res = minimize(neg_log_likelihood_multivariate, [lam1_init, lam2_init], args=(X,), method='Nelder-Mead')
print("Joint lambdas:", res.x)

print("Marginal lambda x1:", stats.boxcox_normmax(x1, method="mle"))
print("Marginal lambda x2:", stats.boxcox_normmax(x2, method="mle"))