改前版本 直接从 word 转的 凑合着看吧

# 代码解析

# Before Problem -- 数据传入与预处理，全局函数设置

# 导入必要的库

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.linear_model import LinearRegression
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.metrics import accuracy_score
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LogisticRegression
from sklearn import preprocessing
import joblib
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA
plt.rcParams['font.sans-serif'] = [u'simHei']   # 显示中文
plt.rcParams['axes.unicode_minus'] = False      # 解决负号问题

# 数据读入

# 从 train.csv 读入 train 数据

1
2

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)

# 数据预处理

id 列对于数据分析没有用处，直接清洗掉这一列

1

train_data = train_data.drop(columns='id')

经检验，原始数据中不存在缺失值与明显的离群点，因此在这里全部保留

将 ' 洪水概率 ' 设置为目标变量，将原始数据划分为 X，y

1
2
3

X=train_data.drop(columns='洪水概率')
y=train_data['洪水概率']
y=np.array(y)

# 数据相关性简要分析

利用不同指标与洪水概率的相关系数来衡量这些指标与洪水概率的关系

1
2
3
4

train_data['洪水概率']=y
heat_corr = train_data.corr()['洪水概率']
heat_corr_sorted=heat_corr.sort_values(ascending=False)
heat_corr=heat_corr[:-1]

可视化相关系数图

1
2
3
4
5
6
7

plt.figure(figsize=(10, 6))
heat_corr.plot(kind='bar', color='skyblue')
plt.title('不同因素对洪水概率的影响')
plt.xlabel('影响因素')
plt.ylabel('corr')
plt.grid(axis='y')
plt.show()

# 建立三个水平的指标集

20： 全指标

10：半指标

5：核心指标

划分依据: 相关系数大小

根据语义与实际情况考虑，先剔除如下的两个属性形成 18 个指标的属性集

湿地损失、无效防灾

进一步建立精简指标集，剔除如下的指标集：

湿地损失、无效防灾、农业实践、排水系统、规划不足、季风强度、政策因素、流域、基础设施恶化

1

X.columns

Index(['季风强度', '地形排水', '河流管理', '森林砍伐', '城市化', '气候变化', '大坝质量', '淤积', '农业实践',
       '侵蚀', '无效防灾', '排水系统', '海岸脆弱性', '滑坡', '流域', '基础设施恶化', '人口得分', '湿地损失',
       '规划不足', '政策因素'],
      dtype='object')

1
2
3

X_COPY=X
eighteen_columns=X.drop(columns=['湿地损失','无效防灾'])
eleven_columns=X.drop(columns=['湿地损失','无效防灾','农业实践','排水系统','规划不足','季风强度','政策因素','流域','基础设施恶化'])

# 利用过滤法获取指标集

1

from sklearn.feature_selection import SelectKBest,f_classif

1
2
3
4
5
6

ten_columns=['淤积', '滑坡', '城市化', '海岸脆弱性', '气候变化', '侵蚀', '人口得分', '大坝质量', '河流管理', '地形排水']
X_selected = X[ten_columns]
selector = SelectKBest(score_func=f_classif, k=4)
X_new = selector.fit_transform(X_selected, y)
selected_features_indices = selector.get_support(indices=True)
selected_features = X_selected.columns[selected_features_indices]

1
2
3

selected_features=list(selected_features)
selected_features.append('森林砍伐')
five_columns=selected_features

# 全局函数设置

建立余弦相似度的计算函数

cosine(a,b)

input : 向量A，B

output : 两个向量的余弦相似度 0-1

1
2
3
4
5
6

def cosine(a,b):
    dot_product = np.dot(a, b)
    norm_a = np.linalg.norm(a)
    norm_b = np.linalg.norm(b)
    cosine_similarity = dot_product / (norm_a * norm_b) 
    return cosine_similarity

# 划分训练集与测试集

1

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Problem ONE -- 全指标集数学模型建立与求解

利用多种数学模型拟合本问题，并得出不同指标的权重
多元线性回归
梯度提升树
决策树
随机森林模型

# 模型一 线性回归

原理：利用标准方程法求解线性方程组

1
2
3
4
5
6
7
8
9
10
11
12
13

from sklearn.metrics import accuracy_score
model_line=LinearRegression()
model_line.fit(X_train,y_train)
y_pred = model_line.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = model_line.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
print(f"模型的系数为: {model_line.coef_}")
print(f"模型的截距为: {model_line.intercept_}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')
joblib.dump(model_line,'../Model/linear.pkl')

Model score: 0.8444528508563269
均方误差为: 0.0004043169737229802
模型的系数为: [0.00560656 0.00564613 0.00565985 0.00567684 0.00565284 0.0056578
 0.00565561 0.00564293 0.0056437  0.0056492  0.00564238 0.0056413
 0.005674   0.00564575 0.0056467  0.00562068 0.00568107 0.00563092
 0.00561738 0.00564302]
模型的截距为: -0.053357024379501294
余弦相似度为0.9992129608069907





['../Model/linear.pkl']

# 模型二 梯度提升树

原理：递归的决策树，能发现更深层的变量关系

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

from sklearn.ensemble import GradientBoostingRegressor
import numpy as np
import matplotlib.pyplot as plt
from sklearn.tree import plot_tree
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_val_score
n_estimators = 2000
gb_model = GradientBoostingRegressor(n_estimators=n_estimators, random_state=42)
cv_scores = cross_val_score(gb_model, X_train, y_train, cv=5)
print("每折交叉验证的评分:", cv_scores)
gb_model.fit(X_train, y_train)
gb_score  = gb_model.score(X_test,y_test)
print(f'Model Score : {gb_score}')
feature_importances = gb_model.feature_importances_
print("特征重要性:", feature_importances)
train_loss = np.zeros((gb_model.n_estimators,), dtype=np.float64)
for i, y_pred in enumerate(gb_model.staged_predict(X_train)):
    train_loss[i] = mean_squared_error(y_train, y_pred)

test_loss = np.zeros((gb_model.n_estimators,), dtype=np.float64)
for i, y_pred in enumerate(gb_model.staged_predict(X_test)):
    test_loss[i] = mean_squared_error(y_test, y_pred)
y_pred = gb_model.predict(X_test)
theta = cosine(y_pred,y_test)
print(f'余弦相似度为{theta}')
mse = mean_squared_error(y_test, y_pred)
print(f'均方误差为{mse}')
plt.figure(figsize=(10, 6))
plt.plot(np.arange(1, gb_model.n_estimators + 1), train_loss, label='Training Loss')
plt.plot(np.arange(1, gb_model.n_estimators + 1), test_loss, label='Test Loss')
plt.xlabel('Number of Iterations')
plt.ylabel('Mean Squared Error')
plt.title('Training and Test Loss')
plt.legend()
plt.show()
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances)), feature_importances, tick_label=[f'{X.columns[i]}' for i in range(0, len(feature_importances))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()
joblib.dump(gb_model,'../Model/gb_tree.pkl')

---------------------------------------------------------------------------

KeyboardInterrupt                         Traceback (most recent call last)

Cell In[35], line 9
      7 n_estimators = 2000
      8 gb_model = GradientBoostingRegressor(n_estimators=n_estimators, random_state=42)
----> 9 cv_scores = cross_val_score(gb_model, X_train, y_train, cv=5)
     10 print("每折交叉验证的评分:", cv_scores)
     11 gb_model.fit(X_train, y_train)


File f:\a\envs\env\lib\site-packages\sklearn\model_selection\_validation.py:515, in cross_val_score(estimator, X, y, groups, scoring, cv, n_jobs, verbose, fit_params, pre_dispatch, error_score)
    512 # To ensure multimetric format is not supported
    513 scorer = check_scoring(estimator, scoring=scoring)
--> 515 cv_results = cross_validate(
    516     estimator=estimator,
    517     X=X,
    518     y=y,
    519     groups=groups,
    520     scoring={"score": scorer},
    521     cv=cv,
    522     n_jobs=n_jobs,
    523     verbose=verbose,
    524     fit_params=fit_params,
    525     pre_dispatch=pre_dispatch,
    526     error_score=error_score,
    527 )
    528 return cv_results["test_score"]


File f:\a\envs\env\lib\site-packages\sklearn\model_selection\_validation.py:266, in cross_validate(estimator, X, y, groups, scoring, cv, n_jobs, verbose, fit_params, pre_dispatch, return_train_score, return_estimator, error_score)
    263 # We clone the estimator to make sure that all the folds are
    264 # independent, and that it is pickle-able.
    265 parallel = Parallel(n_jobs=n_jobs, verbose=verbose, pre_dispatch=pre_dispatch)
--> 266 results = parallel(
    267     delayed(_fit_and_score)(
    268         clone(estimator),
    269         X,
    270         y,
    271         scorers,
    272         train,
    273         test,
    274         verbose,
    275         None,
    276         fit_params,
    277         return_train_score=return_train_score,
    278         return_times=True,
    279         return_estimator=return_estimator,
    280         error_score=error_score,
    281     )
    282     for train, test in cv.split(X, y, groups)
    283 )
    285 _warn_or_raise_about_fit_failures(results, error_score)
    287 # For callabe scoring, the return type is only know after calling. If the
    288 # return type is a dictionary, the error scores can now be inserted with
    289 # the correct key.


File f:\a\envs\env\lib\site-packages\sklearn\utils\parallel.py:63, in Parallel.__call__(self, iterable)
     58 config = get_config()
     59 iterable_with_config = (
     60     (_with_config(delayed_func, config), args, kwargs)
     61     for delayed_func, args, kwargs in iterable
     62 )
---> 63 return super().__call__(iterable_with_config)


File f:\a\envs\env\lib\site-packages\joblib\parallel.py:1051, in Parallel.__call__(self, iterable)
   1048 if self.dispatch_one_batch(iterator):
   1049     self._iterating = self._original_iterator is not None
-> 1051 while self.dispatch_one_batch(iterator):
   1052     pass
   1054 if pre_dispatch == "all" or n_jobs == 1:
   1055     # The iterable was consumed all at once by the above for loop.
   1056     # No need to wait for async callbacks to trigger to
   1057     # consumption.


File f:\a\envs\env\lib\site-packages\joblib\parallel.py:864, in Parallel.dispatch_one_batch(self, iterator)
    862     return False
    863 else:
--> 864     self._dispatch(tasks)
    865     return True


File f:\a\envs\env\lib\site-packages\joblib\parallel.py:782, in Parallel._dispatch(self, batch)
    780 with self._lock:
    781     job_idx = len(self._jobs)
--> 782     job = self._backend.apply_async(batch, callback=cb)
    783     # A job can complete so quickly than its callback is
    784     # called before we get here, causing self._jobs to
    785     # grow. To ensure correct results ordering, .insert is
    786     # used (rather than .append) in the following line
    787     self._jobs.insert(job_idx, job)


File f:\a\envs\env\lib\site-packages\joblib\_parallel_backends.py:208, in SequentialBackend.apply_async(self, func, callback)
    206 def apply_async(self, func, callback=None):
    207     """Schedule a func to be run"""
--> 208     result = ImmediateResult(func)
    209     if callback:
    210         callback(result)


File f:\a\envs\env\lib\site-packages\joblib\_parallel_backends.py:572, in ImmediateResult.__init__(self, batch)
    569 def __init__(self, batch):
    570     # Don't delay the application, to avoid keeping the input
    571     # arguments in memory
--> 572     self.results = batch()


File f:\a\envs\env\lib\site-packages\joblib\parallel.py:263, in BatchedCalls.__call__(self)
    259 def __call__(self):
    260     # Set the default nested backend to self._backend but do not set the
    261     # change the default number of processes to -1
    262     with parallel_backend(self._backend, n_jobs=self._n_jobs):
--> 263         return [func(*args, **kwargs)
    264                 for func, args, kwargs in self.items]


File f:\a\envs\env\lib\site-packages\joblib\parallel.py:263, in <listcomp>(.0)
    259 def __call__(self):
    260     # Set the default nested backend to self._backend but do not set the
    261     # change the default number of processes to -1
    262     with parallel_backend(self._backend, n_jobs=self._n_jobs):
--> 263         return [func(*args, **kwargs)
    264                 for func, args, kwargs in self.items]


File f:\a\envs\env\lib\site-packages\sklearn\utils\parallel.py:123, in _FuncWrapper.__call__(self, *args, **kwargs)
    121     config = {}
    122 with config_context(**config):
--> 123     return self.function(*args, **kwargs)


File f:\a\envs\env\lib\site-packages\sklearn\model_selection\_validation.py:686, in _fit_and_score(estimator, X, y, scorer, train, test, verbose, parameters, fit_params, return_train_score, return_parameters, return_n_test_samples, return_times, return_estimator, split_progress, candidate_progress, error_score)
    684         estimator.fit(X_train, **fit_params)
    685     else:
--> 686         estimator.fit(X_train, y_train, **fit_params)
    688 except Exception:
    689     # Note fit time as time until error
    690     fit_time = time.time() - start_time


File f:\a\envs\env\lib\site-packages\sklearn\ensemble\_gb.py:538, in BaseGradientBoosting.fit(self, X, y, sample_weight, monitor)
    535     self._resize_state()
    537 # fit the boosting stages
--> 538 n_stages = self._fit_stages(
    539     X,
    540     y,
    541     raw_predictions,
    542     sample_weight,
    543     self._rng,
    544     X_val,
    545     y_val,
    546     sample_weight_val,
    547     begin_at_stage,
    548     monitor,
    549 )
    551 # change shape of arrays after fit (early-stopping or additional ests)
    552 if n_stages != self.estimators_.shape[0]:


File f:\a\envs\env\lib\site-packages\sklearn\ensemble\_gb.py:615, in BaseGradientBoosting._fit_stages(self, X, y, raw_predictions, sample_weight, random_state, X_val, y_val, sample_weight_val, begin_at_stage, monitor)
    608     old_oob_score = loss_(
    609         y[~sample_mask],
    610         raw_predictions[~sample_mask],
    611         sample_weight[~sample_mask],
    612     )
    614 # fit next stage of trees
--> 615 raw_predictions = self._fit_stage(
    616     i,
    617     X,
    618     y,
    619     raw_predictions,
    620     sample_weight,
    621     sample_mask,
    622     random_state,
    623     X_csc,
    624     X_csr,
    625 )
    627 # track deviance (= loss)
    628 if do_oob:


File f:\a\envs\env\lib\site-packages\sklearn\ensemble\_gb.py:257, in BaseGradientBoosting._fit_stage(self, i, X, y, raw_predictions, sample_weight, sample_mask, random_state, X_csc, X_csr)
    254     sample_weight = sample_weight * sample_mask.astype(np.float64)
    256 X = X_csr if X_csr is not None else X
--> 257 tree.fit(X, residual, sample_weight=sample_weight, check_input=False)
    259 # update tree leaves
    260 loss.update_terminal_regions(
    261     tree.tree_,
    262     X,
   (...)
    269     k=k,
    270 )


File f:\a\envs\env\lib\site-packages\sklearn\tree\_classes.py:1247, in DecisionTreeRegressor.fit(self, X, y, sample_weight, check_input)
   1218 def fit(self, X, y, sample_weight=None, check_input=True):
   1219     """Build a decision tree regressor from the training set (X, y).
   1220 
   1221     Parameters
   (...)
   1244         Fitted estimator.
   1245     """
-> 1247     super().fit(
   1248         X,
   1249         y,
   1250         sample_weight=sample_weight,
   1251         check_input=check_input,
   1252     )
   1253     return self


File f:\a\envs\env\lib\site-packages\sklearn\tree\_classes.py:379, in BaseDecisionTree.fit(self, X, y, sample_weight, check_input)
    368 else:
    369     builder = BestFirstTreeBuilder(
    370         splitter,
    371         min_samples_split,
   (...)
    376         self.min_impurity_decrease,
    377     )
--> 379 builder.build(self.tree_, X, y, sample_weight)
    381 if self.n_outputs_ == 1 and is_classifier(self):
    382     self.n_classes_ = self.n_classes_[0]


KeyboardInterrupt: 

# 模型三 决策树

原理：利用信息增益实现节点的划分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score
from sklearn.tree import plot_tree
from sklearn.metrics import mean_squared_error
tree_model = DecisionTreeRegressor(random_state=42)
cv_scores = cross_val_score(tree_model, X_train, y_train, cv=5)
print("每折交叉验证的评分:", cv_scores)
print("交叉验证评分的平均值:", np.mean(cv_scores))
tree_model.fit(X_train, y_train)
tree_score = tree_model.score(X_test, y_test)
print("Model Score :", tree_score)
feature_importances = tree_model.feature_importances_
print("特征重要性:", feature_importances)
feature_importances_flat = np.ravel(feature_importances)
# 可视化特征重要性
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances_flat)), feature_importances_flat, tick_label=[f'{X.columns[i]}' for i in range(0, len(feature_importances_flat))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()
# 在测试集上进行预测
y_pred = tree_model.predict(X_test)
theta = cosine(y_pred,y_test)
print(f'余弦相似度为{theta}')
mse = mean_squared_error(y_test, y_pred)
print(f'均方误差为{mse}')
joblib.dump(tree_model,'../Model/ju_tree.pkl')

每折交叉验证的评分: [0.04694866 0.04297524 0.0450139  0.04790429 0.04502323]
交叉验证评分的平均值: 0.045573064737767635
Model Score : 0.05301513744804842
特征重要性: [0.05072873 0.04679774 0.04822921 0.0498514  0.05054892 0.05003494
 0.04871661 0.0509244  0.05185402 0.04976417 0.04994974 0.0515037
 0.05054295 0.05088047 0.04947362 0.04914221 0.04920693 0.0502652
 0.05116479 0.05042025]

余弦相似度为0.9952026345820846
均方误差为0.0024615176549126197





['../Model/ju_tree.pkl']

# 模型四 随机森林

原理： 抽样决策树并使用投票策略

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

rf = RandomForestRegressor(n_estimators=4, random_state=42) 
rf.fit(X_train, y_train)
y_pred = rf.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = rf.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')
feature_importances = rf.feature_importances_
print("特征重要性:", feature_importances)
feature_importances_flat = np.ravel(feature_importances)
# 可视化特征重要性
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances_flat)), feature_importances_flat, tick_label=[f'{X_train.columns[i]}' for i in range(0, len(feature_importances_flat))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()
import joblib
joblib.dump(rf,'../Model/rf_forest.pkl')

Model score: 0.499235508143682
均方误差为: 0.0013016476676680254
余弦相似度为0.9974668593808351
特征重要性: [0.04968638 0.04813227 0.04933601 0.05021542 0.05116166 0.0494662
 0.05073097 0.05094782 0.05060163 0.05058865 0.04939483 0.05095208
 0.05029554 0.05063372 0.04975227 0.04803734 0.04872015 0.05049445
 0.05024486 0.05060773]

['../Model/rf_forest.pkl']

# Problem TWO 风险聚类

# 读入数据

1
2

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)

# 聚类算法：Kmeans

原理: 利用n维空间距离对初始聚类中心进行学习，得到最终聚类中心

# 聚类预测函数

predict_cluster(kmeans,points)
input: kmeans模型 points需要预测的点集
output: 点集的预测结果

1
2
3
4
5

def predict_cluster(kmeans, points):
    points = np.array(points)
    points = points.reshape(-1,1)
    cluster_label = kmeans.predict(points)
    return cluster_label

# 建立 Kmeans 模型

1
2
3
4
5
6
7
8
9
10
11
12

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
from sklearn.cluster import KMeans
y=train_data['洪水概率']
y=np.array(y)
#y=-np.log(1/y-1)
data=np.array(y)
data_2d = data.reshape(-1, 1)
kmeans = KMeans(n_clusters=3)
kmeans.fit(data_2d)
cluster_centers = kmeans.cluster_centers_
print(cluster_centers)

f:\a\envs\env\lib\site-packages\sklearn\cluster\_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(


[[0.56441718]
 [0.44231226]
 [0.50212141]]

# 建立聚类中心到风险的映射字典

1
2
3
4
5
6

temp = np.array(cluster_centers)
indices= np.argsort(temp,axis=0)
dic={}
dic[f'{indices[0]}']='低风险'
dic[f'{indices[1]}']='中风险'
dic[f'{indices[2]}']='高风险'

# 示例

data= [0.55,0.99,0.51,0.42]

1
2

data= [0.55,0.99,0.51,0.42]
print(predict_cluster(kmeans,data))

[1 1 2 0]

# 基于风险预警的预测模型

# 模型建立 MyModel

内核    ：某一个概率预测模型与kmeans聚类模型
predict : 预测输入的X对应的y的风险类别
score   : 利用预测成功的比例衡量得分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

class MyModel:
    def __init__(self,model,kmeans) -> None:
        self.model = model
        self.kmeans = kmeans
    def predict(self,X_test) -> np.ndarray:
        y_pred=self.model.predict(X_test)
        y_pred=np.array(y_pred)
        y_pred=predict_cluster(kmeans,y_pred)
        y_pred=np.array(y_pred)
        return y_pred
    def score(self,X_test,y_test):
        y_pred=self.predict(X_test)
        y_test = np.array(y_test)
        y_true=predict_cluster(kmeans,y_test)
        correct_predictions = np.sum(y_pred == y_true)
        overlap = correct_predictions / y_pred.size
        return overlap

# 数据读入

1
2
3
4
5
6

def map_dic(dic,a):
    a=np.array(a)
    lis=[]
    for i in a:
        lis.append(dic[f'[{i}]'])
    return np.array(lis)

1
2
3
4
5
6
7
8
9
10
11
12
13
14

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data_copy=train_data
train_data_copy['洪水概率'] = predict_cluster(kmeans,np.array(train_data_copy['洪水概率']))
train_data_copy['洪水概率'] = map_dic(dic,train_data_copy['洪水概率'])
train_data_copy.to_csv('../Data/train_dealed.csv',index=False,encoding='utf-8')
train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data=train_data.drop(columns='id')

X=train_data.drop(columns='洪水概率')
y=train_data['洪水概率']
y=np.array(y)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 利用先前模型构建自定义的分类模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import os
models = ['linear.pkl','gb_tree.pkl','ju_tree.pkl','rf_forest.pkl']
path = '../Model'
m_path='../Model/3_means'
for model in models:
    model_name=model.replace('.pkl','')
    print(f'当前模型为{model}')
    model_path=os.path.join(path,model)
    model = joblib.load(model_path)
    model = MyModel(model,kmeans)
    # 这里可以加入模型的预测过程
    y_pred = model.predict(X_test)
    y_true=predict_cluster(kmeans,y_test)
    mse = mean_squared_error(y_true, y_pred)
    score = model.score(X_test,y_test)
    print(f'Model score: {score}')
    print(f"均方误差为: {mse}")
    accuracy = cosine(y_true,y_pred)
    print(f'余弦相似度为{accuracy}')
    if model_name == 'linear':
        print(f'不同指标的权重为')
        print(model.model.coef_)
    else:
        print(f'不同指标的权重为')
        print(model.model.feature_importances_)
    joblib.dump(model,m_path+'/'+model_name+'.pkl')

当前模型为linear.pkl
Model score: 0.7765729680757218
均方误差为: 0.6680161171113177
余弦相似度为0.8397624947557178
不同指标的权重为
[0.00560656 0.00564613 0.00565985 0.00567684 0.00565284 0.0056578
 0.00565561 0.00564293 0.0056437  0.0056492  0.00564238 0.0056413
 0.005674   0.00564575 0.0056467  0.00562068 0.00568107 0.00563092
 0.00561738 0.00564302]
当前模型为gb_tree.pkl
Model score: 0.7688148201130105
均方误差为: 0.6767756240612259
余弦相似度为0.8430676769184631
不同指标的权重为
[0.05053843 0.0513698  0.05097438 0.04985333 0.04968169 0.0498542
 0.05148326 0.05025949 0.04944183 0.04878379 0.04970291 0.04908077
 0.0490136  0.05042491 0.04939282 0.05113171 0.05085744 0.04953673
 0.04897519 0.04964374]
当前模型为ju_tree.pkl
Model score: 0.5086617552392533
均方误差为: 1.1945068306988056
余弦相似度为0.7075274653380997
不同指标的权重为
[0.05072873 0.04679774 0.04822921 0.0498514  0.05054892 0.05003494
 0.04871661 0.0509244  0.05185402 0.04976417 0.04994974 0.0515037
 0.05054295 0.05088047 0.04947362 0.04914221 0.04920693 0.0502652
 0.05116479 0.05042025]
当前模型为rf_forest.pkl
Model score: 0.5983406051069309
均方误差为: 1.0319147414348044
余弦相似度为0.7889832533919223
不同指标的权重为
[0.04968638 0.04813227 0.04933601 0.05021542 0.05116166 0.0494662
 0.05073097 0.05094782 0.05060163 0.05058865 0.04939483 0.05095208
 0.05029554 0.05063372 0.04975227 0.04803734 0.04872015 0.05049445
 0.05024486 0.05060773]

# 逻辑回归

1
2
3
4
5
6
7
8

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data=train_data.drop(columns='id')
X=train_data.drop(columns='洪水概率')
y=train_data['洪水概率']
y=np.array(y)
y=predict_cluster(kmeans,y)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

1
2
3
4

model = LogisticRegression(solver='liblinear')
model.fit(X_train,y_train)
y_pred = model.predict(X_test)

1
2
3
4
5
6
7
8

y_true=y_test
mse = mean_squared_error(y_true, y_pred)
score =np.sum(y_true==y_pred)/len(y_pred)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
accuracy = cosine(y_true,y_pred)
print(f'余弦相似度为{accuracy}')
joblib.dump(model,'../Model/3_means/logistic.pkl')

Model score: 0.7634599337195718
均方误差为: 0.5021147748134372
余弦相似度为0.8741869419358508





['../Model/3_means/logistic.pkl']

# 模型灵敏度分析

利用结果的熵值来衡量模型的灵敏度

1
2
3
4

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data = train_data.drop(columns='id')
X=train_data.drop(columns='洪水概率')

# 设置熵的计算函数 calculate_entropy

input: list
ouput :总熵

1
2
3
4
5
6
7

def calculate_entropy(lst):
    value_counts = np.unique(lst, return_counts=True)
    values, counts = value_counts
    total_count = len(lst)
    probabilities = counts / total_count
    entropy = -np.sum(probabilities * np.log2(probabilities))
    return entropy

# 设置基于属性变化率的计算函数 average_change_rate

input:list
output: 平均相邻变化率

1
2
3
4
5
6
7

def average_change_rate(lst):
    # 检查列表是否至少有两项
    if len(lst) < 2:
        return None
    change_rates = [(lst[i+1] - lst[i]) / lst[i] for i in range(len(lst) - 1)]
    avg_change_rate = sum(change_rates) / len(change_rates)
    return avg_change_rate

# 启动基于属性对概率变化率的灵敏度计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import os
models = ['linear.pkl','gb_tree.pkl','ju_tree.pkl','rf_forest.pkl']
path='../Model'
m_path='../Model/3_means'
entropy_list=[]
for model in models:
    model_copy=model
    train_data_path='../Data/train.csv'
    train_data = pd.read_csv(train_data_path)
    train_data = train_data.drop(columns='id')
    X=train_data.drop(columns='洪水概率')
    model_name=model.replace('.pkl','')
    print(f'当前模型为{model}')
    # 读入模型
    model_path=os.path.join(path,model)
    model = joblib.load(model_path)
    av_list=[]
    # 计算不同模型的变化率
    for each_column in X.columns: 
        temp_X=pd.DataFrame(columns=X.columns)
        target=np.array(X[each_column])
        min=np.min(target)
        max=np.max(target)
        column_list=[]
        for i in range(min,max+1):
            temp_X_copy=X
            temp_X_copy[each_column]=i
            y_pred=model.predict(temp_X_copy)
            column_list.append(y_pred.mean())
        av_list.append(average_change_rate(column_list))
    df = pd.DataFrame({
    '属性': X.columns,
    '变化率': av_list
    })
    df.plot(kind='bar', x='属性', y='变化率')
    plt.title(f'{model_name}属性平均变化率')
    plt.xlabel('属性')
    plt.ylabel('变化率')
    plt.show()

当前模型为linear.pkl

当前模型为gb_tree.pkl

当前模型为ju_tree.pkl

当前模型为rf_forest.pkl

# 启动基于预测结果的熵值的灵敏度计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

import os
models = ['linear.pkl','gb_tree.pklfvv','ju_tree.pkl','rf_forest.pkl','logistic.pkl']
path='../Model'
m_path='../Model/3_means'
entropy_list=[]
for model in models:
    model_copy=model
    train_data_path='../Data/train.csv'
    train_data = pd.read_csv(train_data_path)
    train_data = train_data.drop(columns='id')
    X=train_data.drop(columns='洪水概率')
    model_name=model.replace('.pkl','')
    print(f'当前模型为{model}')
    model_path=os.path.join(m_path,model_copy)
    model = joblib.load(model_path)
    train_data_path='../Data/train.csv'
    train_data = pd.read_csv(train_data_path)
    train_data = train_data.drop(columns='id')
    X=train_data.drop(columns='洪水概率')
    train_data = pd.read_csv(train_data_path)
    temp_X=X
    y_pred=model.predict(X)
    print(y_pred)
    temp_X['洪水概率']=y_pred
    temp_X['id']=train_data['id']
    for i in range(3):    
        temp=temp_X[temp_X['洪水概率']==i]
        temp['洪水概率']=dic[f'[{i}]']
        name=dic[f'[{i}]']
        csv_path='../Data/3_means/'
        temp.to_csv(csv_path+f'{model_name}'+'/'+f'{name}.csv',index=False,encoding='utf-8')
    entropy_list.append(calculate_entropy(y_pred))
df = pd.DataFrame({
    '模型': ['线性回归','梯度提升树','决策树','随机森林','逻辑回归'],
    '熵值': entropy_list
})
df.plot(kind='bar', x='模型', y='熵值')
plt.title(f'不同模型的灵敏度')
plt.xlabel('模型')
plt.ylabel('熵值')
plt.show()
    

当前模型为linear.pkl
[2 2 2 ... 2 2 1]


C:\Users\12780\AppData\Local\Temp\ipykernel_25292\1261529529.py:28: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  temp['洪水概率']=dic[f'[{i}]']


当前模型为gb_tree.pkl
[2 2 2 ... 2 2 1]


C:\Users\12780\AppData\Local\Temp\ipykernel_25292\1261529529.py:28: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  temp['洪水概率']=dic[f'[{i}]']


当前模型为ju_tree.pkl
[2 0 2 ... 2 0 1]


C:\Users\12780\AppData\Local\Temp\ipykernel_25292\1261529529.py:28: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  temp['洪水概率']=dic[f'[{i}]']


当前模型为rf_forest.pkl
[2 2 2 ... 2 2 2]


C:\Users\12780\AppData\Local\Temp\ipykernel_25292\1261529529.py:28: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  temp['洪水概率']=dic[f'[{i}]']


当前模型为logistic.pkl
[2 2 2 ... 2 2 1]


C:\Users\12780\AppData\Local\Temp\ipykernel_25292\1261529529.py:28: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  temp['洪水概率']=dic[f'[{i}]']

# Problem THREE 从不同指标集建立洪水的预测模型

# 利用半指标集重新训练模型

# 数据读入与划分

1
2
3
4
5
6
7
8

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data = train_data.drop(columns='id')
X=train_data.drop(columns='洪水概率')
X=X[ten_columns]
y=train_data['洪水概率']
y=np.array(y)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 多元线性回归

1
2
3
4
5
6
7
8
9
10
11
12

from sklearn.metrics import accuracy_score
model_line=LinearRegression()
model_line.fit(X_train,y_train)
y_pred = model_line.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = model_line.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
print(f"模型的系数为: {model_line.coef_}")
print(f"模型的截距为: {model_line.intercept_}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')

Model score: 0.37461847859805275
均方误差为: 0.0016255673314973993
模型的系数为: [0.00507833 0.00501385 0.00495665 0.00492889 0.00504095 0.00491355
 0.00504961 0.00504281 0.0050641  0.00502683]
模型的截距为: 0.2568960934135129
余弦相似度为0.9968319141241095

# 梯度提升树

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

from sklearn.ensemble import GradientBoostingRegressor
import numpy as np
import matplotlib.pyplot as plt
from sklearn.tree import plot_tree
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_val_score
n_estimators = 1200
gb_model = GradientBoostingRegressor(n_estimators=n_estimators, random_state=42)
gb_model.fit(X_train, y_train)
# 在测试集上评估模型性能
gb_score  = gb_model.score(X_test,y_test)
print(f'Model Score : {gb_score}')
feature_importances = gb_model.feature_importances_
print("特征重要性:", feature_importances)
# 计算每轮迭代的损失值
train_loss = np.zeros((gb_model.n_estimators,), dtype=np.float64)
for i, y_pred in enumerate(gb_model.staged_predict(X_train)):
    train_loss[i] = mean_squared_error(y_train, y_pred)

test_loss = np.zeros((gb_model.n_estimators,), dtype=np.float64)
for i, y_pred in enumerate(gb_model.staged_predict(X_test)):
    test_loss[i] = mean_squared_error(y_test, y_pred)
# 在测试集上进行预测
y_pred = gb_model.predict(X_test)
theta = cosine(y_pred,y_test)
print(f'余弦相似度为{theta}')
# 计算均方误差
mse = mean_squared_error(y_test, y_pred)
print(f'均方误差为{mse}')
# 绘制损失曲线
plt.figure(figsize=(10, 6))
plt.plot(np.arange(1, gb_model.n_estimators + 1), train_loss, label='Training Loss')
plt.plot(np.arange(1, gb_model.n_estimators + 1), test_loss, label='Test Loss')
plt.xlabel('Number of Iterations')
plt.ylabel('Mean Squared Error')
plt.title('Training and Test Loss')
plt.legend()
plt.show()
# 可视化特征重要性
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances)), feature_importances, tick_label=[f'{X_train.columns[i]}' for i in range(0, len(feature_importances))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()

Model Score : 0.378759314749781
特征重要性: [0.10139095 0.10060794 0.09787091 0.09610327 0.09900577 0.09537671
 0.10140741 0.10290733 0.10224735 0.10308238]
余弦相似度为0.996852927605829
均方误差为0.0016148039690650664

# 决策树

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score
from sklearn.tree import plot_tree
from sklearn.metrics import mean_squared_error
tree_model = DecisionTreeRegressor(random_state=42)
tree_model.fit(X_train, y_train)
tree_score = tree_model.score(X_test, y_test)
print("Model Score :", tree_score)
feature_importances = tree_model.feature_importances_
print("特征重要性:", feature_importances)
feature_importances_flat = np.ravel(feature_importances)
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances_flat)), feature_importances_flat, tick_label=[f'{X_train.columns[i]}' for i in range(0, len(feature_importances_flat))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()
y_pred = tree_model.predict(X_test)
theta = cosine(y_pred,y_test)
print(f'余弦相似度为{theta}')
mse = mean_squared_error(y_test, y_pred)
print(f'均方误差为{mse}')

Model Score : -0.37999290530770646
特征重要性: [0.09880106 0.09983387 0.10125522 0.10165829 0.10035667 0.10129578
 0.09900273 0.09949576 0.0994488  0.09885181]

余弦相似度为0.9930189485511745
均方误差为0.003587044560474931

# 随机森林

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

rf = RandomForestRegressor(n_estimators=4, random_state=42)  # 使用100棵树
rf.fit(X_train, y_train)
y_pred = rf.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = rf.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')
feature_importances = rf.feature_importances_
print("特征重要性:", feature_importances)
feature_importances_flat = np.ravel(feature_importances)
# 可视化特征重要性
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances_flat)), feature_importances_flat, tick_label=[f'{X_train.columns[i]}' for i in range(0, len(feature_importances_flat))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()

Model score: 0.15691875212778306
均方误差为: 0.002191438805653956
余弦相似度为0.9957272905379122
特征重要性: [0.09854545 0.10019614 0.10220218 0.1030608  0.1000937  0.10149794
 0.09828828 0.0987665  0.09905775 0.09829126]

# 利用核心指标集重新训练模型

# 数据读入

1
2
3
4
5
6
7
8
9
10

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data = train_data.drop(columns='id')
X=train_data.drop(columns='洪水概率')
# 在这里进行数据指标的选择
X=X[five_columns]
y=train_data['洪水概率']
y=np.array(y)
#y=-np.log(1/y-1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 多元线性回归

1
2
3
4
5
6
7
8
9
10
11
12

from sklearn.metrics import accuracy_score
model_line=LinearRegression()
model_line.fit(X_train,y_train)
y_pred = model_line.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = model_line.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
print(f"模型的系数为: {model_line.coef_}")
print(f"模型的截距为: {model_line.intercept_}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')

Model score: 0.181427445558773
均方误差为: 0.0021277328437480304
模型的系数为: [0.00482241 0.00476873 0.00479914 0.00474271 0.00475476]
模型的截距为: 0.38645331900338376
余弦相似度为0.9958512145383939

# 梯度提升树

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

from sklearn.ensemble import GradientBoostingRegressor
import numpy as np
import matplotlib.pyplot as plt
from sklearn.tree import plot_tree
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import cross_val_score
# 设置迭代次数为100
n_estimators = 1000
# 创建并训练梯度提升树回归模型
gb_model = GradientBoostingRegressor(n_estimators=n_estimators, random_state=42)
gb_model.fit(X_train, y_train)
# 在测试集上评估模型性能
gb_score  = gb_model.score(X_test,y_test)
print(f'Model Score : {gb_score}')
feature_importances = gb_model.feature_importances_
print("特征重要性:", feature_importances)
# 计算每轮迭代的损失值
train_loss = np.zeros((gb_model.n_estimators,), dtype=np.float64)
for i, y_pred in enumerate(gb_model.staged_predict(X_train)):
    train_loss[i] = mean_squared_error(y_train, y_pred)

test_loss = np.zeros((gb_model.n_estimators,), dtype=np.float64)
for i, y_pred in enumerate(gb_model.staged_predict(X_test)):
    test_loss[i] = mean_squared_error(y_test, y_pred)
# 在测试集上进行预测
y_pred = gb_model.predict(X_test)
theta = cosine(y_pred,y_test)
print(f'余弦相似度为{theta}')
# 计算均方误差
mse = mean_squared_error(y_test, y_pred)
print(f'均方误差为{mse}')
# 绘制损失曲线
plt.figure(figsize=(10, 6))
plt.plot(np.arange(1, gb_model.n_estimators + 1), train_loss, label='Training Loss')
plt.plot(np.arange(1, gb_model.n_estimators + 1), test_loss, label='Test Loss')
plt.xlabel('Number of Iterations')
plt.ylabel('Mean Squared Error')
plt.title('Training and Test Loss')
plt.legend()
plt.show()

Model Score : 0.18546794557325652
特征重要性: [0.20028747 0.20177494 0.20162684 0.20300577 0.19330498]
余弦相似度为0.9958717349992869
均方误差为0.002117230287145886

# 决策树

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score
from sklearn.tree import plot_tree
from sklearn.metrics import mean_squared_error
# 创建并训练决策树回归模型
tree_model = DecisionTreeRegressor(random_state=42)
tree_model.fit(X_train, y_train)
# 在测试集上评估模型性能
tree_score = tree_model.score(X_test, y_test)
print("Model Score :", tree_score)
# 输出特征重要性
feature_importances = tree_model.feature_importances_
print("特征重要性:", feature_importances)
# 将特征重要性展平
feature_importances_flat = np.ravel(feature_importances)
# 可视化特征重要性
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances_flat)), feature_importances_flat, tick_label=[f'{X.columns[i]}' for i in range(0, len(feature_importances_flat))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()
# 在测试集上进行预测
y_pred = tree_model.predict(X_test)
theta = cosine(y_pred,y_test)
print(f'余弦相似度为{theta}')
mse = mean_squared_error(y_test, y_pred)
print(f'均方误差为{mse}')

Model Score : 0.06121741625318666
特征重要性: [0.19938786 0.20108056 0.19957446 0.20072828 0.19922885]

余弦相似度为0.9952413117618393
均方误差为0.0024401972992366535

# 随机森林

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

rf = RandomForestRegressor(n_estimators=4, random_state=42) 
rf.fit(X_train, y_train)
y_pred = rf.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = rf.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')
feature_importances = rf.feature_importances_
print("特征重要性:", feature_importances)
feature_importances_flat = np.ravel(feature_importances)
# 可视化特征重要性
plt.figure(figsize=(30, 18))
plt.bar(range(len(feature_importances_flat)), feature_importances_flat, tick_label=[f'{X_train.columns[i]}' for i in range(0, len(feature_importances_flat))])
plt.xlabel('Feature')
plt.ylabel('Importance')
plt.title('Feature Importance')
plt.show()

Model score: 0.08951717139739401
均方误差为: 0.0023666371509472196
余弦相似度为0.9953848954929839
特征重要性: [0.19743656 0.2002398  0.20023184 0.20070336 0.20138843]

# Problem FOUR 预测洪水发生的概率

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import BaggingClassifier
from sklearn.ensemble import BaggingRegressor
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_iris
import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import joblib
class bagging_model:
    a_gb=0.15
    a_ju=0.15
    a_line=0.4
    a_rf=0.2
    def __init__(self) -> None:
        self.gb_tree=GradientBoostingRegressor(n_estimators=1000, random_state=42)
        self.ju_tree=DecisionTreeRegressor(random_state=42)
        self.normal_line=LinearRegression()
        self.rf_forest=RandomForestRegressor(n_estimators=20)
    def step(self,X_train,y_train,X_test,y_test):
        y_test=np.array(y_test)
        y_test=y_test.flatten()
        pre_gb=self.gb_tree.predict(X_test)
        pre_gb=pre_gb.flatten()
        pre_ju=self.ju_tree.predict(X_test)
        pre_ju=pre_ju.flatten()
        pre_line=self.normal_line.predict(X_test)
        pre_line=pre_line.flatten()
        pre_rf=self.rf_forest.predict(X_test)
        pre_rf=pre_rf.flatten()
        df=pd.DataFrame({
            'gb':pre_gb,
            'ju':pre_ju,
            'line':pre_line,
            'rf':pre_rf,
            'target':y_test
        })
        X_t=df.drop(columns='target')
        y_t=df['target']
        X_tr, X_te, y_tr, y_te = train_test_split(X_t, y_t, test_size=0.2, random_state=42)
        MODEL=LinearRegression(fit_intercept=False)
        MODEL.fit(X_tr,y_tr)
        temp=MODEL.coef_
        temp=temp/sum(temp)
        self.a_gb=temp[0]
        self.a_ju=temp[1]
        self.a_line=temp[2]
        self.a_rf=temp[3]
    def fit(self,X_train,y_train):
        X_t=X_train
        y_t=y_train
        # 有放回地抽样4次
        df = pd.merge(X_train,y_train,left_index=True,right_index=True)
        data_list=[df.sample(frac=1, replace=True) for i in range(4)]
        data=data_list[0]
        X_train=data.drop(columns='洪水概率')
        y_train=data['洪水概率']
        self.gb_tree.fit(X_train,y_train)
        data=data_list[1]
        X_train=data.drop(columns='洪水概率')
        y_train=data['洪水概率']
        self.ju_tree.fit(X_train,y_train)
        data=data_list[2]
        X_train=data.drop(columns='洪水概率')
        y_train=data[['洪水概率']]
        self.normal_line.fit(X_train,y_train)
        data=data_list[3]
        X_train=data.drop(columns='洪水概率')
        y_train=data['洪水概率']
        self.rf_forest.fit(X_train,y_train)
        self.step(X_t,y_t,X_test,y_test)
    def predict(self,X_test):
        pre_gb=self.gb_tree.predict(X_test)
        pre_gb=pre_gb.flatten()
        pre_ju=self.ju_tree.predict(X_test)
        pre_ju=pre_ju.flatten()
        pre_line=self.normal_line.predict(X_test)
        pre_line=pre_line.flatten()
        pre_rf=self.rf_forest.predict(X_test)
        pre_rf=pre_rf.flatten()
        y_pred=pre_gb*self.a_gb+pre_ju*self.a_ju+pre_line*self.a_line+pre_rf*self.a_rf
        return y_pred
    def score(self,X_test,y_test):
        score_gb=self.gb_tree.score(X_test,y_test)
        score_ju=self.ju_tree.score(X_test,y_test)
        score_line=self.normal_line.score(X_test,y_test)
        score_rf=self.rf_forest.score(X_test,y_test)
        score = score_gb*self.a_gb+score_ju*self.a_ju+score_line*self.a_line+score_rf*self.a_rf
        return score

# 训练集成模型

1
2
3
4
5
6
7
8

train_data_path='../Data/train.csv'
train_data = pd.read_csv(train_data_path)
train_data = train_data.drop(columns='id')
X=train_data.drop(columns='洪水概率')
y=train_data['洪水概率']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
model=bagging_model()
model.fit(X_train,y_train)

# 集成模型输出

1
2
3
4
5
6
7

y_pred=model.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
score = model.score(X_test,y_test)
print(f'Model score: {score}')
print(f"均方误差为: {mse}")
accuracy = cosine(y_test,y_pred)
print(f'余弦相似度为{accuracy}')

Model score: 0.182629512994938
均方误差为: 0.0021174324247682824
余弦相似度为0.9958713234970328

1
2
3
4
5
6
7
8
9
10
11
12
13

# 读取测试数据并进行预测
test_data_path = '../Data/test.csv'
test_data = pd.read_csv(test_data_path)
test_ids = test_data['id']
X=test_data.drop(columns='id')
X_test=X.drop(columns='洪水概率')

predictions = model.predict(X_test)

# 将预测结果写入submit.csv
submit_data = pd.DataFrame({'id': test_ids, '洪水概率': predictions})
submit_data.to_csv('../Data/submit.csv', index=False)
print("预测结果已保存到 submit.csv")

# 通过绘制图形（直方图和 QQ 图）、统计检验（Shapiro-Wilk 检验、Kolmogorov-Smirnov 检验、Anderson-Darling 检验、D'Agostino's K^2 检验）、描述性统计（偏度（Skewness）和峰度（Kurtosis））

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

import pandas as pd
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns

# 加载数据
data = pd.read_csv('/mnt/data/submit.csv')

# 提取洪水概率
flood_probabilities = data['洪水概率']

# 绘制直方图和Q-Q图
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
sns.histplot(flood_probabilities, kde=True)
plt.title('洪水概率的直方图')

plt.subplot(1, 2, 2)
stats.probplot(flood_probabilities, dist="norm", plot=plt)
plt.title('洪水概率的Q-Q图')

plt.show()

# Shapiro-Wilk检验
shapiro_test = stats.shapiro(flood_probabilities)
print(f"Shapiro-Wilk 检验: W={shapiro_test[0]}, p值={shapiro_test[1]}")
if shapiro_test[1] > 0.05:
    print("Shapiro-Wilk 检验结果: 数据服从正态分布")
else:
    print("Shapiro-Wilk 检验结果: 数据不服从正态分布")

# Kolmogorov-Smirnov检验
ks_test = stats.kstest(flood_probabilities, 'norm', args=(np.mean(flood_probabilities), np.std(flood_probabilities)))
print(f"Kolmogorov-Smirnov 检验: D={ks_test[0]}, p值={ks_test[1]}")
if ks_test[1] > 0.05:
    print("Kolmogorov-Smirnov 检验结果: 数据服从正态分布")
else:
    print("Kolmogorov-Smirnov 检验结果: 数据不服从正态分布")

# Anderson-Darling检验
ad_test = stats.anderson(flood_probabilities, dist='norm')
print(f"Anderson-Darling 检验: A2={ad_test.statistic}, 临界值={ad_test.critical_values}, 显著性水平={ad_test.significance_level}")
if ad_test.statistic < ad_test.critical_values[2]:  # Typically, the 5% significance level is used
    print("Anderson-Darling 检验结果: 数据服从正态分布")
else:
    print("Anderson-Darling 检验结果: 数据不服从正态分布")

# 偏度和峰度
skewness = stats.skew(flood_probabilities)
kurtosis = stats.kurtosis(flood_probabilities)
print(f"偏度: {skewness}")
print(f"峰度: {kurtosis}")
if abs(skewness) < 0.5 and abs(kurtosis) < 3:
    print("偏度和峰度结果: 数据服从正态分布")
else:
    print("偏度和峰度结果: 数据不服从正态分布")

# 显示Anderson-Darling检验结果
ad_test_results = pd.DataFrame({
    '临界值': ad_test.critical_values,
    '显著性水平(%)': ad_test.significance_level
})

print(ad_test_results)

# 绘制洪水概率的折线图
plt.figure(figsize=(10, 6))
plt.plot(data['id'], flood_probabilities, marker='o')
plt.title('洪水概率的折线图')
plt.xlabel('ID')
plt.ylabel('洪水概率')
plt.grid(True)
plt.show()

# 洪水灾害的数据分析与预测：基于多模型集成学习的方法研究

# 摘要

本研究围绕洪水灾害的分析与预测展开，利用多模型集成学习方法对洪水发生概率及其影响因素进行了系统性和科学性的探讨。研究通过对附件提供的 train.csv 数据进行预处理，确保数据的准确性和完整性，并基于多种相关性分析和预测模型，提出了有效的洪水预测和防洪减灾策略。

针对问题一，本文采用皮尔逊相关系数和斯皮尔曼等级相关系数对 20 个指标与洪水发生概率之间的相关性进行了详细量化分析。结果显示，基础设施恶化、季风强度、大坝质量、地形排水、河流管理等指标与洪水发生概率具有显著相关性。基于这些分析，研究提出了一系列合理的防洪措施建议，如加强基础设施建设、优化河流管理、提升大坝质量等，以有效降低洪水发生的风险。

针对问题二，本文利用 K-means 聚类方法将洪水发生概率划分为高、中、低三个风险等级，并深入分析了不同风险等级下的指标特征差异。通过对高、中、低风险类别中各项影响因素的详细探讨，建立了四个具有高度准确性的预警评价模型，明确了不通类别洪水风险的主要特征，为洪水风险管理提供了科学依据。

针对问题三，研究基于前述相关性分析结果，分别构建了多元线性回归模型、梯度提升树模型、决策树模型、随机森林模型和逻辑回归模型，系统评估了各模型的预测性能。为了提高预测的准确性和稳定性，本文通过 Bagging 技术整合多种模型的预测结果，提出了一个集成学习模型。该模型在对 test.csv 数据进行验证时，显示出优异的预测性能，进一步提升了洪水发生概率预测的准确性。

针对问题四，本文利用集成学习模型对 test.csv 中的数据进行了洪水发生概率的预测，并将预测结果填入 submit.csv 中。通过绘制直方图和折线图，对预测结果的分布进行了统计分析，并通过 Quantile-Quantile 图、Shapiro-Wilk 检验、Kolmogorov-Smirnov 测试和 Anderson-Darling 检验，验证了预测结果的分布特性。结果表明，模型预测的洪水发生概率分布较为合理，但不完全符合正态分布，具有一定的重尾现象，进一步验证了模型的有效性和可靠性。

总体而言，本研究不仅为洪水灾害的预测提供了有效的数学模型和工具，也为制定科学合理的防洪减灾策略提供了重要参考和借鉴。通过对数据的深入分析和模型的多样化应用，本文为相关领域的研究和实践提供了新的思路和方法。本文的研究成果不仅具有重要的理论价值，也在实际应用中具有广泛的推广潜力。

* 关键词 *：洪水灾害，数据分析，多模型集成学习，皮尔逊相关系数，斯皮尔曼等级相关系数，K-means 聚类，风险分类，Bagging 技术，预警评价模型，正态分布分析，防洪减灾措施

# 目录

* 洪水灾害的数据分析与预测：基于多模型集成学习的方法研究 *

* 一、问题背景与问题重述 *

1.1 问题背景

1.2 已知条件

1.3 问题重述

* 二、 模型假设 *

* 三、 符号说明 *

* 四、 研究思路 *

* 五、* * 问题一的建模与求解 *

5.1 数据预处理

5.2 问题一的模型建立与求解

5.2.1 用皮尔逊相关系数和斯皮尔曼等级相关系数研究各指标与洪水发生概率的相关性。

5.2.2 通过多元线性回归、递归提升树、随机森林对指标的优劣进行分析。

5.2.3 结果分析

5.3 建议与措施

* 六、问题二的建模与求解 *

6.1 问题 2 的分析

6.1.1 数据预处理

6.1.2 聚类方法选择

6.2 k-means 建模实现与分析

6.2.1 模型建立

6.2.2 高中低风险特征分析

6.2.3 建立发生洪水不同风险预警模型

6.3 模型灵敏度分析

6.3.1 基于特征变化率的灵敏度计算

6.3.2 基于预测结果熵值的灵敏度计算

七、 问题三的建模与求解

7.1 指标选取

7.1.2 基于过滤法选取关键指标

7.2 各模型结果分析

7.3 最终模型建立

* 八、问题四的建模与求解 *

8.1 结果预测

8.2 结果可视化

8.3 判断洪水发生概率是否符合正态分布

8.3.2 Shapiro-Wilk 检验

8.3.3 Kolmogorov-Smirnov 测试

8.3.4 Anderson-Darling 检验

* 九、** 模型分析 *

9.1 模型的优点

9.2 模型的缺点

9.3 模型的推广

* 十、** 总结 *

* 十一、** 参考文献 *

* 十二、** 附录 *

#

#

**
**

# * 一、问题背景与问题重述 *

# *1.1 问题背景 *

洪水作为一种由暴雨、急剧融冰化雪、风暴潮等自然因素引发的自然灾害，对全球范围内的社会经济系统构成了重大威胁。从先秦时期《尚书・尧典》到现代全球频发的洪水事件，无不彰显着洪水灾害的广泛性与严重性。进入 21 世纪以后，随气候环境的严重变化和人类活动加剧，洪水灾害的频率和强度显著提高，对人类社会造成了前所未有的挑战，本论文基于此，以附件里的数据为参考，对于洪水发生概率及其影响因素进行了深入的分析。

# *1.2 已知条件 *

​ 为了研究与洪水发生概率有关的各个指标，本题提供了 train.csv、test.csv、submit.csv 三个附件为研究依据。

1.1048574 个包含洪水事件 id、季风强度、地形排水、河流管理、森林砍伐、城市化、气候变化、大坝质量、淤积、农业实践、侵蚀、无效防灾、排水系统、海岸脆弱性、滑坡、流域、基础设施恶化、人口得分、湿地损失、规划不足、政策因素和发生洪水的概率的事件训练表格。（train.csv）

2.745305 个包含洪水事件 id 和上述 20 个指标得分，缺少洪水概率的测试表格。(test.csv)

3. 包含 test.csv 的洪水事件 id，缺少洪水概率的事件提交表格。(submit.csv)

为了更好的对洪水预防提供建议，本小组选取以下文章作为参考文献：

 1.选取对于河流管理因素治理的参考[1]

 2.选取对于水利工程等基础设施因素治理的参考[2][3]

 3.选取对于气候变化因素治理的参考[4][5]

 4.选取对于降水预测治理的参考[6]

# *1.3 问题重述 *

​ 问题 1 通过分析 train.csv 中的数据，分析并可视化上述 20 个指标与洪水发生的关联的强弱性。并根据可能的原因针对洪水预防提出建议和措施。

​ 问题 2 将附件 train.csv 中洪水发生的概率聚类成不同的类别，分析具有高中低风险的洪水事件的指标特征。

​ 问题 3 基于问题 1 中的指标分析结果建立洪水事件发生的概率预测模型。从 20 个指标中选取合适指标，预测洪水发生的概率，并验证模型准确性。并且尝试分析如果只用五个指标应该如何改进模型。

​ 问题 4 基于问题 3 建立的洪水发生概率预测模型，预测 test.csv 中事件发生洪水概率并填入 submit.csv 中，并绘制发生洪水概率的直方图和折线图，分析结果是否呈现正态分布。

# 二、 * 模型假设 *

\1. 假设题目所给的数据真实可靠；

\2. 假设洪水事件由且仅由这 20 个指标影响；

\3. 假设洪水发生概率与其他事件呈线性关系。

# 三、 * 符号说明 *

# 四、 * 研究思路 *


# 五、 * 问题一的建模与求解 *

在分析该题目的时候，本文首先使用相关系数判定各数据的相关性，然后使用多元回归算法建立这 20 个指标关于洪水概率的数学模型，最后使用梯度提升树、随机森林回归模型预测洪水发生的概率，并进行了 K 折交叉验证、模型评估、特征重要性分析以及损失曲线的绘制。

# *5.1 数据预处理 *

（1）针对给出的超过 100 万条的洪水事件数据，本文对 train.csv 进行基本的数据清洗，处理缺失值、异常值。

（2）由于 id 列对于数据分析无直接用处，本文对该列进行清洗处理。

# *5.2 问题一的模型建立与求解 *

# *5.2.**1 用皮尔逊相关系数和斯皮尔曼等级相关系数研究各指标与洪水发生概率的相关性。*

​ 首先，我们根据式一计算各个指标对于洪水概率的皮尔逊相关系数：

 （1）

​ 斯皮尔曼等级相关系数和！[img](file:///C:\Users\JackDu\AppData\Local\Temp\ksohtml23228\wps106.jpg) 的计算公式如式二三所示：

 （2）

 （3）

本文使用.core () 函数，通过皮尔逊相关系数公式和斯皮尔曼等级相关系数来对各个指标和洪水概率的相关系数进行计算，最终得到相关系数表格一：

表 1 相关系数前 10 的指标

最终可视化图片为：


图 1 不同因素对洪水概率的影响 图 2 不同因素与洪水概率的斯皮尔曼相关系数

为了便于分析相关性差异，我们将皮尔逊相关系数减去 0.175、斯皮尔曼等级相关系数减去 0.17 后再次可视化，得到结果为：


图 3 调整后的相关系数柱状图

# *5.2.**2 通过多元线性回归、递归提升树、随机森林对指标的优劣进行分析。*

先在模型一中，通过多元线性回归拟合模型里各指标的系数，通过计算预测值！[img](file:///C:\Users\JackDu\AppData\Local\Temp\ksohtml23228\wps113.jpg) 和真实值！[img](file:///C:\Users\JackDu\AppData\Local\Temp\ksohtml23228\wps114.jpg) 之间的均方误差！[img](file:///C:\Users\JackDu\AppData\Local\Temp\ksohtml23228\wps115.jpg)、在测试集的分数！[img](file:///C:\Users\JackDu\AppData\Local\Temp\ksohtml23228\wps116.jpg) 以及余弦相似度用以评估模型性能，表达式如下：

 (4)

 (5)

 (6)

# *5.2.**2.1 线性回归模型 *

多元线性回归模型是一种常用的统计分析方法，该模型表达式为：

 (7)

​ 首先初始化回归系数，然后使用梯度下降或正规方程法来迭代更新回归系数，最小化损失函数（均方误差），训练完成后，输出不同特征的回归系数，这些系数反映了每个特征对预测目标的影响程度。

使用训练数据（X_train 和 y_train）训练线性回归模型，并使用模型对测试集 X_test 进行预测，返回预测值 y_pred，并且得到性能指标如下：

表 2 线性回归模型性能指标

最终得到公式：

(8)

# *5.2.**2.2 梯度提升树模型 *

本文使用梯度提升树（Gradient Boosting Trees, GBT）模型来预测洪水发生概率，并通过一系列步骤来评估模型的性能、特征重要性以及训练过程中的损失变化。该模型通过组合多个弱学习器（决策树）形成一个强学习器，算法通过递归逐步提升优化模型的预测能力，在每一轮迭代中，计算当前模型的残差并基于该残差训练新的决策树，而新训练的决策树会被加入到模型中，与之前的模型加权融合以提高预测能力。

本文设置迭代次数为 500，即构建 500 棵决策树拟合数据。并使用 5 折交叉验证帮助评估模型在不同数据子集上的性能，估计模型泛化能力，其输出的每折交叉评分为 R2 的值。本文使用完整的训练集（X_train 和 y_train）训练梯度提升树模型，并在测试集上评估模型性能。同时输出特征重要性并使用 matplotlib 库绘制特征重要性的条形图，以直观地展示哪些特征对模型预测的贡献最大。

本文使用 staged_predict 方法获取模型训练集和测试集上每轮迭代后的预测结果，计算并存储每轮迭代后的均方误差（MSE），以观察训练集和测试集上损失值的变化情况，进而使用 matplotlib 库绘制训练集和测试集上的损失值（MSE）随迭代次数变化的曲线，帮助观察模型是否过拟合（即训练集上的损失持续下降，而测试集上的损失在某个点后开始上升）。

以下为该模型的性能指标：

表 3 梯度提升树模型性能指标

最终可视化结果如下图：


图 4 梯度提升树各指标的重要性

# *5.2.**2.3 随机森林模型 *

随机森林模型是多个决策树的集成，通过对每棵树的预测结果进行投票或加权平均来得到最终的预测结果。随机森林通过随机选择样本和特征来构建每棵树，从而增加模型的鲁棒性和泛化能力。整个过程线通过自助法（Bootstrap）从训练数据中有放回地抽样，生成多个不同的子集，然后对每个子集，随机选择部分特征，训练一棵决策树，最后将所有树的预测结果进行投票（分类问题）或加权平均（回归问题），得到最终的预测结果，从而输出不同指标的权重（特征重要性），表示每个特征在所有树中的分割点选择过程中所贡献的纯净度提升的总和。

本文设置随机森林的决策树数量为 200，通过抽取部分训练数据训练模型并使用投票策略以加快模型训练速度并提高模型精度。该模型的性能指标如表 4 所示：

表 4 随机森林模型性能指标

随机森林得到的特征重要性如下图：


图 5 随机森林各指标的重要性

# *5.2.3* * 结果分析 *

根据皮尔逊相关系数和斯皮尔曼等级相关系数结果分析，可以发现 **** 基础设施恶化、季风强度、地形排水 * 三个指标与洪水发生概率的相关性最强，而 * 海岸脆弱性、排水系统、侵蚀 **** 三个指标与洪水发生概率的相关性最弱。

针对指标优劣性的分析，本文采用了线性回归模型、递归提升树回归模型和随机森林模型。本文对三种模型得出的指标系数进行降序排列，并记前十项指标得分为 1，后十项指标得分为 - 1，将三种模型的指标得分进行累加，通过比较各指标得分情况判断指标的重要性。最终得到各指标得分条形图。


图 6 各指标得分条形图

由图可以发现，经过三个模型的综合分析，森林砍伐与洪水发生概率的相关性最强，而淤积、滑坡、城市化、海岸脆弱性、气候变化、侵蚀、人口得分、大坝质量、河流管理、地形排水十个指标与洪水发生概率的相关性较强，农业实践、排水系统、规划不足、季风强度、政策因素、流域、基础设施恶化七个指标与洪水发生概率的相关性较弱 **，湿地损失、无效防灾 ** 两个指标与洪水发生概率的相关性最弱。

尽管三个模型得出的指标相关性排序结果各有差异，但观察数据可以发现各指标重要性数值差异并不大，因此在考虑洪水发生概率及应对洪水灾害时，* 以上 20 个指标都需要考虑 *。

综上，本文认为 **** 森林砍伐、淤积、气候变化、人口得分、大坝质量、河流管理、地形排水 * 与洪水发生概率的相关性 * 最大 *，而 * 无效防灾、湿地损失 * 与洪水发生概率的相关性 * 不大 ****。

对于各指标与洪水发生的相关性的原因，本文作出以下分析：森林通过树木和植被可以有效保持土壤，减少地表径流，并且树木的根系可以固土。而当森林遭到砍伐且暴雨发生时，土壤失去固定易被冲刷，导致泥沙进入河流，削弱河道蓄水能力，导致水量激增进而增加洪水发生的概率。同时砍伐森林会降低蒸腾作用，导致削弱大量水分返回大气的效果，进而导致地表径流增加，增加洪水发生概率。淤积会减少河流、湖泊和水库的蓄水容量，同时改变河道的水流速度和流向，增加洪水发生的复杂性和频率，因此降水量较大时，水体更易溢出河道，导致洪水发生。针对气候变化，现世界范围气候变化多端，导致极端天气事件（如暴雨、飓风）频率增加，降水在短时间内过于集中，导致洪水灾害发生。同时，人口增长和城市化导致土地利用变化，如水泥地面增加减少地表渗透，加剧地表径流，且人口密度高的地区，排水系统和其他基础设施在暴雨时容易超负荷，进而增加洪水发生的概率。大坝在控制水流和调节河流水位方面至关重要，若大坝质量不佳，可能在洪水期间受损坍塌，导致严重泄洪。因此良好的河流管理（如疏浚河道、修建防洪堤坝）可以有效地预防和控制洪水，若河流管理不恰当，即难以预测灾害发生。良好的地形排水设计能够快速排除积水，进而减少洪水风险，且低洼地区和排水不通畅的区域更容易发生洪水。因此地形排水与洪水发生概率的相关性也较大。[2]

针对低相关性的指标，虽然无效的防灾措施可能会增加洪水的破坏程度，但它们并不会直接增加洪水发生的概率，无效防灾主要影响的是洪水造成的损失，而不是洪水发生的本身。而湿地确实具有调节水流、减缓洪水的功能，但其影响相对于森林和河流管理可能较小，湿地损失更多影响的是水生态系统的健康和水质，而不是洪水发生的直接原因。

# *5.3** 建议与措施 *

结合相关文献，本文针对以下方面给出提前预防洪水的建议和措施：

（1）加强森林保护和恢复：实施大规模的植树造林项目，恢复被破坏的森林，扩大森林范围，提高土壤的吸水和固土能力。同时制定并严格执行森林保护法律法规，防止非法乱砍和乱伐现象。

（2）改善河流管理和治理：定期清理河道中的淤积物，保持河道的蓄水和排水能力 [1]。同时需要修建防洪堤坝、水闸等设施，并定期检查和维护，确保其在涨洪期间可以正常运行 [2]。

（3）提高大坝质量和管理水平：先需要对现有的大坝进行定期检测和维护，确保其结构安全。同时采用先进的技术和材料，科学设计和建造新大坝，提高大坝的防洪能力 [2]。

（4）应对气候变化：利用气象监测和数据分析技术，建立暴雨早期预警系统，及时发布洪水预警信息，提前做好预防。需要加强对气候变化对区域降水模式和洪水风险影响的研究，为防洪决策提供科学依据 [3]。

（5）改善城市和农村排水系统：升级和改造城市和农村的排水系统，确保在暴雨期间能迅速排除积水 [4]。除此之外，还需要定期清理排水管道和沟渠，防止淤积和堵塞。

（6）加强公众教育和应急准备：除了进行生态治理，还需要通过社区活动和媒体宣传，提高公众的防灾意识和应急响应能力。同时制定详细的洪水应急预案，并定期进行演练，确保居民在洪水来临时能够迅速、安全撤离。

# * 六、问题二的建模与求解 *

# *6.1* * 问题 2 的分析 *

问题 2 要求将附件 train.csv 中的洪水发生概率进行聚类，并分析具有高、中、低风险的洪水事件的指标特征。然后选取合适的指标，计算不同指标的权重，建立不同风险的预警评价模型，最后进行模型的灵敏度分析。

# *6.1.1* * 数据预处理 *

在进行聚类分析之前，首先需要对数据进行预处理：首先检查数据中是否存在缺失值，并选择合适的方法进行填补。例如，可以使用均值填充法或者插值法填补缺失值。

# *6.1.2* * 聚类方法选择 *

选用 K 均值聚类（K-means clustering）方法对洪水发生概率进行聚类。K 均值聚类是一种非监督学习方法，通过迭代不断调整聚类中心的位置，使得每个数据点到其最近聚类中心的距离之和最小。本文使用该算法将洪水发生概率分为高、中、低三个风险等级。该算法通过迭代的方式寻找 k 个簇的划分方案，使得聚类结果对应的代价函数最小，代价函数通常定义为各个样本点到其所属簇中心点的误差平方和。该算法利用 N 维空间距离对初始聚类中心进行学习，得到最终聚类中心，步骤如下：

\1. 首先选择要将数据集分成 3 个簇，然后随机选择 3 个数据点作为初始簇中心（即为高、中、低三个风险等级）。

\2. 将每个数据点分配到距离其最近的簇中心。

\3. 根据分配的数据点更新簇中心点，即通过计算属于每个簇的数据点的平均值实现。

\4. 重复步骤 2 和 3 直至簇中心点不再发生变化。

\5. 输出 3 个簇和每个簇的中心点。

# *6.2* *k-means 建模实现与分析 *

# *6.2.1* * 模型建立 *

首先，我们以公式（9）计算数据点到聚类中心的距离，然后对于每个数据点，用公式（10）找到最近的聚类中心。最后重新计算每个聚类中心的位置，使其等于簇内所有数据点的均值。

 (9)

 (10)

 (11)

​ 重复上述步骤，直到聚类中心不再发生改变，我们可以得到各风险层级聚类中心表格。

表 5 各风险层级聚类中心

洪水发生概率更靠近哪个聚类中心则被划分为哪个聚类。例如，当洪水发生概率为 0.55 时为高风险，洪水发生概率为 0.51 时为中风险，洪水发生概率为 0.42 时为低风险。

​ 我们将区分开后的 train.csv 分别保存为高风险.csv、中风险.csv 和低风险.csv 以便于后续特征分析。

# *6.2.2* * 高中低风险特征分析 *

​ 对于不同风险等级的数据，我们首先分别计算其 20 个指标的平均值和方差。然后通过式（12）进行标准化数据，即将不同量纲的数据转换到一个相同的尺度上，使其均值为 0，标准差为 1。这样可以消除不同特征之间的量纲差异，使得它们在同一尺度上进行比较。

 (12)

​ 据此公式，将 2032 个指标的数据标准化和特征可视化于图像。


图 7 不同风险类别的标准化特征平均值条形图

​ 通过对图 7 的分析，我们可以知道高中险聚类中，大坝质量、河流管理两个指标的影响最大，而流域、人口得分、规划不足三个指标的影响最小。在中风险聚类中，河流管理、侵蚀、无效防灾、海洋脆弱性、湿地损失、规划不足六个指标的影响最大，而季风强度、基础设施恶化两个指标的影响最小。在低风险聚类中，河流管理、大坝质量两个指标的影响极为突出，而流域、人口得分、规划不足三个指标的影响最小。由此可以得出，洪水发生高风险地区主要受设施建设影响，而中风险地区主要受地质因素影响，低风险地区也主要受设施建设影响。


图 8 不同风险类别的标准化特征方差条形图

​ 通过对图 8 的分析，我们可以知道在三个聚类中，地形排水、城市化、海岸脆弱性、政策因素的不确定性和波动性都很高，而季风强度、森林砍伐、气候变化、基础设施恶化四个指标的不确定性和波动性较低，也可以看出其对于各聚类洪水发生事件都有影响，进而体现其相关性强。除以上指标外，在高风险聚类中城市化、大坝质量、流域等指标的不确定性和波动性较强，即其在高风险聚类中的相关性并不强。

​ 然后，我们利用差值图，通过计算高风险和低风险类别的特征差异，来直观展示哪些特征在不同风险类别中存在显著差异，为了突出差异，我们将最后的结果均减去 0.8，以此来筛选具有显著差异的因素。

 (13)

​ 

图 9 高风险与低风险类别特征平均值差异条形图

​ 最后我们使用热力图，通过热力图结合层次聚类，直观展示不同风险类别中各特征的表现，识别特征之间的相似性和差异性，帮助分析洪水风险的主要影响因素。对标准化后的数据进行层次聚类分析，将相似的特征聚合在一起。根据聚类结果对数据进行重新排序，以确保相似的特征相邻。使用热力图将重新排序后的数据可视化，用颜色表示数据值的大小。最终得到该不同风险与指标间的热力图。在该热力图中，颜色越淡则代表该指标与该聚类中的洪水发生相关性越弱，而颜色越浓则表示其相关性越强。因此，可以直观分析得到在 **** 高风险聚类 * 中，* 河流管理、大坝质量、湿地损失、基础设施恶化 * 四个指标的相关性最强，而 * 低风险聚类 * 中 * 地形排水、淤积、人口得分、流域、滑坡 * 等指标相关性最强，在 * 中风险聚类 * 则 * 季风强度、气候变化、基础设施恶化 **** 三个指标的相关性最强。通过以上分析也可以更好帮助我们在不同风险聚类更好制定预防措施并分析风险来源。


图 10 各指标层次聚类热力图

# *6.2.**3* * 建立发生洪水不同风险预警模型 *

​ 为了建立一个有效的洪水风险预警评价模型，我们分别尝试了选取 20、问题 1 中前 15 个、前 10 个相关性最大的指标进行建模。结果显示当选取 20 个指标的情况下模型拟合程度 **** 最高 ****。故我们选取全部指标进行模型的建立。

​ 为了构建洪水风险预警评价模型，我们选择了几种常见的机器学习模型，包括逻辑斯特回归、线性回归、梯度提升树、决策树和随机森林。这些模型在处理回归和分类问题时具有较好的表现，并且能够提供特征重要性分析。我们通过问题 1 中训练的模型，得出不同指标的权重，并且在本问题中增加了 **** 决策树 * 做法和逻辑斯特回归，使得预警模型更加的全面。并且我们定义了一个 * MyModel**** 类，用于调用预警模型，并且以公式（14）来获取模型预测正确率（Model Score）。

 (14)

​

​ 我们通过分析建立建立的五个模型 (线性回归:linear、梯度提升树：gb_tree、决策树：ju_tree、随机森林：rf_forest，逻辑斯特回归：logistic)，以 MyModel 调用模型得出结果，最终得到模型的性能可视化结果为：


图 11 各模型性能对比条形图

对比这三个图，我们可以发现，线性回归模型得到的模型性能最优。

# *6.3* * 模型灵敏度分析 *

灵敏度分析的目标是评估不同模型对输入特征变化的响应程度。这一过程有助于识别哪些特征对模型的预测结果影响最大。为了深入探究线性回归、梯度提升树、决策树和随机森林四个模型的灵敏度。我们采用两种方法进行灵敏度分析：基于特征变化率的灵敏度计算和基于预测结果熵值的灵敏度计算。

而逻辑斯特模型不能使用特征变化率的灵敏度计算，因为逻辑斯特模型的输出是通过对线性组合后的特征进行 Sigmoid 函数变换而得：

 (15)

 (16)

特征变化率假设特征与目标变量之间是线性关系，但逻辑斯蒂回归是非线性的。这使得特征值的简单变化并不能直接线性地反映在输出概率的变化上，并且 Sigmoid 函数具有非线性性质，在极端值处（接近 0 或 1）变化很小，而在中间值（接近 0.5）变化较大。这意味着特征值的相同变化在不同范围内对输出概率的影响不同，导致特征变化率无法简单地反映这种非线性变化。这种非线性会导致特征变化率难以稳定地捕捉特征变化对输出的影响。

# *6.3.1* * 基于特征变化率的灵敏度计算 *

​ 基于特征变化率的灵敏度计算通过评估输入特征变化对模型预测结果的影响，来衡量每个特征对模型预测的敏感度。对每个特征逐个进行变化，从最小值逐步增加到最大值，每次变化后记录模型的预测结果的平均值。然后利用式（17）计算变化率。最后计算所有变化率的平均值，作为该特征的灵敏度度量。变化率越高代表模型的灵敏度越高，效果越好。

 (17)

 (18)


​ 通过公式（18）计算得到四个模型的基于特征变化率的灵敏度。

​  

图 12 线性回归模型属性平均变化率 图 13 递归提升树模型属性平均变化率

 

图 14 决策树模型属性平均变化率 图 15 随机森林模型属性平均变化率

由上图可以看出，不同模型的特征变化率各有差异。图 12 表明在线性回归模型中，季风强度、地形排水、河流管理等特征的变化率较高，说明这些特征对模型的灵敏度较大。图 13 表明在梯度提升树模型中，季风强度和地形排水的变化率较高，其他特征的影响相对较小。图 14 表明决策树模型中，特征变化率分布较为均匀，无明显突出的特征。图 15 表明随机森林模型中整体变化率较低，说明模型对单个特征的依赖程度较低。

# *6.3.2* * 基于预测结果熵值的灵敏度计算 *

​ 基于预测结果熵值的灵敏度计算通过评估模型预测结果的分布不确定性，来衡量模型的灵敏度。熵值越高，表示模型的预测结果越不确定，灵敏度越高。我们运用公式（19）得到熵值：

 (19)


在图中，我们展示了五种不同模型的灵敏度（通过熵值衡量）。熵值越高，表示模型的预测结果越不确定，预测的稳定性越差；熵值越低，表示模型的预测结果越确定，模型的稳定性越高。以下是针对五个模型的具体分析：

1. 线性回归模型：该模型熵值最高，接近 1.6，表示该模型的预测结果具有极高的不确定性。尽管线性回归模型简单且易于理解，但在处理复杂性较高的非线性关系时，表现不如其他模型。且高熵值意味着线性回归模型在特征变化时，预测结果的波动较大，对输入数据的变化非常敏感。

2. 梯度提升树模型：该模型熵值较高但低于线性回归模型，这表明该模型在预测结果的确定性方面表现比线性回归模型稍好。梯度提升树模型通过集成多棵决策树的预测结果，可以捕捉到更复杂的特征关系，因此在稳定性上有所提升。

3. 决策树模型：该模型熵值接近梯度提升树模型，但略高。这表明单棵决策树模型在处理特征变化时的稳定性相对较差，预测结果波动较大。且该模型容易过拟合，特别是在训练数据较少的情况下，熵值较高恰好反映了这点。

4. 随机森林模型：该模型熵值接近 1.4，略低于梯度提升树和决策树模型。随机森林模型通过集成多棵决策树，增强了预测结果的稳定性。较低的熵值表明随机森林模型在特征变化时的预测结果相对稳定，具有较好的鲁棒性。

5. 逻辑回归模型：该模型熵值与随机森林模型相近，逻辑回归模型在处理分类问题时具有较好的表现，但在特征变化率计算上存在局限性。熵值反映了逻辑回归模型在预测结果上的确定性，表现出一定的稳定性。

综上，线性回归模型的预测结果不确定性最高，灵敏度也较高。梯度提升树和决策树模型的熵值稍低于线性回归模型，但仍较高，说明这些模型在处理特征变化时有较大波动。随机森林模型和逻辑回归模型的熵值较低，即这些模型的预测结果稳定性表现较好，灵敏度较低。整体而言，随机森林模型在灵敏度分析中表现最佳，具有较高的预测稳定性和确定性。因此，随机森林模型可能是处理洪水概率预测的最佳选择。

# 七、* 问题三的建模与求解 *

*7.1 指标选取 *
*7.1.1** 初步选取指标 *

本文基于问题 1 的结论，将 20 个指标分成三组逐步分析，最终选取拟合度最高的指标组。首先将 20 个指标全部进行考虑，再根据问题 1 的结论剔除相关性最弱的湿地损失、无效防灾两个指标，最后剔除湿地损失、无效防灾、农业实践、排水系统、规划不足、季风强度、政策因素、流域、基础设施恶化 9 个相关性不强的指标，保留剩余 11 个相关性较强的指标。本文通过以上三组指标选取重新训练已有模型，并对各模型计算测试集分数、均方误差、余弦相似度以评估模型性能。

# *7.1.2* * 基于过滤法选取关键指标 *

过滤法是一种特征选择的方法，用于选择对模型预测结果影响最大的特征。其主要优点是不依赖于特定的模型，通过评估各个特征与目标变量之间的统计关系来进行特征选择。该方法可以帮助我们仅从题中 20 个指标的角度挖掘分析其价值高低，从而实现特征排序和选择。
在本文中，我们基于上述方法，从 11 个相关性较强的指标中筛选出 5 个关键影响指标。这 11 个指标包括森林砍伐、淤积、滑坡、城市化、海岸脆弱性、气候变化、侵蚀、人口得分、大坝质量、河流管理、地形排水。由于森林砍伐与洪水发生概率的相关性最强，因此我们首先保留该指标，并从剩余 10 个指标中筛选出 4 个最具影响力的指标，最终选取森林砍伐、淤积、大坝质量、河流管理、地形排水这 5 个关键指标。

使用这些关键指标，我们对模型进行重新训练并进行评估，以验证所选指标的有效性。

# *7.2 各模型 ** 结果分析 *

以上 4 组指标选取后各模型评估结果如下表：

表 6 各指标分组下模型的性能评估

通过横向对比各个模型，我们可以发现预测的准确度随着指标的增加而增加，这表明我们需要选择所有指标作为参数才能得到最优的预测模型。而如果只能选择五个数据作为参数的话，经过实验表明，如果只能选取五个指标，那么通过滤波法选取的指标表现是最优的。

# *7.3 最终模型建立 *

在本文中，我们采用 Bagging (Bootstrap Aggregating) 技术来建立最终的集成学习模型。Bagging 是一种集成学习方法，通过结合多个模型的预测结果来提高模型的准确性和稳定性。

Bagging 技术的基本思想是通过对数据进行有放回抽样生成多个不同的训练子集，然后在每个子集上训练一个模型，最后将所有模型的预测结果进行结合（例如通过投票或平均）来得到最终的预测结果。Bagging 的优点在于可以有效减少模型的方差，提高预测的准确性和稳定性。

在实现集成模型的过程中，我们首先创建了一个包含多个模型并实现了 fit、predict 和 score 方法的类 BaggingModel。在 fit 方法中，我们通过对训练数据进行有放回抽样，生成多个数据子集，并分别在每个子集上训练一个模型。随后，在 step 方法中，我们调整各个模型的权重，使得它们的和为 1。在 predict 方法中，我们使用各个模型对测试数据进行预测，并根据调整后的权重计算最终的加权预测结果。

最终我们得到了训练的最优解模型 (二十个指标) 和如果只能选取五个参数后优化的集成预测模型。

# * 八、问题四的建模与求解 *

# *8.1 结果预测 *

我们采用 7.3 中建立的集成学习模型，通过对 test.csv 模型里参数的读取和预测，最终得到各个 id 的洪水发生概率，并且保存在 submit.csv 中。

# *8.2* * 结果可视化 *

本文使用 matplotlib 库绘制直方图和折线图。


图 16 洪水概率预测直方图 图 17 洪水概率预测 Q-Q 图！[img](file:///C:\Users\JackDu\AppData\Local\Temp\ksohtml23228\wps148.jpg)

图 18 洪水概率预测折线图

*8.3 判断洪水发生概率是否符合正态分布 *

满足正态分布的直方图呈现钟形，即中间高、两边低，且直方图左右对称，左右两侧的频数大致相等，正态分布图像只有一个峰值，数据集集中在某一个值附近。由此可见，可以通过直方图 16 初步判断洪水发生概率符合正态分布。

但仅凭直方图并不能得出最终结论，因此本文结合其他方法进行综合分析。*8.3.**1* *Quantile-Quantile 图 *

Q-Q 图将样本数据从小到大排列，根据个数 n 计算每个数据点的理论分位数。对于正态分布，该理论分位数根据正态分布的累积分布函数计算得到，第 i 个数据点对应分位数常用方法如式（20）。计算后在图上绘制实际数据与理论分布分位数的点。如图 17 绘制了洪水概率预测 Q-Q 图，数据点偏离了对角线，特别是尾部部分，由此说明数据并不完全服从正态分布且存在重尾现象。

 (20)

# *8.3.**2* *Shapiro-Wilk 检验 *

该方法将样本数据按从小到大排列，按排列后的顺序计算权重系数，这些系数依赖于样本大小 n 和正态分布的期望与方差，再由式（21）计算统计量 W 值，根据 W 值和样本大小 n 计算 p 值。本模型计算得到 p 值远小于 0.05，拒绝原假设，数据不服从正态分布。

 (21)

# *8.3.**3* *Kolmogorov-Smirnov 测试 *

首先从每个数据点中减去平均值并除以标准差，将数据转换为均值为 0，标准差为 1 的分布，进行数据标准化。再使用 Kolmogorov-Smirnov 测试，比较该样本与正态分布的相似性。最终由式（22）返回样本分布与参考分布之间的最大差异 D，并输出 p 值即表示观察到的数据或更极端的数据在零假设下出现的概率。该模型最终预测概率的 p 值为 0.0，小于显著性水平（0.05），故认为数据不服从正态分布。

 (22)

 (23)


# *8.3.**4* *Anderson-Darling 检验 *

该方法将样本数据从小到大排序，计算该理论分布在数据点处的累加分布函数，并计算统计量，如式（24）所示。计算所得 A2 值远大于所有显著水平对应的临界值，拒绝原假设，数据不服从正态分布。

(24)

综上，预测洪水发生概率不服从正态分布。

# * 九、** 模型分析 *

# *9.1* * 模型的优点 *

​ 本模型具有多个显著优点。首先，模型在数据预处理阶段对附件中的数据进行了全面的清理，确保了数据的准确性和完整性。这包括处理缺失值和异常值，从而提升了模型的可靠性。其次，本文通过详细的相关性分析，使用了皮尔逊系数和斯皮尔曼相关系数来评估各指标对洪水发生概率的影响，并可视化结果，使得相关性差异一目了然。此外，梯度提升树模型在拟合优度检验中表现出色，说明其能较为准确地预测各影响指标与洪水发生概率之间的关系。本文使用多个模型（如多元线性回归、梯度提升树、随机森林等）进行对比验证，确保了模型结果的准确性和可靠性。最后，本文在建立最终洪水发生概率预测模型时采用了 Bagging 技术，提高了模型的预测精度和稳定性。

# *9.2* * 模型的缺点 *

​ 尽管本模型在多个方面表现出色，但仍存在一些不足。首先，各模型训练后计算的各指标系数的相关性并不完全一致，表明在某些参数调整上还有提升空间。其次，尽管模型采用了多种技术手段进行优化，但在实际应用中，模型的泛化能力和鲁棒性仍需进一步验证。

# *9.3* * 模型的推广 *

​ 本模型具备广泛的应用前景。首先，本文所建立的模型经过充分训练，可以随时调用和调整选用的指标重新训练，以适应不同的实际需求。其次，针对洪水指标的影响分析和预防措施的提出，对政府、相关部门和公众采取防洪减灾措施具有积极指导意义。此外，本文使用的多元线性回归模型、梯度提升树模型、随机森林算法和 K-means 聚类算法适用于各领域的问题分析与预测，具有很强的通用性和借鉴意义。本文采用的模型性能评估指标和方法，如余弦相似度、均方误差、K 折交叉验证及测试集分数，可以有效评估其他模型的性能，值得推广使用。该模型可以用于洪水风险评估与管理，将洪水风险评估纳入城市发展计划，识别并减轻潜在的洪水危害，指导关键基础设施的建设和维护，以更好应对洪水灾害。利用该模型评估森林砍伐等指标对洪水风险的影响，可以规划再造林项目及其他相关措施以降低洪水风险。最后，该模型为分析其他自然灾害的影响因素提供了思路和指导，能够帮助预测和分析其他自然灾害的发生概率。

# * 十、** 总结 *

本文通过对已知条件中的数据进行分析，研究了与洪水发生概率相关的指标，并通过数据集训练模型来预测洪水发生概率。在模型的构建和分析过程中，本文综合运用了多种方法，如相关性分析、多元线性回归、梯度提升树、随机森林和 Bagging 技术等。通过对比各个模型的性能，发现梯度提升树模型在预测洪水发生概率方面表现最佳。此外，本文筛选出森林砍伐、淤积、大坝质量、河流管理、地形排水五个关键指标，并验证了这些指标在预测洪水发生概率中的重要性。本研究不仅为洪水预测提供了有力工具，还提出了切实可行的防洪措施，能够有效降低洪水灾害带来的损失，为相关领域的研究和实际应用提供了宝贵的参考和借鉴。

# * 十一、** 参考文献 *

[1] 耿晓君，杨晓茹，李爱花，等。现代化防洪减灾体系中河道治理思路探讨 [J]. 中国水利，2023 (14):11-14.

[2] 刘影。水利工程对洪水的影响研究 [J]. 水利技术监督，2016,24 (05):15-16.

[3] 冯燕。气候变化背景下塔里木河流域洪水灾害的经济社会影响及应对 [J]. 环境社会学，2023 (01):83-105+197.

[4] 原文林，杨逸凡，赵小棚，等。降雨监测与预报技术在防洪减灾中的应用进展 [J/OL]. 人民长江：1-8 [2024-07-06].http://kns.cnki.net/kcms/detail/42.1202.TV.20240607.1034.004.html.

[5] 杨羽菲。中山市防洪减灾思路与对策研究 [J]. 广东水利水电，2023 (10):60-64.

[6] 刘欣。水利部将全面增强海河流域防洪减灾能力 [N]. 法治日报，2023-08-25 (005).DOI:10.28241/n.cnki.nfzrb.2023.004825.