【数学建模】 机器学习与统计模型

机器学习与统计模型

1. 统计分布与假设检验

1.1 统计量与常见统计分布

常见统计分布

1. 正态分布（Normal Distribution）：正态分布是一种连续概率分布，广泛用于自然科学和社会科学。其概率密度函数（PDF）定义如下：
   $$f(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$
   其中，$\mu$ 是均值，$\sigma$ 是标准差。正态分布的图像呈钟形曲线，均值$\mu$决定曲线的中心位置，标准差$\sigma$决定曲线的宽度。

2. 二项分布（Binomial Distribution）：二项分布是离散概率分布，表示在 $n$ 次独立试验中成功次数的分布，每次试验成功的概率为 $p$。其概率质量函数（PMF）定义如下：
   $$P(X=k)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^k(1-p{)}^{n-k}$$
   其中，$k$ 是成功次数，$n$ 是试验次数，$p$ 是成功的概率，$\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ 是二项系数。

3. 泊松分布（Poisson Distribution）：泊松分布是离散概率分布，用于描述单位时间内随机事件发生次数。其概率质量函数（PMF）定义如下：
   $$P(X=k)=\frac{{\lambda }^k{e}^{-\lambda }}{k!}$$
   其中，$\lambda$ 是单位时间内事件发生的平均次数，$k$ 是事件发生的次数。

Python代码示例

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm, binom, poisson

# 正态分布
mu, sigma = 0, 1
x = np.linspace(-5, 5, 1000)
plt.plot(x, norm.pdf(x, mu, sigma), label='Normal Distribution')

# 二项分布
n, p = 10, 0.5
x_binom = np.arange(0, n + 1)
plt.bar(x_binom, binom.pmf(x_binom, n, p), alpha=0.5, label='Binomial Distribution')

# 泊松分布
lambda_ = 3
x_poisson = np.arange(0, 15)
plt.bar(x_poisson, poisson.pmf(x_poisson, lambda_), alpha=0.5, label='Poisson Distribution')

plt.legend()
plt.xlabel('Value')
plt.ylabel('Probability')
plt.title('Common Probability Distributions')
plt.show()

在上面的代码中，我们绘制了正态分布、二项分布和泊松分布的概率分布图。正态分布使用 scipy.stats.norm.pdf 计算概率密度函数，二项分布和泊松分布分别使用 scipy.stats.binom.pmf 和 scipy.stats.poisson.pmf 计算概率质量函数。

1.2 正态性检验

正态性检验用于判断样本是否来自正态分布。常用的方法包括Shapiro-Wilk检验和Kolmogorov-Smirnov检验。

Shapiro-Wilk检验

Shapiro-Wilk检验用于检验样本数据是否符合正态分布。其检验假设为：

• $H_0$：样本来自正态分布
• $H_1$：样本不来自正态分布

Python代码示例

from scipy.stats import shapiro

data = np.random.normal(loc=0, scale=1, size=100)
stat, p = shapiro(data)
print('Statistics=%.3f, p=%.3f' % (stat, p))

if p > 0.05:
    print('Sample looks Gaussian (fail to reject H0)')
else:
    print('Sample does not look Gaussian (reject H0)')

在这段代码中，我们生成了一组正态分布的数据，并使用 Shapiro-Wilk 检验进行正态性检验。如果 p 值大于 0.05，我们无法拒绝原假设，认为样本数据符合正态分布；否则，拒绝原假设，认为样本数据不符合正态分布。

1.3 独立性检验

独立性检验用于判断两个变量是否独立，常用方法有卡方检验。

卡方检验

卡方检验用于判断两个分类变量之间是否存在关联。其检验假设为：

• $H_0$：变量独立
• $H_1$：变量不独立

Python代码示例

import pandas as pd
from scipy.stats import chi2_contingency

# 构造列联表
data = [[10, 20, 30], [6, 9, 17]]
df = pd.DataFrame(data, columns=['A', 'B', 'C'])
stat, p, dof, expected = chi2_contingency(df)
print('Statistics=%.3f, p=%.3f' % (stat, p))

if p > 0.05:
    print('Variables are independent (fail to reject H0)')
else:
    print('Variables are not independent (reject H0)')

在这段代码中，我们创建了一个列联表，并使用卡方检验判断变量之间是否独立。如果 p 值大于 0.05，我们无法拒绝原假设，认为变量之间独立；否则，拒绝原假设，认为变量之间不独立。

1.4 两组样本的差异性检验

常见的方法有t检验（包括独立样本t检验和配对样本t检验）。

独立样本t检验

独立样本t检验用于比较两组独立样本的均值是否相等。其检验假设为：

• $H_0$：两组样本均值相等
• $H_1$：两组样本均值不相等

Python代码示例

from scipy.stats import ttest_ind

# 两组样本
data1 = np.random.normal(0, 1, 100)
data2 = np.random.normal(0.5, 1, 100)

stat, p = ttest_ind(data1, data2)
print('Statistics=%.3f, p=%.3f' % (stat, p))

if p > 0.05:
    print('Same distributions (fail to reject H0)')
else:
    print('Different distributions (reject H0)')

在这段代码中，我们生成了两组正态分布的数据，并使用独立样本t检验比较两组样本的均值。如果 p 值大于 0.05，我们无法拒绝原假设，认为两组样本均值相等；否则，拒绝原假设，认为两组样本均值不相等。

1.5 方差分析与事后多重比较

方差分析（ANOVA）用于比较三组或更多组的均值是否相等。

方差分析

方差分析用于检验三组或更多组样本的均值是否相等。其检验假设为：

• $H_0$：各组均值相等
• $H_1$：至少有一组均值不等

Python代码示例

from scipy.stats import f_oneway

# 三组样本
data1 = np.random.normal(0, 1, 100)
data2 = np.random.normal(0.5, 1, 100)
data3 = np.random.normal(1, 1, 100)

stat, p = f_oneway(data1, data2, data3)
print('Statistics=%.3f, p=%.3f' % (stat, p))

if p > 0.05:
    print('Same means (fail to reject H0)')
else:
    print('Different means (reject H0)')

在这段代码中，我们生成了三组正态分布的数据，并使用方差分析检验三组样本的均值是否相等。如果 p 值大于 0.05，我们无法拒绝原假设，认为各组样本均值相等；否则，拒绝原假设，认为至少有一组样本均值不相等。

事后多重比较

若ANOVA显著，可进行事后多重比较，如Tukey HSD检验。

Python代码示例

import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 构造数据
data = pd.DataFrame({
    'value': np.concatenate([data1, data2, data3]),
    'group': np.repeat(['A', 'B

', 'C'], 100)
})

# ANOVA
model = ols('value ~ C(group)', data=data).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

# 事后多重比较
tukey = pairwise_tukeyhsd(endog=data['value'], groups=data['group'], alpha=0.05)
print(tukey)

在这段代码中，我们首先使用 ANOVA 检验三组样本的均值是否相等，然后进行 Tukey HSD 检验进行事后多重比较，找出具体哪组样本之间的均值存在显著差异。

1.6 相关系数

相关系数用于衡量两个变量之间的线性关系。常见的相关系数包括皮尔逊相关系数和斯皮尔曼相关系数。

皮尔逊相关系数

皮尔逊相关系数用于衡量两个连续变量之间的线性关系。其值范围为 [-1, 1]，计算公式为：
$$r=\frac{\sum (x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum (x_i-\bar{x}{)}^2\sum (y_i-\bar{y}{)}^2}}$$

Python代码示例

from scipy.stats import pearsonr

# 两组样本
x = np.random.normal(0, 1, 100)
y = x + np.random.normal(0, 0.1, 100)

corr, p = pearsonr(x, y)
print('Pearson correlation=%.3f, p=%.3f' % (corr, p))

在这段代码中，我们生成了两组正态分布的数据，并使用皮尔逊相关系数计算两组样本之间的相关系数。如果 p 值小于 0.05，表明两组样本之间存在显著的线性关系；否则，认为不存在显著的线性关系。

2. 回归不止用于预测

2.1 逻辑回归

逻辑回归是一种广泛应用于分类问题的回归分析方法。它用于预测二分类结果（例如，成功/失败，0/1等）。逻辑回归模型的输出是一个概率，表示某一事件发生的可能性。

逻辑回归模型

逻辑回归模型的基本形式是：
$$P(y=1|x)=\frac{1}{1+\exp (-({\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_n{x}_n))}$$
其中，$P(y=1|x)$ 是事件 $y=1$ 发生的概率，${\beta }_0,{\beta }_1,\dots ,{\beta }_n$ 是模型参数，$x_1,x_2,\dots ,x_n$ 是自变量。

Python代码示例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix, classification_report

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 2)
y = (X[:, 0] + X[:, 1] > 1).astype(int)

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 逻辑回归模型
model = LogisticRegression()
model.fit(X_train, y_train)

# 预测
y_pred = model.predict(X_test)

# 评价模型
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))

# 绘制决策边界
xx, yy = np.meshgrid(np.linspace(0, 1, 500), np.linspace(0, 1, 500))
Z = model.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
plt.contourf(xx, yy, Z, alpha=0.8)
plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Logistic Regression Decision Boundary')
plt.show()

在这段代码中，我们生成了一个简单的模拟数据集，并使用逻辑回归模型进行二分类预测。然后，我们使用混淆矩阵和分类报告来评估模型性能，并绘制了决策边界。

2.2 分位数回归

分位数回归是一种用于估计因变量的条件分位数的回归分析方法。与传统的最小二乘回归不同，分位数回归估计的是条件分布的特定分位数（例如，中位数、四分位数等）。

分位数回归模型

分位数回归模型的基本形式与最小二乘回归相似，但其目标是最小化不同分位数的加权绝对误差。对于分位数 $\tau$，其目标函数为：
$$\underset{\beta }{\min }\sum _{i=1}^n{\rho }_{\tau }(y_i-x_i^T\beta )$$
其中，${\rho }_{\tau }$ 是检查函数，定义为：
$${\rho }_{\tau }(u)=u(\tau -I(u<0))$$
$I(\cdot )$ 是指示函数。

Python代码示例

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# 生成模拟数据
np.random.seed(0)
n = 100
X = np.random.rand(n)
y = 2 + 3*X + np.random.randn(n)

# 添加常数项
X = sm.add_constant(X)

# 分位数回归模型
model = sm.QuantReg(y, X)

# 拟合不同分位数
quantiles = [0.1, 0.5, 0.9]
models = [model.fit(q=q) for q in quantiles]
lines = [model.params[0] + model.params[1] * X[:, 1] for model in models]

# 绘制回归结果
plt.scatter(X[:, 1], y, alpha=0.5)
for i, q in enumerate(quantiles):
    plt.plot(X[:, 1], lines[i], label=f'Quantile {q}')
plt.xlabel('Feature')
plt.ylabel('Target')
plt.legend()
plt.title('Quantile Regression')
plt.show()

3. 使用scipy与statsmodels完成统计建模

3.1 scipy进行统计模型构建

Scipy.stats模块是用于统计计算的模块。该模块提供了大量的概率分布函数，包括连续分布（如正态分布、指数分布、泊松分布等）和离散分布（如二项分布、泊松分布、超几何分布等）。这些概率分布函数可以进行概率密度函数（PDF）和累积分布函数（CDF）的计算，以及生成随机数。此外，stats模块还包含了许多常用的统计函数，如描述性统计、假设检验、非参数统计、拟合优度检验等。这些函数可以帮助用户进行各种统计分析任务，如计算均值、中位数、方差、标准差等描述性统计指标，进行t检验、卡方检验、F检验等假设检验，以及计算样本的拟合优度等。它的常用功能可以列举如下：

• 生成特定分布的概率密度。SciPy的stats模块提供了许多概率分布函数，包括连续分布和离散分布。连续分布包括：正态分布norm（norm是分布的类名，下同）；指数分布expon，对数正态分布lognorm，柯西分布cauchy，拉普拉斯分布laplace，贝塔分布beta，帕累托分布pareto，指数分布等。离散分布包括：二项分布binom，泊松分布poisson，超几何分布hypergeom，几何分布geom，离散均匀分布uniform等。除此以外，还可以实现多元正态分布multivariate_normal、高斯过程分布：gaussian_process_gp等。这些概率分布函数可以用于计算概率密度函数（PDF）、累积分布函数（CDF）、分位数函数（quantile function）等。
• 根据分布生成随机数：通过上述类的rvs方法可以生成对应的样本。如果已经给定了一个离散分布律，可以使用numpy库中的np.random.choice函数来生成样本。
• 进行一些常见的假设检验：在前面的学习中我们已经接触过一些常用的假设检验函数，包括单样本t检验ttest_1samp，配对t检验ttest_rel，独立样本t检验ttest_ind，方差齐性检验levene ，方差检验f_oneway，卡方检验chi2_contingency。使用这些函数的一般步骤是：首先，导入相应的函数；然后，准备数据并指定要进行的假设检验的类型；最后，调用该函数并传入数据，获取结果。
• 概率密度进行拟合：可以使用fit方法拟合特定分布，下面是一个拟合正态分布的例子，通过这个例子熟悉fit方法的使用。

import numpy as np  
from scipy import stats  
import matplotlib.pyplot as plt  
# 生成样本数据  
x = np.random.normal(loc=0, scale=1, size=1000)  
# 对正态分布进行拟合  
dist = stats.norm  
fit_result = dist.fit(x)  
# 打印拟合参数  
print("拟合均值：", fit_result[0])  
print("拟合标准差：", fit_result[1])  
# 绘制概率密度图  
x_fit = np.linspace(dist.ppf(0.01), dist.ppf(0.99), 100)  
y_fit = dist.pdf(x_fit, *fit_result[:2])  
plt.plot(x_fit, y_fit, 'r-', lw=5)  
# 绘制样本数据直方图  
plt.hist(x, bins=30, density=True, alpha=0.5, color='g')  
plt.show()

运行结果：

拟合均值： -0.02407312279382517
拟合标准差： 0.981194371899449

3.2 statsmodels进行统计模型构建

Statsmodels是一个Python库，用于拟合多种统计模型、执行统计测试以及进行数据探索和可视化。它提供了丰富的统计模型、统计检验和描述性统计分析的功能，可以用于回归分析、时间序列分析等多种统计任务。Statsmodels库包含线性模型、广义线性模型和鲁棒线性模型等线性模型，以及线性混合效应模型等。此外，它还包含方差分析（ANOVA）方法、时间序列过程和状态空间模型、广义的矩量法等，可以进行假设检验、置信区间计算、效应量计算等统计测试，并提供了概率分布的计算和可视化等功能。

它的常用功能可以列举如下：

• 方差分析（ANOVA）是统计学中用于检验多个组间均值差异显著性的方法，可以通过Statsmodels中的anova_lm函数来实现。此函数专门针对线性最小二乘（OLS）模型进行方差分析，使用户能够判断不同组之间的均值差异是否统计上显著，为数据分析和决策提供重要参考。

• 线性回归是数据分析中的基础技术，Statsmodels的OLS类提供了执行此类分析的功能。OLS能够处理包括简单线性回归、多元线性回归和多项式回归在内的多种回归模型，通过估计因变量与自变量之间的线性关系，帮助预测因变量的值。无论是探索自变量的影响、控制混杂变量，还是进行统计推断，OLS都是经济学、社会科学、医学和生物学等领域中不可或缺的分析工具。

• 时间序列分析是研究时间序列数据以揭示其内在规律和特征的统计方法。Statsmodels的tsa模块提供了包括ARIMA模型在内的多种时间序列分析工具，如季节性自回归积分滑动平均模型（SARIMA）、向量自回归模型（VAR）等。这些方法不仅支持时间序列的建模和预测，还能进行平稳性检验、趋势分析等，使用户能够更深入地理解并分析时间序列数据。

线性回归

首先，我们进行一个简单的线性回归分析，使用 Statsmodels 的 OLS 类。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 1)
y = 2 + 3 * X + np.random.randn(100, 1)

# 添加常数项
X = sm.add_constant(X)

# 线性回归模型
model = sm.OLS(y, X).fit()

# 打印回归结果
print(model.summary())

# 绘制回归结果
plt.scatter(X[:, 1], y, label='Data')
plt.plot(X[:, 1], model.predict(X), color='red', label='Fitted Line')
plt.xlabel('X')
plt.ylabel('y')
plt.legend()
plt.title('Linear Regression')
plt.show()

方差分析（ANOVA）

接下来，我们进行一个方差分析，使用 Statsmodels 的 anova_lm 函数。

import statsmodels.formula.api as smf

# 生成模拟数据
np.random.seed(0)
data = pd.DataFrame({
    'value': np.random.rand(30),
    'group': np.repeat(['A', 'B', 'C'], 10)
})

# 方差分析模型
model = smf.ols('value ~ C(group)', data=data).fit()
anova_results = sm.stats.anova_lm(model, typ=2)

# 打印方差分析结果
print(anova_results)

时间序列分析

最后，我们进行一个时间序列分析，使用 Statsmodels 的 ARIMA 模型。

from statsmodels.tsa.arima.model import ARIMA

# 生成模拟时间序列数据
np.random.seed(0)
data = np.cumsum(np.random.randn(100))

# 拟合 ARIMA 模型
model = ARIMA(data, order=(1, 1, 1)).fit()

# 打印模型摘要
print(model.summary())

# 绘制预测结果
forecast = model.get_forecast(steps=10)
forecast_index = np.arange(len(data), len(data) + 10)
forecast_values = forecast.predicted_mean
conf_int = forecast.conf_int()

plt.plot(data, label='Observed')
plt.plot(forecast_index, forecast_values, label='Forecast', color='red')
plt.fill_between(forecast_index, conf_int.iloc[:, 0], conf_int.iloc[:, 1], color='red', alpha=0.3)
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.title('ARIMA Model Forecast')
plt.show()

4. 机器学习简介

4.1 机器学习中的重要概念

机器学习是一种通过数据训练模型来进行预测或决策的技术。以下是一些关键概念：

1. 数据集（Dataset）：数据集是机器学习的基础，通常分为训练集、验证集和测试集。训练集用于训练模型，验证集用于调参，测试集用于评估模型性能。

2. 特征（Feature）和标签（Label）：特征是用于描述数据的属性或变量，标签是数据的目标值或分类。

3. 监督学习（Supervised Learning）：监督学习是通过已知的输入输出对进行训练，目的是预测未知数据的输出。常见算法有线性回归、逻辑回归、支持向量机、决策树、随机森林等。

4. 无监督学习（Unsupervised Learning）：无监督学习是从没有标签的数据中发现模式。常见算法有聚类（如K-means）、降维（如主成分分析PCA）等。

5. 半监督学习（Semi-Supervised Learning）：半监督学习结合了少量有标签数据和大量无标签数据来训练模型。

6. 强化学习（Reinforcement Learning）：强化学习是通过与环境交互来学习策略，以最大化累积奖励。

7. 模型评估（Model Evaluation）：常用的评估指标有准确率、精确率、召回率、F1分数、AUC-ROC等。

8. 过拟合（Overfitting）和欠拟合（Underfitting）：过拟合是指模型在训练集上表现很好，但在测试集上表现不佳。欠拟合是指模型在训练集和测试集上都表现不佳。

4.2 机器学习算法的区分

机器学习算法可以根据不同的任务和方法进行分类：

1. 回归算法（Regression Algorithms）：用于预测连续值。

   • 线性回归（Linear Regression）
   • 岭回归（Ridge Regression）
   • 套索回归（Lasso Regression）

2. 分类算法（Classification Algorithms）：用于预测离散的类别标签。

   • 逻辑回归（Logistic Regression）
   • 支持向量机（SVM）
   • 决策树（Decision Tree）
   • 随机森林（Random Forest）
   • K近邻（K-Nearest Neighbors, KNN）
   • 朴素贝叶斯（Naive Bayes）

3. 聚类算法（Clustering Algorithms）：用于将数据分组。

   • K-means聚类
   • 层次聚类（Hierarchical Clustering）
   • DBSCAN

4. 降维算法（Dimensionality Reduction Algorithms）：用于减少特征数量。

   • 主成分分析（PCA）
   • t-SNE
   • LDA

5. 强化学习算法（Reinforcement Learning Algorithms）：用于学习策略以最大化奖励。

   • Q学习（Q-Learning）
   • 深度Q网络（DQN）

Python代码示例

我们将展示一些常见的机器学习算法的示例代码。

线性回归

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 1) * 10
y = 2.5 * X + np.random.randn(100, 1) * 2

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# 线性回归模型
model = LinearRegression()
model.fit(X_train, y_train)

# 预测
y_pred = model.predict(X_test)

# 绘制结果
plt.scatter(X_test, y_test, label='True Data')
plt.plot(X_test, y_pred, color='red', label='Predictions')
plt.xlabel('X')
plt.ylabel('y')
plt.legend()
plt.title('Linear Regression')
plt.show()

逻辑回归

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 2)
y = (X[:, 0] + X[:, 1] > 1).astype(int)

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 逻辑回归模型
model = LogisticRegression()
model.fit(X_train, y_train)

# 预测
y_pred = model.predict(X_test)

# 评价模型
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))

# 绘制决策边界
xx, yy = np.meshgrid(np.linspace(0, 1, 500), np.linspace(0, 1, 500))
Z = model.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
plt.contourf(xx, yy, Z, alpha=0.8)
plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Logistic Regression Decision Boundary')
plt.show()

K-means聚类

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 2)

# K-means模型
kmeans = KMeans(n_clusters=3, random_state=0)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)

# 绘制聚类结果
plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')
centers = kmeans.cluster_centers_
plt.scatter(centers[:, 0], centers[:, 1], c='red', s=200, alpha=0.75, marker='X')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('K-means Clustering')
plt.show()

主成分分析（PCA）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 5)

# PCA模型
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# 绘制降维结果
plt.scatter(X_pca[:, 0], X_pca[:, 1])
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.title('PCA')
plt.show()

5. 基于距离的KNN模型

K最近邻（K-Nearest Neighbors，KNN）是一种简单且直观的监督学习算法，广泛应用于分类和回归问题。KNN的基本思想是：给定一个新的数据点，通过计算它与训练数据集中所有数据点之间的距离，选择距离最近的K个点，并通过这K个点的标签来预测新数据点的标签。

5.1 KNN模型的基本原理

KNN模型的基本原理可以分为以下几个步骤：

1. 选择K值：选择一个正整数K，表示在做出预测时要考虑的最近邻居的数量。
2. 计算距离：对于每个待预测数据点，计算它与训练数据集中所有数据点的距离。常用的距离度量包括欧几里得距离、曼哈顿距离和余弦相似度等。对于欧几里得距离，计算公式为：
   $$d(x_i,x_j)=\sqrt{\sum _{k=1}^n(x_{ik}-x_{jk}{)}^2}$$
   其中，$x_i$ 和 $x_j$ 是两个数据点的特征向量，$n$ 是特征的维数。
3. 寻找最近邻：根据计算出的距离，选择距离待预测数据点最近的K个数据点。
4. 预测：对于分类问题，通过K个最近邻的数据点的类别标签进行投票，以少数服从多数的方式决定待预测数据点的类别。对于回归问题，使用K个最近邻的数据点的均值作为预测值。

5.2 KNN的代码实现

以下是一个KNN分类和回归的Python代码实现示例，使用scikit-learn库。

KNN分类

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 2)
y = (X[:, 0] + X[:, 1] > 1).astype(int)

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# KNN分类模型
knn = KNeighborsClassifier(n_neighbors=3)
knn.fit(X_train, y_train)

# 预测
y_pred = knn.predict(X_test)

# 评价模型
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))

# 绘制决策边界
xx, yy = np.meshgrid(np.linspace(0, 1, 500), np.linspace(0, 1, 500))
Z = knn.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
plt.contourf(xx, yy, Z, alpha=0.8)
plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('KNN Classification Decision Boundary')
plt.show()

• KNN分类：我们生成了一组二维模拟数据，并使用KNN算法进行分类。在这段代码中，我们将数据分为训练集和测试集，使用3个最近邻进行分类，然后打印混淆矩阵和分类报告，最后绘制决策边界。
  

KNN回归

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsRegressor
import matplotlib.pyplot as plt

# 生成模拟数据
np.random.seed(0)
X = np.random.rand(100, 1) * 10
y = 2.5 * X + np.random.randn(100, 1) * 2

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# KNN回归模型
knn_reg = KNeighborsRegressor(n_neighbors=3)
knn_reg.fit(X_train, y_train)

# 预测
y_pred = knn_reg.predict(X_test)

# 绘制结果
plt.scatter(X_test, y_test, label='True Data')
plt.scatter(X_test, y_pred, color='red', label='Predictions')
plt.xlabel('X')
plt.ylabel('y')
plt.legend()
plt.title('KNN Regression')
plt.show()

• KNN回归：我们生成了一组一维模拟数据，并使用KNN算法进行回归。同样，我们将数据分为训练集和测试集，使用3个最近邻进行回归，然后绘制预测结果与真实数据的散点图。

6. 聚类算法

聚类是一种无监督学习方法，其目标是将一组数据点分成若干组，使得同一组内的数据点相似度高，而不同组间的数据点相似度低。以下是几种常见的聚类算法。

6.1 基于原型的聚类算法：KMeans算法

KMeans算法是一种迭代聚类算法，其核心思想是通过不断迭代将数据点划分到最近的聚类中心，直到聚类中心不再变化或变化很小。

KMeans算法通过迭代将数据点分配到最近的聚类中心，直到聚类中心不再变化。其核心公式如下：

1. 距离计算：
   $$d(x_i,{\mu }_j)=\sqrt{\sum _{k=1}^n(x_{ik}-{\mu }_{jk}{)}^2}$$
   其中，$x_i$ 是数据点，${\mu }_j$ 是第j个聚类中心，$n$ 是特征的维数。

2. 更新聚类中心：
   $${\mu }_j=\frac{1}{|C_j|}\sum _{x_i\in C_j}{x}_i$$
   其中，$C_j$ 是属于第j个聚类的所有数据点的集合。

步骤：

1. 随机初始化K个聚类中心。
2. 计算每个数据点到聚类中心的距离，并将数据点分配到最近的聚类中心。
3. 更新聚类中心为当前簇中所有数据点的均值。
4. 重复步骤2和步骤3，直到聚类中心不再变化或变化很小。

代码实现：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

# 生成模拟数据
np.random.seed(0)
X = np.random.randn(300, 2)

# KMeans聚类
kmeans = KMeans(n_clusters=3)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)

# 绘制聚类结果
plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')
centers = kmeans.cluster_centers_
plt.scatter(centers[:, 0], centers[:, 1], c='red', s=200, alpha=0.75, marker='X')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('KMeans Clustering')
plt.show()

6.2 基于层次的聚类算法：层次聚类法

层次聚类算法通过构建一个层次结构的树（树状图）来实现聚类，主要分为两类：自下而上的凝聚层次聚类和自上而下的分裂层次聚类。

层次聚类法按照聚类的顺序可以分为自底向上和自顶向下两种类型。

• 自底向上的层次聚类（凝聚层次聚类）：先单独考虑每个对象，然后逐渐合并这些原子簇，直到所有对象都在一个簇中或满足终止条件。这种方法的优点是不需要预先设定簇的数量，适合发现任意形状的簇，并且对异常值不敏感。缺点是计算量大，容易受到初始值的影响，可能陷入局部最优解。
• 自顶向下的层次聚类（分裂层次聚类）：开始时将所有对象视作一个簇，然后逐渐分裂这个簇，直到每个对象都是一个簇或满足终止条件。优点是可以发现任意形状的簇，能够处理大规模数据集。缺点是需要预先设定簇的数量，容易受到初始值的影响，可能陷入局部最优解。

这两种方法没有孰优孰劣之分，只是在实际应用的时候要根据数据特点以及你想要的“类”的个数，来考虑是自上而下更快还是自下而上更快。判断”类”的方法就是：最短距离法、最长距离法、中间距离法、类平均法等，其中类平均法往往被认为是最常用也最好用的方法，一方面因为其良好的单调性，另一方面因为其空间扩张/浓缩的程度适中。

6.2.1 凝聚层次聚类

分裂层次聚类中，关键步骤是如何选择簇进行分裂和如何进行分裂。这里介绍两种常见的度量：

1. 簇的直径：簇内点对之间的最大距离。
   $$\text{diam}(C)=\underset{x,y\in C}{\max }d(x,y)$$
   其中，$d(x,y)$ 是数据点$x$和$y$之间的距离。

2. 分裂方法：通常使用KMeans或其他聚类算法进行分裂。

   • 例如使用KMeans将一个簇分成两个子簇：
     $$d(x_i,{\mu }_j)=\sqrt{\sum _{k=1}^n(x_{ik}-{\mu }_{jk}{)}^2}$$
     其中，$x_i$ 是数据点，${\mu }_j$ 是第j个聚类中心，$n$ 是特征的维数。

步骤(分裂层次聚类)

1. 初始状态：将所有数据点视为一个簇。
2. 选择簇：选择一个簇进行分裂，通常选择具有最大直径或其他度量最大值的簇。
3. 分裂簇：将选择的簇分裂成两个或多个子簇，常见的方法是应用KMeans或其他聚类算法。
4. 更新簇：将分裂产生的子簇加入簇集合中。
5. 重复步骤2和步骤3，直到达到所需的簇数量或每个簇中的数据点个数达到某个阈值。

代码实现

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.datasets import make_blobs

# 生成模拟数据
np.random.seed(0)
X, _ = make_blobs(n_samples=300, centers=4, cluster_std=0.60, random_state=0)


# 分裂层次聚类函数
def divisive_clustering(X, n_clusters):
    clusters = [X]  # 初始簇
    while len(clusters) < n_clusters:
        # 选择最大簇进行分裂
        largest_cluster_index = np.argmax([len(cluster) for cluster in clusters])
        largest_cluster = clusters.pop(largest_cluster_index)

        # 使用KMeans将最大簇分裂成两个子簇
        kmeans = KMeans(n_clusters=2, random_state=0).fit(largest_cluster)
        sub_clusters = [largest_cluster[kmeans.labels_ == i] for i in range(2)]

        # 添加分裂后的子簇
        clusters.extend(sub_clusters)

    return clusters


# 执行分裂层次聚类
n_clusters = 4
clusters = divisive_clustering(X, n_clusters)

# 创建一个距离矩阵用于绘制树状图
linked = linkage(X, 'ward')

# 绘制树状图
plt.figure(figsize=(10, 7))
dendrogram(linked, truncate_mode='level', p=5)
plt.title('Divisive Hierarchical Clustering Dendrogram')
plt.xlabel('Sample index')
plt.ylabel('Distance')
plt.show()

# 绘制聚类结果
plt.figure(figsize=(10, 7))
colors = ['red', 'green', 'blue', 'purple']
for i, cluster in enumerate(clusters):
    plt.scatter(cluster[:, 0], cluster[:, 1], c=colors[i], label=f'Cluster {i + 1}')
plt.title('Divisive Hierarchical Clustering')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.legend()
plt.show()

6.2.2 分裂层次聚类

层次聚类通过构建一个层次结构的树（树状图）来实现聚类。这里介绍常见的凝聚层次聚类（自下而上）。

1. 初始状态：每个数据点作为一个单独的簇。

2. 距离计算：

   • 最小距离（single linkage）：
     $$d(C_i,C_j)=\underset{x\in C_i,y\in C_j}{\min }d(x,y)$$
   • 最大距离（complete linkage）：
     $$d(C_i,C_j)=\underset{x\in C_i,y\in C_j}{\max }d(x,y)$$
   • 平均距离（average linkage）：
     $$d(C_i,C_j)=\frac{1}{|C_i||C_j|}\sum _{x\in C_i}\sum _{y\in C_j}d(x,y)$$

3. 合并簇：选择距离最小的两个簇合并，重复此步骤直到达到所需的簇数量或合并所有簇。

步骤（凝聚层次聚类）：

1. 将每个数据点视为一个簇。
2. 计算所有簇之间的距离，并合并距离最小的两个簇。
3. 重复步骤2，直到所有数据点合并为一个簇。

代码实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.datasets import make_blobs

# 生成模拟数据
np.random.seed(0)
X, _ = make_blobs(n_samples=300, centers=4, cluster_std=0.60, random_state=0)

# 层次聚类
linked = linkage(X, 'ward')

# 绘制树状图
plt.figure(figsize=(10, 7))
dendrogram(linked, truncate_mode='level', p=5)
plt.title('Hierarchical Clustering Dendrogram')
plt.xlabel('Sample index')
plt.ylabel('Distance')
plt.show()

6.3 基于密度的聚类算法：DBSCAN算法

DBSCAN（Density-Based Spatial Clustering of Applications with Noise）是一种基于密度的聚类算法，能够发现任意形状的簇，并识别噪声点。

DBSCAN通过密度的方式发现簇，其核心公式和概念如下：

1. ε邻域：对于一个点$p$，ε邻域包含所有与$p$的距离小于ε的点：
   N ϵ ( p ) = { q ∈ D ∣ d ( p , q ) ≤ ϵ } N_\epsilon(p) = \{ q \in D | d(p, q) \leq \epsilon \} Nϵ​(p)={q∈D∣d(p,q)≤ϵ}

2. 核心点：如果一个点的ε邻域包含至少MinPts个点，则该点为核心点：
   $$|N_{ϵ}(p)|\ge \text{MinPts}$$

3. 直接密度可达：如果点$q$在点$p$的ε邻域内，则$q$是直接密度可达的：
   $$q\in N_{ϵ}(p)$$

4. 密度可达：如果存在一系列点$p_1,p_2,\dots ,p_n$，使得$p_{i+1}$是直接密度可达于$p_i$，则$p_n$是密度可达的。

步骤：

1. 对每个点，计算其ε邻域内的点数。
2. 如果某个点的ε邻域内的点数大于等于MinPts，则将其标记为核心点。
3. 对于每个核心点，将其ε邻域内的所有点（包括核心点和边界点）归为一个簇。
4. 重复步骤2和步骤3，直到所有核心点都被处理。

代码实现：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_blobs

# 生成模拟数据
np.random.seed(0)
X, _ = make_blobs(n_samples=300, centers=4, cluster_std=0.60, random_state=0)

# DBSCAN聚类
# eps参数指定了两个样本被认为是邻居的最大距离，而min_samples参数指定了一个样本点被视为核心点所需的最小邻居数。通过调整这些参数，我们可以控制聚类的结果。
dbscan = DBSCAN(eps=0.5, min_samples=5)
y_dbscan = dbscan.fit_predict(X)

# 绘制聚类结果
plt.scatter(X[:, 0], X[:, 1], c=y_dbscan, s=50, cmap='viridis')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('DBSCAN Clustering')
plt.show()

6.4 基于模型的聚类算法：高斯混合聚类法

高斯混合模型（Gaussian Mixture Model，GMM）是一种基于概率模型的聚类算法，假设数据是由多个高斯分布组成。

高斯混合模型假设数据是由多个高斯分布组成，其核心公式包括：

1. 高斯分布概率密度函数：
   $$\mathcal{N}(x|\mu ,\Sigma )=\frac{1}{(2\pi {)}^{d/2}|\Sigma {|}^{1/2}}\exp \left(-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )\right)$$
   其中，$\mu$ 是均值向量，$\Sigma$ 是协方差矩阵。

2. 混合高斯模型：
   $$p(x|\lambda )=\sum _{j=1}^K{\pi }_j\mathcal{N}(x|{\mu }_j,{\Sigma }_j)$$
   其中，${\pi }_j$ 是第j个高斯分布的混合系数，满足${\sum }_{j=1}^K{\pi }_j=1$。

3. E步（Expectation）：计算每个数据点属于每个高斯分布的后验概率（即责任度）：
   $$\gamma (z_{ij})=\frac{{\pi }_j\mathcal{N}(x_i|{\mu }_j,{\Sigma }_j)}{{\sum }_{k=1}^K{\pi }_k\mathcal{N}(x_i|{\mu }_k,{\Sigma }_k)}$$

4. M步（Maximization）：更新参数，使对数似然函数最大化：
   $${\pi }_j=\frac{1}{N}\sum _{i=1}^N\gamma (z_{ij})$$
   $${\mu }_j=\frac{{\sum }_{i=1}^N\gamma (z_{ij})x_i}{{\sum }_{i=1}^N\gamma (z_{ij})}$$
   $${\Sigma }_j=\frac{{\sum }_{i=1}^N\gamma (z_{ij})(x_i-{\mu }_j)(x_i-{\mu }_j{)}^T}{{\sum }_{i=1}^N\gamma (z_{ij})}$$

步骤：

1. 初始化每个高斯分布的参数（均值、协方差矩阵和混合系数）。
2. E步：计算每个数据点属于每个高斯分布的后验概率。
3. M步：更新高斯分布的参数，使得对数似然函数最大化。
4. 重复E步和M步，直到参数收敛。

代码实现：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture

# 生成模拟数据
np.random.seed(0)
X, _ = make_blobs(n_samples=300, centers=4, cluster_std=0.60, random_state=0)

# 高斯混合模型聚类
gmm = GaussianMixture(n_components=4)
gmm.fit(X)
y_gmm = gmm.predict(X)

# 绘制聚类结果
plt.scatter(X[:, 0], X[:, 1], c=y_gmm, s=50, cmap='viridis')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Gaussian Mixture Model Clustering')
plt.show()

6.5 聚类算法总结

通过以上示例，我们展示了几种常见聚类算法的原理和实现：

• KMeans算法：通过迭代寻找聚类中心，将数据点分配到最近的聚类中心。
• 层次聚类法：通过构建树状图，自下而上或自上而下地合并或分裂簇。
• DBSCAN算法：通过密度的方式发现簇，可以识别噪声点。
• 高斯混合聚类法：通过概率模型，假设数据是由多个高斯分布组成。

每种聚类算法都有其优缺点和适用场景，可以根据具体的数据特点和任务需求选择合适的算法。

一些疑问

高斯混合聚类为什么要更新参数，使对数似然函数最大化

对数似然函数

对数似然函数是指模型参数 θ = { μ k , Σ k , π k } \theta = \{\mu_k, \Sigma_k, \pi_k\} θ={μk​,Σk​,πk​} 下，观测数据的对数概率。对数似然函数的形式为：

$$\log L(\theta )=\sum _{i=1}^N\log \left(\sum _{k=1}^K{\pi }_k\mathcal{N}(x_i|{\mu }_k,{\Sigma }_k)\right)$$

其中，$N$ 是数据点的总数，$K$ 是高斯分布的个数，$\mathcal{N}(x_i|{\mu }_k,{\Sigma }_k)$ 表示第 $k$ 个高斯分布对数据点 $x_i$ 的概率密度函数。

对数似然函数的目的是衡量在给定参数 $\theta$ 下，观测数据出现的概率。我们希望通过调整参数 $\theta$ 来最大化这个概率，即使得对数似然函数最大化。

EM算法步骤

1. 初始化：随机初始化参数 θ = { μ k , Σ k , π k } \theta = \{\mu_k, \Sigma_k, \pi_k\} θ={μk​,Σk​,πk​}。

2. E步（期望步）：计算每个数据点属于每个高斯分布的责任度（posterior probability），即 ${\gamma }_{ik}$：

$${\gamma }_{ik}=\frac{{\pi }_k\mathcal{N}(x_i|{\mu }_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_j\mathcal{N}(x_i|{\mu }_j,{\Sigma }_j)}$$

3. M步（最大化步）：根据责任度更新参数 $\theta$：

$${\mu }_k=\frac{{\sum }_{i=1}^N{\gamma }_{ik}{x}_i}{{\sum }_{i=1}^N{\gamma }_{ik}}$$

$${\Sigma }_k=\frac{{\sum }_{i=1}^N{\gamma }_{ik}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T}{{\sum }_{i=1}^N{\gamma }_{ik}}$$

$${\pi }_k=\frac{{\sum }_{i=1}^N{\gamma }_{ik}}{N}$$

4. 重复：重复E步和M步直到对数似然函数收敛。

对数似然函数的作用

对数似然函数在GMM和其他统计模型中具有重要作用：

1. 参数估计：通过最大化对数似然函数，可以得到最优的模型参数，使得观测数据的概率最大化。
2. 模型评估：对数似然函数的值可以用于比较不同模型的拟合优度。对数似然值越大，模型对数据的拟合越好。
3. 优化目标：在EM算法中，对数似然函数提供了明确的优化目标，通过迭代更新参数使其逐步收敛到最优值。

学习心得

通过这次学习了解了正态性检验和独立性检验的检验方法。
当自变量是离散变量而因变量是连续变量的时候可以通过方差分析来分析自变量操作的不同对因变量分组均值水平的影响。反过来，自变量是连续变量而因变量是离散变量的时候如果做线性回归则容易出现斜率暴增的现象，引入了逻辑回归。分位数回归可以提供更加全面和准确的预测，并且能够更好地处理异常值和复杂的非线性关系。掌握了使用scipy和statsmodels两个库对数据进行统计处理。
介绍了机器学习的基本概念，如：数据集、特征、训练、验证、测试、偏差、方差等。
并且罗列了机器学习算法的不同分类：回归算法、分类算法、聚类、降维、强化学习。
详细介绍了KNN、XGBoost、K-Mean、神经网络、贝叶斯算法。
其中K最近邻（K-Nearest Neighbors，KNN）是一种简单且直观的监督学习算法，广泛应用于分类和回归问题。KNN的基本思想是：给定一个新的数据点，通过计算它与训练数据集中所有数据点之间的距离，选择距离最近的K个点，并通过这K个点的标签来预测新数据点的标签。
集成学习是一种机器学习方法，通过结合多个学习器的预测结果来提高整体的预测精度。XGBoost对梯度进行了优化，提高了预测的准确度。与传统的机器学习算法相比，XGBoost通过优化梯度下降的方式，更加精确地逼近真实函数，从而在许多数据集上取得了更高的预测精度。其次，XGBoost能够自主处理数据缺失的问题，这极大地削弱了我们预处理的困难
k-means 算法是根据给定的 n 个数据对象的数据集，构建 k 个划分聚类的方法，每个划分聚类即为一个簇。该方法将数据划分为 n 个簇，每个簇至少有一个数据对象，每个数据对象必须属于而且只能属于一个簇。同时要满足同一簇中的数据对象相似度高，不同簇中的数据对象相似度较小。