机器学习实战|采用聚类和回归探究蓝莓产量情况

最新推荐文章于 2024-06-07 19:59:11 发布

暴躁的秋秋

最新推荐文章于 2024-06-07 19:59:11 发布

阅读量1.7k

点赞数 21

文章标签：机器学习人工智能

本文链接：https://blog.csdn.net/m0_67431719/article/details/135565777

版权

1.项目背景

蓝莓在全球范围内备受欢迎，其独特的风味和丰富的营养价值令消费者为之倾倒。蓝莓生长对适宜气候和土壤的依赖，该数据集包括蓝莓的生长情况、花期内的温度以及蜜蜂的密度（表示授粉），本项目通过三个模型对该数据进行了分析，得到适合蓝莓的生长环境，以及预测了蓝莓的产量，该项目能够帮助读者更好的理解数据分析的流程，以及建立一些常见的预测模型和聚类模型。

2.数据说明

字段名描述
id 蓝莓唯一标识
Clonesize 蓝莓克隆平均大小,单位：m²
Honeybee 蜜蜂密度（单位：蜜蜂/m²/分钟）
Bumbles 大型蜜蜂密度（单位：大型蜜蜂/m²/分钟）
Andrena 安德烈纳蜂密度（单位：安德烈纳蜂/m²/分钟）
Osmia 钥匙蜂密度（单位：钥匙蜂/m²/分钟）
MaxOfUpperTRange 花期内最高温带日平均气温的最高记录, 单位: °C
MinOfUpperTRange 花期内最高温带日平均气温的最低记录, 单位: °C
AverageOfUpperTRange 花期内最高温带日平均气温, 单位: °C
MaxOfLowerTRange 花期内最低温带日平均气温的最高记录, 单位: °C
MinOfLowerTRange 花期内最低温带日平均气温的最低记录, 单位: °C
AverageOfLowerTRange 花期内最低温带日平均气温, 单位: °C
RainingDays 花期内降雨量大于 0 的日数总和, 单位: 天
AverageRainingDays 花期内降雨日数的平均值, 单位: 天
fruitset 果实集
fruitmass 果实质量
seeds 种子数
yield 产量

3.Python库导入及数据读取

In [1]:
# 导入需要的库
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import scipy.stats as stats
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import RandomizedSearchCV
In [2]:
# 读取数据
train_data = pd.read_csv("/home/mw/input/data5876/train.csv")
test_data = pd.read_csv("/home/mw/input/data5876/test.csv")
4.数据预览及数据处理

4.1数据预览

In [3]:
# 查看数据维度
print('训练集纬度：')
print(train_data.shape)
print('-'*50)
print('测试集纬度:')
print(test_data.shape)
训练集纬度：
(15289, 18)
--------------------------------------------------
测试集纬度:
(10194, 17)
In [4]:
# 查看数据信息
print('训练集信息：')
print(train_data.info())
print('-'*50)
print('测试集信息:')
print(test_data.info())
训练集信息：
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 15289 entries, 0 to 15288
Data columns (total 18 columns):
id                      15289 non-null int64
clonesize               15289 non-null float64
honeybee                15289 non-null float64
bumbles                 15289 non-null float64
andrena                 15289 non-null float64
osmia                   15289 non-null float64
MaxOfUpperTRange        15289 non-null float64
MinOfUpperTRange        15289 non-null float64
AverageOfUpperTRange    15289 non-null float64
MaxOfLowerTRange        15289 non-null float64
MinOfLowerTRange        15289 non-null float64
AverageOfLowerTRange    15289 non-null float64
RainingDays             15289 non-null float64
AverageRainingDays      15289 non-null float64
fruitset                15289 non-null float64
fruitmass               15289 non-null float64
seeds                   15289 non-null float64
yield                   15289 non-null float64
dtypes: float64(17), int64(1)
memory usage: 2.1 MB
None
--------------------------------------------------
测试集信息:
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10194 entries, 0 to 10193
Data columns (total 17 columns):
id                      10194 non-null int64
clonesize               10194 non-null float64
honeybee                10194 non-null float64
bumbles                 10194 non-null float64
andrena                 10194 non-null float64
osmia                   10194 non-null float64
MaxOfUpperTRange        10194 non-null float64
MinOfUpperTRange        10194 non-null float64
AverageOfUpperTRange    10194 non-null float64
MaxOfLowerTRange        10194 non-null float64
MinOfLowerTRange        10194 non-null float64
AverageOfLowerTRange    10194 non-null float64
RainingDays             10194 non-null float64
AverageRainingDays      10194 non-null float64
fruitset                10194 non-null float64
fruitmass               10194 non-null float64
seeds                   10194 non-null float64
dtypes: float64(16), int64(1)
memory usage: 1.3 MB
None
In [5]:
# 查看各列缺失值
print('训练集信息缺失情况：')
print(train_data.isna().sum())
print('-'*50)
print('测试集信息缺失情况:')
print(test_data.isna().sum())
训练集信息缺失情况：
id                      0
clonesize               0
honeybee                0
bumbles                 0
andrena                 0
osmia                   0
MaxOfUpperTRange        0
MinOfUpperTRange        0
AverageOfUpperTRange    0
MaxOfLowerTRange        0
MinOfLowerTRange        0
AverageOfLowerTRange    0
RainingDays             0
AverageRainingDays      0
fruitset                0
fruitmass               0
seeds                   0
yield                   0
dtype: int64
--------------------------------------------------
测试集信息缺失情况:
id                      0
clonesize               0
honeybee                0
bumbles                 0
andrena                 0
osmia                   0
MaxOfUpperTRange        0
MinOfUpperTRange        0
AverageOfUpperTRange    0
MaxOfLowerTRange        0
MinOfLowerTRange        0
AverageOfLowerTRange    0
RainingDays             0
AverageRainingDays      0
fruitset                0
fruitmass               0
seeds                   0
dtype: int64
In [6]:
# 查看重复值
print('训练集信息重复情况：')
print(train_data.duplicated().sum())
print('-'*50)
print('测试集信息重复情况:')
print(test_data.duplicated().sum())
训练集信息重复情况：
0
--------------------------------------------------
测试集信息重复情况:
0
4.2数据处理

In [7]:
# 删除ID列
train_data.drop(['id'],axis=1,inplace=True)
test_data.drop(['id'],axis=1,inplace=True)
5.相关性分析

In [8]:
# 计算相关系数矩阵
correlation_matrix = train_data.corr()
# 绘制热图来可视化相关性
plt.figure(figsize=(20,15))
sns.heatmap(correlation_matrix, annot=True, cmap='Blues', fmt=".2f")
plt.title("Correlation Matrix of Variables")
plt.show()
1.与产量（yield）高度相关的变量：fruitset、seeds 和 fruitmass 与产量高度正相关。这表明果实集、种子数和果实质量在蓝莓产量的决定中扮演重要角色。特别是，果实集和种子数与产量的相关系数非常高，表明这些因素对产量影响显著。
2.蜜蜂密度变量的相关性：蜜蜂密度相关的变量（如 honeybee、bumbles、andrena、osmia）与产量的相关性相对较低。这可能意味着虽然这些因素对蓝莓的授粉有一定作用，但它们对最终产量的直接影响可能不如果实集和种子数那么显著。
3.气温和降雨量的影响：气温和降雨量相关的变量（如 MaxOfUpperTRange、MinOfUpperTRange、RainingDays 等）与产量的相关性也不高。这可能表明在该数据集中，气温和降雨对蓝莓产量的影响不如果实的生物学特性那么直接或显著。
4.可能的多重共线性问题：fruitset、seeds 和 fruitmass 之间存在高度的相关性。这种强相关性可能导致多重共线性问题。

6.聚类分析

6.1数据处理

In [9]:
# 选择所有变量进行聚类
x_cluster = train_data.copy()
# 对数据进行标准化
scaler = StandardScaler()
x_scaled = scaler.fit_transform(x_cluster)
6.2确定聚类数

In [10]:
# 使用肘部法则来确定最佳聚类数
inertia = []
silhouette_scores = []
k_range = range(2, 11)
for k in k_range:
    kmeans = KMeans(n_clusters=k, random_state=10).fit(x_scaled)
    inertia.append(kmeans.inertia_)
    silhouette_scores.append(silhouette_score(x_scaled, kmeans.labels_))
In [11]:
plt.figure(figsize=(15,5))

plt.subplot(1, 2, 1)
plt.plot(k_range, inertia, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Inertia')
plt.title('Elbow Method For Optimal k')

plt.subplot(1, 2, 2)
plt.plot(k_range, silhouette_scores, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Silhouette Score')
plt.title('Silhouette Score For Each k')

plt.tight_layout()
plt.show()
图一为肘部法则图，当聚类数从2增加到3时，Inertia快速下降；当聚类数继续增加时，Inertia的下降速率减缓，尤其是在3或4个聚类之后。这表明3或4可能是一个合适的聚类数。
图二为轮廓系数图，最高的轮廓系数出现在2个聚类的时候，但随着聚类数的增加，轮廓系数迅速下降，并在4个聚类后稳定下来。
结合两个图，可以得知3或4个聚类可能是比较合适的聚类数目。虽然2个聚类有最高的轮廓系数，但从肘部法则图来看，3或4个聚类能够提供更加细分的聚类，且Inertia下降率变缓，意味着增加更多聚类数目带来的信息增益较小，所以我这里选择了4类。

6.3建立k均值聚类模型

In [12]:
# 执行K-均值聚类，选择4个聚类
kmeans_final = KMeans(n_clusters=4, random_state=15)
kmeans_final.fit(x_scaled)
# 获取聚类标签
cluster_labels = kmeans_final.labels_
# 将聚类标签添加到原始数据中以进行分析
train_data_clustered = train_data.copy()
train_data_clustered['Cluster'] = cluster_labels
# 查看每个聚类的统计数据
cluster_summary = train_data_clustered.groupby('Cluster').mean()
cluster_summary
Out[12]:

clonesize honeybee bumbles andrena osmia MaxOfUpperTRange MinOfUpperTRange AverageOfUpperTRange MaxOfLowerTRange MinOfLowerTRange AverageOfLowerTRange RainingDays AverageRainingDays fruitset fruitmass seeds yield
Cluster
0 17.328080 0.346109 0.296424 0.498191 0.612562 74.767492 45.189268 62.475436 53.893253 26.075593 44.225397 11.916085 0.228554 0.555068 0.469118 38.871086 6961.936706
1 22.270245 0.437966 0.275738 0.484365 0.560719 90.494726 54.712893 75.610762 65.240384 31.567928 53.465477 26.255506 0.432669 0.436473 0.418616 32.748693 4863.325850
2 17.650084 0.351646 0.295031 0.498170 0.612210 89.643195 54.205118 74.910835 64.626788 31.272423 52.962712 12.246797 0.232205 0.559047 0.477165 39.367173 6975.032587
3 22.004267 0.429347 0.278173 0.488991 0.580457 72.474424 43.797838 60.539772 52.218265 25.271977 42.855192 25.545684 0.422139 0.449389 0.415332 33.048957 5113.166299

聚类0：代表了平均克隆大小最小，蜜蜂密度最低，但有最高的钥匙蜂密度，气温范围适中，降雨天数最少，果实集、果实质量和种子数都很高，产量最高。
聚类1：特征为克隆大小和蜜蜂密度较高，气温范围最高，降雨天数最多，果实集、果实质量和种子数相对最低，产量最低。
聚类2：拥有较小的克隆大小，蜜蜂密度相对较低，气温范围高，降雨天数较少，果实集、果实质量和种子数最高，产量也非常高。
聚类3：具有较大的克隆大小和较高的蜜蜂密度，气温范围和降雨天数适中，果实集和种子数介于中等至较高，果实质量相对较低，产量中等。

通过对聚类特征的分析，我们可以推测可能的种植策略和环境条件：
聚类0和聚类2具有相似的果实特性和产量，但聚类0的气温范围较低，降雨天数少，表明在较干燥和凉爽的环境中，通过优化克隆大小和蜜蜂密度可能能获得高产量。
聚类1在所有聚类中气温范围最高且降雨最多，但产量最低，这可能说明过高的温度和过多的降雨对蓝莓的生产不利。
聚类3有着适中的气温和降雨，较大的克隆大小和高蜜蜂密度，其产量和果实质量表明这可能是一个相对平衡的种植环境。

7.多元线性回归模型

7.1检测多重共线性

In [13]:
x = train_data.drop(['yield'], axis=1)  # 使用除了产量以外的所有列作为特征
x = sm.add_constant(x)
# 计算每个特征的VIF值
vif_data = pd.DataFrame()
vif_data["feature"] = x.columns
vif_data["VIF"] = [variance_inflation_factor(x.values, i) for i in range(x.shape[1])]
vif_data
/opt/conda/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)
Out[13]:

feature VIF
0 const 452.702544
1 clonesize 1.376754
2 honeybee 1.120483
3 bumbles 1.151751
4 andrena 1.186737
5 osmia 1.200157
6 MaxOfUpperTRange 3414.130675
7 MinOfUpperTRange 559.970976
8 AverageOfUpperTRange 30443.904985
9 MaxOfLowerTRange 1105.083282
10 MinOfLowerTRange 24965.147210
11 AverageOfLowerTRange 23062.819227
12 RainingDays 55.514009
13 AverageRainingDays 56.054862
14 fruitset 11.967936
15 fruitmass 12.726585
16 seeds 9.797380

高VIF值通常意味着相关变量之间存在强烈的线性关系，这可能会影响线性回归模型的准确性和稳定性。在这种情况下，可以采取主成分分析来处理多重共线性。

7.2主成分分析

In [14]:
# 数据标准化
scaler = StandardScaler()
x_scaled = scaler.fit_transform(x.drop('const', axis=1))  # 去除常数项
# 执行PCA
pca = PCA()
x_pca = pca.fit_transform(x_scaled)
# 计算主成分的方差贡献率
explained_variance = pca.explained_variance_ratio_
print('方差贡献率:')
print(explained_variance)
方差贡献率:
[3.77418076e-01 2.44286924e-01 8.51790842e-02 8.34677450e-02
 7.26665536e-02 5.73632445e-02 3.90780709e-02 3.21118075e-02
 4.43381609e-03 3.24392256e-03 5.63408861e-04 1.27603202e-04
 4.17748878e-05 1.46630704e-05 1.89651779e-06 1.40930816e-06]
In [15]:
# 计算累积方差贡献率
cumulative_variance = np.cumsum(explained_variance)
# 确定累积方差贡献率达到95%的主成分数量
n_components = np.where(cumulative_variance >= 0.95)[0][0] + 1
print(f'前{n_components}个特征已经达到了95%的累计方差贡献率。')
前7个特征已经达到了95%的累计方差贡献率。
因为前7项的特征已经达到了95%以上的贡献率，所以选择前7项进行建模。

7.3建立多元线性回归模型

In [16]:
# 使用前7个主成分作为特征集
x_pca_reduced = x_pca[:, :n_components]
# 由于PCA是无监督的，我们需要重新获取目标变量'y'
y = train_data['yield']
# 分割数据集
x_train_pca, x_test_pca, y_train, y_test = train_test_split(x_pca_reduced, y, test_size=0.3, random_state=15) # 37分
In [17]:
# 创建多元线性回归模型
regression_model = LinearRegression()
# 拟合训练数据
regression_model.fit(x_train_pca, y_train)
Out[17]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)
In [18]:
# 预测测试数据
y_pred = regression_model.predict(x_test_pca)
# 计算模型性能指标
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"Mean Squared Error (MSE): {mse}")
print(f"R-squared (R2): {r2}")
Mean Squared Error (MSE): 393493.04628487333
R-squared (R2): 0.7870589957164607
7.4残差项检验

In [19]:
# 计算残差
residuals = y_test - y_pred

# 绘制残差序列图
plt.figure(figsize=(10, 6))
plt.plot(residuals, marker='o', linestyle='')
plt.title("Residuals Time Series Plot")
plt.xlabel("Observation")
plt.ylabel("Residuals")
plt.axhline(y=0, color='r', linestyle='-')
plt.show()
残差在零线周围波动，没有显示出明显的趋势或周期性模式，看起来是随机分布的，可以认为残差项具有独立性，这意味着模型的每个预测值的误差是独立的，没有显示出依赖于其他预测值的误差，符合模型的假设，随机误差项是独立的。

In [20]:
# 绘制残差的直方图
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.hist(residuals, bins=30, edgecolor='black')
plt.title("Histogram of Residuals")
plt.xlabel("Residuals")
plt.ylabel("Frequency")

# 绘制残差的Q-Q图
plt.subplot(1, 2, 2)
stats.probplot(residuals, dist="norm", plot=plt)
plt.title("Normal Q-Q Plot of Residuals")

plt.tight_layout()
plt.show()
直方图：残差的分布看起来接近于正态分布，但似乎在中心附近稍微偏离，这可能表明残差分布略微偏斜。
Q-Q图：大多数点似乎沿着直线排列，这表明残差的分布接近正态分布。然而，图的两端有一些点偏离直线，这可能表明残差在尾部有一些偏离正态分布。
综合来看，残差大致呈正态分布，但存在轻微的偏差。这种偏差在实际应用中是常见的，尤其是当样本量较大时。尽管存在一些偏离，模型的残差还是在很大程度上符合正态分布的假设。

In [21]:
# 绘制残差与预测值的散点图以检查同方差性
plt.figure(figsize=(10, 6))
plt.scatter(y_pred, residuals, alpha=0.5)
plt.title("Residuals vs Predicted Values")
plt.xlabel("Predicted Values")
plt.ylabel("Residuals")
plt.axhline(y=0, color='r', linestyle='-')
plt.show()
在理想情况下，如果满足同方差性假设，残差应该随机分布在水平线（红色线）周围，没有任何明显的模式。例如，残差不应该随着预测值的增加而系统地增大或减小，也不应该呈现出漏斗形状。
从图中可以看出来：残差在不同的预测值水平上是随机分布的，没有明显的趋势，也没有漏斗状的分布模式，可以认为残差满足同方差性的假设。

7.5模型预测

In [22]:
# 使用与训练模型时相同数量的主成分进行预测
x_test_scaled = scaler.transform(test_data)
x_test_pca = pca.transform(x_test_scaled)
x_test_pca_reduced = x_test_pca[:, :n_components]
# 使用模型进行预测
test_predictions = regression_model.predict(x_test_pca_reduced)
test_data_with_predictions = test_data.copy()
test_data_with_predictions['predicted_yield'] = test_predictions
test_data_with_predictions.head()
Out[22]:

clonesize honeybee bumbles andrena osmia MaxOfUpperTRange MinOfUpperTRange AverageOfUpperTRange MaxOfLowerTRange MinOfLowerTRange AverageOfLowerTRange RainingDays AverageRainingDays fruitset fruitmass seeds predicted_yield
0 25.0 0.25 0.25 0.25 0.25 86.0 52.0 71.9 62.0 30.0 50.8 24.0 0.39 0.399367 0.408088 31.394569 4403.677514
1 12.5 0.25 0.25 0.75 0.63 94.6 57.2 79.0 68.2 33.0 55.9 1.0 0.10 0.488048 0.442866 36.846956 6141.142129
2 12.5 0.25 0.25 0.63 0.63 86.0 52.0 71.9 62.0 30.0 50.8 16.0 0.26 0.583379 0.487057 40.037644 7279.014537
3 25.0 0.50 0.38 0.38 0.63 86.0 52.0 71.9 62.0 30.0 50.8 16.0 0.26 0.433014 0.422847 33.116091 5133.609279
4 37.5 0.75 0.25 0.25 0.25 94.6 57.2 79.0 68.2 33.0 55.9 24.0 0.39 0.360996 0.388860 29.558019 3598.218257

在建立多元线性回归模型中，主要针对多重共线性的特征采用主成分分析降维，再对训练集进行建模，得到了模型的均方误差和决定系数，并且对模型进行了检验，发现符合多元线性回归模型的假设，即：随机误差项独立，同方差，符合正态分布，最后通过该模型对测试集进行预测，得到预测结果。

8.随机森林

8.1建立模型

In [23]:
# 划分数据
x = train_data.drop('yield', axis=1)
y = train_data['yield']
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=10, stratify=y) #37分
In [24]:
# 模型建立
rf_clf = RandomForestRegressor(random_state=15)
rf_clf.fit(x_train, y_train)
/opt/conda/lib/python3.6/site-packages/sklearn/ensemble/forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)
Out[24]:
RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
                      max_features='auto', max_leaf_nodes=None,
                      min_impurity_decrease=0.0, min_impurity_split=None,
                      min_samples_leaf=1, min_samples_split=2,
                      min_weight_fraction_leaf=0.0, n_estimators=10,
                      n_jobs=None, oob_score=False, random_state=15, verbose=0,
                      warm_start=False)
In [25]:
y_pred = rf_clf.predict(x_test)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"Mean Squared Error (MSE): {mse}")
print(f"R-squared (R2): {r2}")
Mean Squared Error (MSE): 398982.342838073
R-squared (R2): 0.776309763673141
8.2参数优化

In [26]:
# 定义随机搜索的参数范围
param_dist = {
    'n_estimators': [100, 200, 300, 400, 500],
    'max_depth': [10, 20, 30, 40, 50, None],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4],
    'max_features': ['auto', 'sqrt']
}

# 创建随机森林模型
rf = RandomForestRegressor(random_state=15)

# 设置随机搜索
random_search = RandomizedSearchCV(estimator=rf, param_distributions=param_dist, 
                                   n_iter=10, cv=5, verbose=2, random_state=17, n_jobs=-1)

# 执行随机搜索
random_search.fit(x_train, y_train)
Fitting 5 folds for each of 10 candidates, totalling 50 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=-1)]: Done  37 tasks      | elapsed:  2.6min
[Parallel(n_jobs=-1)]: Done  50 out of  50 | elapsed:  3.2min finished
Out[26]:
RandomizedSearchCV(cv=5, error_score='raise-deprecating',
                   estimator=RandomForestRegressor(bootstrap=True,
                                                   criterion='mse',
                                                   max_depth=None,
                                                   max_features='auto',
                                                   max_leaf_nodes=None,
                                                   min_impurity_decrease=0.0,
                                                   min_impurity_split=None,
                                                   min_samples_leaf=1,
                                                   min_samples_split=2,
                                                   min_weight_fraction_leaf=0.0,
                                                   n_estimators='warn',
                                                   n_jobs=None, oob_score=False,
                                                   random_state=15, verbose=0,
                                                   warm_start=False),
                   iid='warn', n_iter=10, n_jobs=-1,
                   param_distributions={'max_depth': [10, 20, 30, 40, 50, None],
                                        'max_features': ['auto', 'sqrt'],
                                        'min_samples_leaf': [1, 2, 4],
                                        'min_samples_split': [2, 5, 10],
                                        'n_estimators': [100, 200, 300, 400,
                                                         500]},
                   pre_dispatch='2*n_jobs', random_state=17, refit=True,
                   return_train_score=False, scoring=None, verbose=2)
随机搜索比网格搜索快，所以这里使用随机搜索。

In [27]:
# 最佳参数和评分
best_params = random_search.best_params_
print("最佳参数:", best_params)
最佳参数: {'n_estimators': 200, 'min_samples_split': 5, 'min_samples_leaf': 1, 'max_features': 'auto', 'max_depth': 10}
In [28]:
# 使用最佳参数创建随机森林模型
best_rf_model = RandomForestRegressor(**best_params,random_state=15)
best_rf_model.fit(x_train, y_train)
y_pred = best_rf_model.predict(x_test)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f"Mean Squared Error (MSE): {mse}")
print(f"R-squared (R2): {r2}")
Mean Squared Error (MSE): 339064.7802002252
R-squared (R2): 0.8099026631765387
In [29]:
# 获取特征重要性
feature_importances = best_rf_model.feature_importances_
# 创建特征重要性的DataFrame
features = x_train.columns
importances_df = pd.DataFrame({'Feature': features, 'Importance': feature_importances})
# 按重要性排序
importances_df.sort_values(by='Importance', ascending=False, inplace=True)
importances_df
Out[29]:

Feature Importance
13 fruitset 0.804402
15 seeds 0.149529
14 fruitmass 0.019532
4 osmia 0.003723
3 andrena 0.003491
12 AverageRainingDays 0.003144
11 RainingDays 0.002704
2 bumbles 0.002239
1 honeybee 0.001963
10 AverageOfLowerTRange 0.001536
0 clonesize 0.001513
6 MinOfUpperTRange 0.001435
8 MaxOfLowerTRange 0.001392
9 MinOfLowerTRange 0.001185
5 MaxOfUpperTRange 0.001141
7 AverageOfUpperTRange 0.001072

这表明果实集（fruitset）是最重要的特征，可能是因为：果实集是产量形成的直接因素。如果果实集比率高，意味着更多的花朵被成功受粉并开始形成果实，最终导致更高的产量，与相关性分析得到的结果一致，果实集对产量具有高强度正相关。

8.3模型预测

In [30]:
# 使用训练好的模型对测试数据进行预测
test_predictions_rf = best_rf_model.predict(test_data)
rf_test_data_with_predictions = test_data.copy()
rf_test_data_with_predictions['predicted_yield'] = test_predictions_rf
rf_test_data_with_predictions.head()
Out[30]:

clonesize honeybee bumbles andrena osmia MaxOfUpperTRange MinOfUpperTRange AverageOfUpperTRange MaxOfLowerTRange MinOfLowerTRange AverageOfLowerTRange RainingDays AverageRainingDays fruitset fruitmass seeds predicted_yield
0 25.0 0.25 0.25 0.25 0.25 86.0 52.0 71.9 62.0 30.0 50.8 24.0 0.39 0.399367 0.408088 31.394569 4275.553904
1 12.5 0.25 0.25 0.75 0.63 94.6 57.2 79.0 68.2 33.0 55.9 1.0 0.10 0.488048 0.442866 36.846956 6026.631688
2 12.5 0.25 0.25 0.63 0.63 86.0 52.0 71.9 62.0 30.0 50.8 16.0 0.26 0.583379 0.487057 40.037644 7210.325302
3 25.0 0.50 0.38 0.38 0.63 86.0 52.0 71.9 62.0 30.0 50.8 16.0 0.26 0.433014 0.422847 33.116091 4571.323946
4 37.5 0.75 0.25 0.25 0.25 94.6 57.2 79.0 68.2 33.0 55.9 24.0 0.39 0.360996 0.388860 29.558019 3968.547011

在建立随机森林模型时候，主要采用的是回归模型，与分类模型有一些差异，随机森林模型具有高准确性、抗过拟合、较好的鲁棒性，所以这里也选择用随机森林模型进行预测，并且通过随机搜索后，选择了参数范围内的最优参数，优化后的随机森林模型在均方误差和绝对系数上都优于多元线性回归模型，最后采用优化后的随机森林模型预测测试集。

9.总结

本项目主要采用了三个模型，一个k均值聚类模型，主要通过聚类得出不同蓝莓的聚合情况，主要将蓝莓分成4类，其中：
聚类0和聚类2具有相似的果实特性和产量，但聚类0的气温范围较低，降雨天数少，表明在较干燥和凉爽的环境中，通过优化克隆大小和蜜蜂密度可能能获得高产量。
聚类1在所有聚类中气温范围最高且降雨最多，但产量最低，这可能说明过高的温度和过多的降雨对蓝莓的生产不利。
聚类3有着适中的气温和降雨，较大的克隆大小和高蜜蜂密度，其产量和果实质量表明这可能是一个相对平衡的种植环境。
第二个模型是多元线性回归模型，通过相关性分析可以得知，产量与其他变量存在线性关系，所以可以采取多元线性回归模型，但是其他变量存在多重共线性，可以采用主成分分析降维，虽然这样做会降低模型的可解释性，但是能够保障建立出来的预测模型更加准确，最后将得到的模型进行残差项检验，发现随机误差项（残差）是符合基本假设的，模型比较好。
第三个模型就是随机森林模型，这是老朋友了，在回归和分类问题上，我比较喜欢使用它，虽然它的解释性比较差，有时候运行比较久，但是得到的模型是比较不错的，预测也比较准确，优化前的随机森林模型与多元线性回归模型差异不大，但是优化后的随机森林模型就比多元线性回归模型好一些，最后还输出了模型中重要特征：发现果实集（fruitset）是最重要的特征，可能是因为：果实集是产量形成的直接因素。如果果实集比率高，意味着更多的花朵被成功受粉并开始形成果实，最终导致更高的产量，与相关性分析得到的结果一致，果实集对产量具有高强度正相关。

字段名	描述
id	蓝莓唯一标识
Clonesize	蓝莓克隆平均大小,单位：m²
Honeybee	蜜蜂密度（单位：蜜蜂/m²/分钟）
Bumbles	大型蜜蜂密度（单位：大型蜜蜂/m²/分钟）
Andrena	安德烈纳蜂密度（单位：安德烈纳蜂/m²/分钟）
Osmia	钥匙蜂密度（单位：钥匙蜂/m²/分钟）
MaxOfUpperTRange	花期内最高温带日平均气温的最高记录, 单位: °C
MinOfUpperTRange	花期内最高温带日平均气温的最低记录, 单位: °C
AverageOfUpperTRange	花期内最高温带日平均气温, 单位: °C
MaxOfLowerTRange	花期内最低温带日平均气温的最高记录, 单位: °C
MinOfLowerTRange	花期内最低温带日平均气温的最低记录, 单位: °C
AverageOfLowerTRange	花期内最低温带日平均气温, 单位: °C
RainingDays	花期内降雨量大于 0 的日数总和, 单位: 天
AverageRainingDays	花期内降雨日数的平均值, 单位: 天
fruitset	果实集
fruitmass	果实质量
seeds	种子数
yield	产量

	clonesize	honeybee	bumbles	andrena	osmia	MaxOfUpperTRange	MinOfUpperTRange	AverageOfUpperTRange	MaxOfLowerTRange	MinOfLowerTRange	AverageOfLowerTRange	RainingDays	AverageRainingDays	fruitset	fruitmass	seeds	yield
Cluster
0	17.328080	0.346109	0.296424	0.498191	0.612562	74.767492	45.189268	62.475436	53.893253	26.075593	44.225397	11.916085	0.228554	0.555068	0.469118	38.871086	6961.936706
1	22.270245	0.437966	0.275738	0.484365	0.560719	90.494726	54.712893	75.610762	65.240384	31.567928	53.465477	26.255506	0.432669	0.436473	0.418616	32.748693	4863.325850
2	17.650084	0.351646	0.295031	0.498170	0.612210	89.643195	54.205118	74.910835	64.626788	31.272423	52.962712	12.246797	0.232205	0.559047	0.477165	39.367173	6975.032587
3	22.004267	0.429347	0.278173	0.488991	0.580457	72.474424	43.797838	60.539772	52.218265	25.271977	42.855192	25.545684	0.422139	0.449389	0.415332	33.048957	5113.166299

	feature	VIF
0	const	452.702544
1	clonesize	1.376754
2	honeybee	1.120483
3	bumbles	1.151751
4	andrena	1.186737
5	osmia	1.200157
6	MaxOfUpperTRange	3414.130675
7	MinOfUpperTRange	559.970976
8	AverageOfUpperTRange	30443.904985
9	MaxOfLowerTRange	1105.083282
10	MinOfLowerTRange	24965.147210
11	AverageOfLowerTRange	23062.819227
12	RainingDays	55.514009
13	AverageRainingDays	56.054862
14	fruitset	11.967936
15	fruitmass	12.726585
16	seeds	9.797380

	clonesize	honeybee	bumbles	andrena	osmia	MaxOfUpperTRange	MinOfUpperTRange	AverageOfUpperTRange	MaxOfLowerTRange	MinOfLowerTRange	AverageOfLowerTRange	RainingDays	AverageRainingDays	fruitset	fruitmass	seeds	predicted_yield
0	25.0	0.25	0.25	0.25	0.25	86.0	52.0	71.9	62.0	30.0	50.8	24.0	0.39	0.399367	0.408088	31.394569	4403.677514
1	12.5	0.25	0.25	0.75	0.63	94.6	57.2	79.0	68.2	33.0	55.9	1.0	0.10	0.488048	0.442866	36.846956	6141.142129
2	12.5	0.25	0.25	0.63	0.63	86.0	52.0	71.9	62.0	30.0	50.8	16.0	0.26	0.583379	0.487057	40.037644	7279.014537
3	25.0	0.50	0.38	0.38	0.63	86.0	52.0	71.9	62.0	30.0	50.8	16.0	0.26	0.433014	0.422847	33.116091	5133.609279
4	37.5	0.75	0.25	0.25	0.25	94.6	57.2	79.0	68.2	33.0	55.9	24.0	0.39	0.360996	0.388860	29.558019	3598.218257

	Feature	Importance
13	fruitset	0.804402
15	seeds	0.149529
14	fruitmass	0.019532
4	osmia	0.003723
3	andrena	0.003491
12	AverageRainingDays	0.003144
11	RainingDays	0.002704
2	bumbles	0.002239
1	honeybee	0.001963
10	AverageOfLowerTRange	0.001536
0	clonesize	0.001513
6	MinOfUpperTRange	0.001435
8	MaxOfLowerTRange	0.001392
9	MinOfLowerTRange	0.001185
5	MaxOfUpperTRange	0.001141
7	AverageOfUpperTRange	0.001072