第十周周四

在搞数学建模竞赛，用python完成了数据分析。

import os
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np
import plotly.graph_objs as go
from concurrent.futures import ProcessPoolExecutor
from multiprocessing import freeze_support
from statsmodels.tsa.statespace.sarimax import SARIMAX
from datetime import datetime, timedelta

# 定义函数读取单个文件的数据
def read_file(file_path):
    return pd.read_csv(file_path, encoding='utf-8', delimiter=",")

def main():
    # 定义第一个文件夹路径
    folder_path_1 = r"C:\Users\86155\Desktop\建模竞赛\2024校赛题目\2024校赛题目--C\附件2：2014年01月-2024年03月全国主要城市空气质量数据\all\all1"

    # 初始化一个空 DataFrame
    all_data_1 = pd.DataFrame()

    # 定义并行处理的进程数
    num_processes = os.cpu_count()  # 使用CPU核心数作为进程数
    i = 0
    # 使用进程池并行读取所有文件的数据
    with ProcessPoolExecutor(max_workers=num_processes) as executor:
        # 提交所有文件读取任务
        futures = [executor.submit(read_file, os.path.join(folder_path_1, file_name)) for file_name in
                   os.listdir(folder_path_1) if file_name.endswith(".csv")]

        # 获取所有任务的结果
        for future in futures:
            i += 1
            print(i)
            data = future.result()
            all_data_1 = pd.concat([all_data_1, data], ignore_index=True)

    # 对第一个文件数据进行清洗和预处理（这里假设无需清洗和预处理）

    # 将日期列和小时列合并为日期时间类型
    all_data_1['datetime'] = pd.to_datetime(
        all_data_1['date'].astype(str) + ' ' + all_data_1['hour'].astype(str).str.zfill(2), format='%Y%m%d %H')
    # 删除原始的日期和小时列
    all_data_1.drop(columns=['date', 'hour'], inplace=True)

    # 用副本解决 'SettingWithCopyWarning' 问题
    shijiazhuang_data = all_data_1[['datetime', '石家庄', 'type']].copy()

    # 将日期设置为索引
    shijiazhuang_data.set_index('datetime', inplace=True)

    # 提取AQI指数的数据
    aqi_data = shijiazhuang_data[shijiazhuang_data['type'] == 'AQI']

    # 定义第二个文件路径
    file_path_2 = r"C:\Users\86155\Desktop\建模竞赛\2024校赛题目\2024校赛题目--C\附件1：2014年01月-2024年02月石家庄天气数据.xlsx"

    # 读取第二个文件的数据
    weather_factors = pd.read_excel(file_path_2)
    if weather_factors.empty:
        print("文件未打开")

    # 将日期列和小时列合并为日期时间类型
    weather_factors['datetime'] = pd.to_datetime(weather_factors['年'].astype(str) + '-' +
                                                 weather_factors['月'].astype(str).str.zfill(2) + '-' +
                                                 weather_factors['日'].astype(str).str.zfill(2) + ' ' +
                                                 weather_factors['时'].astype(str).str.zfill(2) + ':00')
    # 将缺失数据（9999）替换为 NaN
    weather_factors.replace(9999, pd.NA, inplace=True)

    # 合并两个数据集
    merged_data = pd.merge(aqi_data, weather_factors, on='datetime')

    # 划分特征和目标变量
    X = merged_data[['气温（因子10）℃', '露点温度（因子10）℃', '海平面气压（因子10）Pa', '风向（由北顺时针）°',
                     '风速（因子10）m/s', '天空状况总覆盖代码', '液体沉积深度（因子10，周期6小时）mm']]
    y = merged_data['石家庄']

    # 确保在划分训练集和测试集之前用平均值填充NaN
    X.fillna(X.mean(), inplace=True)
    y.fillna(y.mean(), inplace=True)

    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

    # 建立线性回归模型
    model = LinearRegression()

    # 在划分和训练模型之前解决 'Input contains NaN, infinity or a value too large for dtype' 问题
    X_train.fillna(X_train.mean(), inplace=True)
    y_train.fillna(y_train.mean(), inplace=True)
    X_train = X_train.astype(float).fillna(5)
    y_train = y_train.astype(float).fillna(5)
    # 在处理过的训练集上训练模型
    model.fit(X_train, y_train)

    # 获取特征的系数
    feature_coef = pd.Series(model.coef_, index=X.columns).abs().sort_values(ascending=False)

    # 打印出对 AQI 影响最大的因素
    print("对 AQI 影响最大的因素:", feature_coef.index[0])

    X_test = X_test.astype(float).fillna(5)
    # 在测试集上进行预测
    y_pred = model.predict(X_test)

    # 计算均方误差
    mse = mean_squared_error(y_test, y_pred)
    print("均方误差:", mse)

    # 计算均方根误差（RMSE）
    rmse = np.sqrt(mse)
    print("均方根误差(RMSE):", rmse)

    # 计算决定系数（R^2）
    r2 = r2_score(y_test, y_pred)
    print("决定系数(R^2):", r2)

    # 创建交互式图表
    fig = go.Figure()

    # 添加原始数据曲线
    fig.add_trace(go.Scatter(x=merged_data['datetime'], y=merged_data['石家庄'], mode='lines', name='石家庄原始AQI数据'))

    # 处理预测特征集前的缺失值
    weather_factors.fillna(method='ffill', inplace=True)  # 以前向填充方式处理缺失值

    # 确保特征匹配（假设这些是训练模型时使用的特征）
    selected_features = ['气温（因子10）℃', '露点温度（因子10）℃', '海平面气压（因子10）Pa',
                         '风向（由北顺时针）°', '风速（因子10）m/s', '天空状况总覆盖代码',
                         '液体沉积深度（因子10，周期6小时）mm']
    weather_factors_for_prediction = weather_factors[selected_features]

    # 再次确认无 NaN 或 NA
    weather_factors_for_prediction.fillna(weather_factors_for_prediction.mean(), inplace=True)

    # 添加预测AQI指数曲线
    # SARIMA模型预测
    train = merged_data[['datetime', '石家庄']].copy()  # 使用历史AQI数据作为时间序列
    train.set_index('datetime', inplace=True)
    order = (1, 1, 1)  # ARIMA模型的阶数
    seasonal_order = (1, 1, 1, 12)  # 季节性阶数
    model_sarima = SARIMAX(train, order=order, seasonal_order=seasonal_order, enforce_stationarity=False, enforce_invertibility=False)
    results_sarima = model_sarima.fit()


    # 定义当前日期为2024年3月1日
    current_date = pd.to_datetime('2024-03-01')

    # 计算未来两三个月后的日期
    future_date = current_date + pd.DateOffset(months=3)  # 或者 pd.DateOffset(months=3) 如果是未来三个月

    # 将日期转换为字符串形式
    end_date_str = future_date.strftime('%Y-%m-%d')

    # 将字符串日期转换为 pandas 的 datetime 格式
    end_date = pd.to_datetime(end_date_str)

    print(len(train))

    forecast = results_sarima.predict(start='2024-03-01', end='2024-06-30', dynamic=False)  # 预测未来60天，每天24个小时
    fig.add_trace(go.Scatter(x=forecast.index, y=forecast, mode='lines', name='石家庄预测AQI数据'))

    # 为每个特征创建曲线并添加到图表中
    for feature in selected_features:
        # 获取特征的预测值
        feature_pred = model.predict(weather_factors_for_prediction)
        fig.add_trace(go.Scatter(x=merged_data['datetime'], y=merged_data[feature], mode='lines', name=feature + '_原始值', xaxis='x1'))
        fig.add_trace(go.Scatter(x=merged_data['datetime'], y=feature_pred, mode='lines', name=feature + '_预测值', xaxis='x1'))

    # 设置图表布局
    fig.update_layout(
        title='石家庄AQI指数预测',
        xaxis_title='日期',
        yaxis_title='AQI指数'
    )

    # 显示图表
    fig.show()


if __name__ == '__main__':
    freeze_support()
    main()