多元统计分析 · 大题练习 8题

import pandas as pd
import scipy.cluster.hierarchy as sch

dat9 = pd.read_clipboard()

# 1. 欧式距离
dat9_D = sch.distance.pdist(dat9)

# 2. 离差平方和法（ward）系统聚类
dat9_H = sch.linkage(dat9_D, method='ward')

# 3. 树图
import matplotlib.pyplot as plt
sch.dendrogram(dat9_H)
plt.show()
# 根据树图 → 分3类

# 4. 碎石图确定K
from sklearn.cluster import KMeans
k_list = [1,2,3,4,5]
sse_list = []
for k in k_list:
    km = KMeans(n_clusters=k, n_init=10).fit(dat9)
    sse_list.append(km.inertia_)
plt.plot(k_list, sse_list, '-o')
plt.xlabel('K'); plt.ylabel('SSE'); plt.show()
# 根据碎石图 → K=2或3

# 5. K=3 K-means
dat9_KM = KMeans(n_clusters=3, n_init=10).fit(dat9).predict(dat9)

# 6. 输出结果，检查纽约与加州
dat9_class = pd.DataFrame({'类别': dat9_KM+1}, index=dat9.index)
print(dat9_class[dat9_class.类别==1])  # 纽约与加州同在类别1

dat35 = pd.read_clipboard()

# 1. 欧式距离
dat35_D = sch.distance.pdist(dat35)

# 2. 离差平方和法系统聚类
dat35_H = sch.linkage(dat35_D, method='ward')

# 3. 树图
sch.dendrogram(dat35_H)
plt.title("Dat35 系统聚类树图"); plt.show()

# 4. 碎石图
k_list = [1,2,3,4,5]
sse_list = []
for k in k_list:
    km = KMeans(n_clusters=k, n_init=10).fit(dat35)
    sse_list.append(km.inertia_)
plt.plot(k_list, sse_list, '-o')
plt.xlabel("K"); plt.ylabel("SSE"); plt.show()

# 5. K=3 K-means
dat35_KM = KMeans(n_clusters=3, n_init=10).fit(dat35).predict(dat35)

# 6. 北京所在类别
dat35_class = pd.DataFrame({'类别': dat35_KM+1}, index=dat35.index)
bj_class = dat35_class.loc['北京', '类别']
print("与北京同类的地区:", dat35_class[dat35_class.类别==bj_class].index.tolist())

d93 = pd.read_clipboard()

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as lda
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis as qda

X = d93[['性别','年龄','月收入']]
Y = d93.y

# 1. Fisher线性判别
d93_lda = lda().fit(X, Y)

# 2. 回判
d93['ld_G'] = d93_lda.predict(X)
d93_ld_M = d93.pivot_table('序号','y','ld_G', aggfunc=len, margins=True)
print('混淆矩阵:\n', d93_ld_M)
print(f'Fisher正确率: {sum(d93.y==d93.ld_G)/len(d93):.4f}')  # 0.75

# 3. 非线性判别
d93_qda = qda().fit(X, Y)

# 4. 回判
d93['qd_G'] = d93_qda.predict(X)
d93_qd_M = d93.pivot_table('序号','y','qd_G', aggfunc=len, margins=True)
print(d93_qd_M)
print(f'非线性正确率: {sum(d93.y==d93.qd_G)/len(d93):.4f}')  # 0.8214

# 5. 结论：非线性判别更适合，正确率更高（82.14% > 75%）

dat121 = pd.read_clipboard()

X = dat121[['sepal.length','sepal.width','petal.length','petal.width']]
Y = dat121.Species

# 1. Fisher线性判别
dat121_lda = lda().fit(X, Y)

# 2. 回判
dat121['ld_G'] = dat121_lda.predict(X)
dat121_ld_M = dat121.pivot_table('Species','Species','ld_G', aggfunc=len, margins=True)
print('混淆矩阵:\n', dat121_ld_M)
correct = sum(dat121['Species'] == dat121['ld_G'])
print(f'Fisher正确率: {correct}/{len(dat121)} = {correct/len(dat121):.4f}')

# 3. 非线性判别
dat121_qda = qda().fit(X, Y)

# 4. 回判
dat121['qd_G'] = dat121_qda.predict(X)
dat121_qd_M = dat121.pivot_table('Species','Species','qd_G', aggfunc=len, margins=True)
print(dat121_qd_M)
correct_qda = sum(dat121['Species'] == dat121['qd_G'])
print(f'非线性正确率: {correct_qda}/{len(dat121)} = {correct_qda/len(dat121):.4f}')

# 5. 比较两个正确率，更高的更适合

from sklearn.decomposition import PCA
import numpy as np

d61 = pd.read_clipboard()

# 主成分评价函数
def PCscores(X, m=2):
    Z = (X - X.mean()) / X.std()
    p = Z.shape[1]
    pca = PCA(n_components=p).fit(Z)
    Vi = pca.explained_variance_
    Wi = pca.explained_variance_ratio_
    Vars = pd.DataFrame({'Variances': Vi}, index=[f'Comp{i+1}' for i in range(p)])
    Vars['Explained'] = Wi * 100
    Vars['Cumulative'] = np.cumsum(Wi) * 100
    print("\n方差贡献:\n", round(Vars, 4))
    Compi = [f'Comp{i+1}' for i in range(m)]
    loadings = pd.DataFrame(pca.components_[:m].T, columns=Compi, index=X.columns)
    print("\n主成分负荷:\n", round(loadings, 4))
    scores = pd.DataFrame(pca.fit_transform(Z)).iloc[:, :m]
    scores.index = X.index; scores.columns = Compi
    scores['Comp'] = scores.dot(Wi[:m])
    scores['Rank'] = scores['Comp'].rank(ascending=False).astype(int)
    return scores

# 1. 主成分分析
d61_pcs = PCscores(d61)
print(d61_pcs)

# 2. Comp1 49.38% + Comp2 21.74% + Comp3 13.01% = 84.14% > 80% → 选3个

# 3. 第2主成分方差贡献率: 21.74%

# 4. 第二主成分表达式（见负荷表Comp2列）

# 5. 信息图
def Scoreplot(Scores):
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    plt.plot(Scores.iloc[:,0], Scores.iloc[:,1], '*')
    plt.xlabel(Scores.columns[0]); plt.ylabel(Scores.columns[1])
    plt.axhline(y=0, ls=':'); plt.axvline(x=0, ls=':')
    for i in range(len(Scores)):
        plt.text(Scores.iloc[i,0], Scores.iloc[i,1], Scores.index[i])

Scoreplot(d61_pcs); plt.show()
# 得分较高：长江投资、铁龙物流、交运股份

# 6. 前2维仅71.12%，不足80%，说服力不够充分

dat52 = pd.read_clipboard()

# 1. 主成分分析
dat52_pcs = PCscores(dat52)
print(dat52_pcs)

# 2. 观察 Cumulative 列，找到超过80%的最小主成分数

# 3. 第二主成分表达式（写Comp2负荷系数）

# 4. 信息图
Scoreplot(dat52_pcs); plt.show()

# 5. 前2维累计贡献率若不够80%，需要更多主成分

# 6. 查看方差贡献表中 Comp2 的 Explained 值

dat52 = pd.read_clipboard()

import factor_analyzer as fa
from factor_analyzer import FactorAnalyzer as FA

# 1. KMO和Bartlett检验
kmo = fa.calculate_kmo(dat52)
print(f'KMO: {kmo[1]:.4f}')  # 0.7499，接近1，适合
chisq = fa.calculate_bartlett_sphericity(dat52)
print(f'卡方值={chisq[0]:.4f}, p值={chisq[1]:.4f}')  # p=0.0000，变量关系显著
# → 适合做因子分析

# 2. 8因子最大似然法
Fm1 = FA(n_factors=8, method='ml', rotation=None).fit(dat52.values)

# 3. 共同度
communalities = 1 - Fm1.get_uniquenesses()
print('共同度:', communalities)
# x7共同度仅68.6%，未被很好解释

# 4. 方差贡献
def Factors(fa): return [f'F{i}' for i in range(1, fa.n_factors+1)]
Vars_names = ['方差','贡献率','累计贡献率']
Fm1_Vars = pd.DataFrame(Fm1.get_factor_variance(), Vars_names, Factors(Fm1))
print(Fm1_Vars)
# F1累计贡献率86.5%>80%，选1个因子即可

# 5. 2因子varimax旋转
Fp2 = FA(2, method='principal', rotation='varimax').fit(dat52.values)
Fp2_load = pd.DataFrame(Fp2.loadings_, dat52.columns, Factors(Fp2))
print('旋转后载荷:\n', Fp2_load)
# F1: X1-X6,X8（高载荷）; F2: X7（高载荷）

# 6. 因子解释：F1代表综合经营效益，F2代表固定资产周转

# 因子得分函数
def FArank(Vars, Scores):
    Vi = Vars.values[0]
    Wi = Vi / sum(Vi)
    Fi = Scores.dot(Wi)
    Ri = Fi.rank(ascending=False).astype(int)
    return pd.DataFrame({'因子得分': Fi, '因子排名': Ri})

# 7. 因子得分
Fp2_scores = pd.DataFrame(Fp2.transform(dat52.values), dat52.index, Factors(Fp2))
Fp2_Vars = pd.DataFrame(Fp2.get_factor_variance(), Vars_names, Factors(Fp2))
print(FArank(Fp2_Vars, Fp2_scores))
# → 江苏排名最高，应投江苏

d61 = pd.read_clipboard()

import factor_analyzer as fa
from factor_analyzer import FactorAnalyzer as FA

# 1. KMO和Bartlett检验
kmo = fa.calculate_kmo(d61)
print(f'KMO: {kmo[1]:.4f}')
chisq = fa.calculate_bartlett_sphericity(d61)
print(f'卡方值={chisq[0]:.4f}, p值={chisq[1]:.4f}')
# KMO接近1且p<0.05 → 适合做因子分析

# 2. 4因子最大似然法
Fm4 = FA(n_factors=4, method='ml', rotation=None).fit(d61.values)

# 3. 共同度
communalities = 1 - Fm4.get_uniquenesses()
print('共同度:', communalities)

# 4. 方差贡献
Vars_df = pd.DataFrame(Fm4.get_factor_variance(),
    ['方差','贡献率','累计贡献率'],
    [f'F{i}' for i in range(1,5)])
print('\n方差贡献:\n', Vars_df)
# 选累计贡献率>80%的因子数

# 5. 2因子varimax旋转
F2 = FA(2, method='principal', rotation='varimax').fit(d61.values)
F2_load = pd.DataFrame(F2.loadings_, d61.columns, ['F1','F2'])
print('\n旋转后载荷:\n', F2_load)

# 6. 根据载荷解释因子含义

# 7. 因子得分
F2_scores = pd.DataFrame(F2.transform(d61.values), d61.index, ['F1','F2'])
F2_vars = pd.DataFrame(F2.get_factor_variance(), ['方差','贡献率','累计贡献率'], ['F1','F2'])
print('\n因子得分与排名:\n', FArank(F2_vars, F2_scores))
# 得分最高者推荐投资