数据探索-工业蒸汽量预测-阿里云天池大赛

包含几个大的模块:
1、变量箱型图
2、采用模型预测的形式找出异常样本
3、绘制训练数据集中所有变量的直方图和Q-Q图(查看变量是否符合正态分布)
4、绘制KDE分布图,可以查看并对比训练集和测试集中特征变量的分布情况,发现两个数据集中分布不一致的特征变量
5、计算变量与target之间的相关性系数并用热力图的形式显示
6、根据相关系数筛选特征变量
7、做Box-Cox变换,使变量分布更接近正态分布

 

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

from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error

warnings.filterwarnings("ignore")

pd.set_option('display.max_columns', 1000)

# 读取数据
train_data_file = ".//data//zhengqi_train.txt"
test_data_file = ".//data//zhengqi_test.txt"
train_data = pd.read_csv(train_data_file, sep='\t', encoding='utf-8')
test_data = pd.read_csv(test_data_file, sep='\t', encoding='utf-8')

print("训练集数据信息:")
# 查看训练集数据信息
print(train_data.info())
# 训练集统计信息
print(train_data.describe())
# 查看训练集数据字段信息
print(train_data.head())

print("测试集数据基本信息:")
# 查看测试集数据基本信息
print(test_data.info())
# 测试集统计信息
print(test_data.describe())
# 查看测试集数据字段信息
print(test_data.head())

####################################################################
# 可视化数据分布
print("可视化数据分布:")

# 变量箱型图
# 绘制38个特征变量V0~V37的箱形图
column = train_data.columns.tolist()[:39]  # 列表头
print("column:", column)
fig = plt.figure(figsize=(80, 60), dpi=75)  # 指定绘图对象宽度和高度
# 从图中发现数据存在许多偏离较大的异常,可以考虑移除
for i in range(38):
    plt.subplot(7, 8, i + 1)  # #7行8列的第i个子图
    sns.boxplot(train_data[column[i]], orient="v", width=0.5)  # 箱式图
    plt.ylabel(column[i], fontsize=36)
plt.savefig("./images/变量箱型图.png")
plt.show()

####################################################################
# 采用模型预测的形式找出异常样本

# function to detect outliers based on the predictions of a model
def find_outliers(model, X, y, sigma=3):

    # predict y values using model
    try:
        y_pred = pd.Series(model.predict(X), index=y.index)
    # if predicting fails, try fitting the model first
    except:
        model.fit(X, y)
        y_pred = pd.Series(model.predict(X), index=y.index)
        print("y_pred:\n", type(y_pred), '\n', y_pred)

    # calculate residuals between the model prediction and true y values
    resid = y - y_pred
    mean_resid = resid.mean()
    std_resid = resid.std()

    # calculate z statistic, define outliers to be where |z| > sigma
    z = (resid - mean_resid) / std_resid
    outliers = z[abs(z) > sigma].index

    # print and plot the results
    print('R2=', model.score(X, y))
    print("mse=", mean_squared_error(y, y_pred))
    print('---------------------------------------')

    print('mean of residuals:', mean_resid)
    print('std of residuals:', std_resid)
    print('---------------------------------------')

    print(len(outliers), 'outliers index:')
    print(outliers.tolist())

    plt.figure(figsize=(15, 5))
    ax_131 = plt.subplot(1, 3, 1)
    plt.plot(y, y_pred, '.')
    plt.plot(y.loc[outliers], y_pred.loc[outliers], 'ro')
    plt.legend(['Accepted', 'Outlier'])
    plt.xlabel('y')
    plt.ylabel('y_pred')

    ax_132 = plt.subplot(1, 3, 2)
    plt.plot(y, y - y_pred, '.')
    plt.plot(y.loc[outliers], y.loc[outliers] - y_pred.loc[outliers], 'ro')
    plt.legend(['Accepted', 'Outlier'])
    plt.xlabel('y')
    plt.ylabel('y - y_pred')

    ax_133 = plt.subplot(1, 3, 3)
    # bins设置长条形的数目
    z.plot.hist(bins=50, ax=ax_133)
    z.loc[outliers].plot.hist(color='r', bins=50, ax=ax_133)
    plt.legend(['Accepted', 'Outlier'])
    plt.xlabel('z')

    plt.savefig('./images/outliers.png')
    plt.show()

    return outliers

# 通过岭回归模型找出异常值,并绘制其分布
# 也可以采用其它回归模型代替岭回归模型
X_train = train_data.iloc[:, 0:-1]
y_train = train_data.iloc[:, -1]
outliers = find_outliers(Ridge(), X_train, y_train)

###############################################################################
# Q-Q图是指数据的分位数和正态分布的分位数对比参照的图,如果数据符合正态分布,则所有的点都会落在直线上。
# Q-Q图是通过比较数据和正态分布的分位数是否相等来判断数据是不是符合正态分布

# 绘制训练数据集中所有变量的直方图和Q-Q图
# 设置图的行列数
train_cols = 6
train_rows = len(train_data.columns)
plt.figure(figsize=(4 * train_cols, 3 * train_rows))

i = 0
for col in train_data.columns:
    i += 1
    ax = plt.subplot(train_rows, train_cols, i)
    # fit:设定函数图像,与原图进行比较
    sns.distplot(train_data[col], fit=stats.norm)

    i += 1
    ax = plt.subplot(train_rows, train_cols, i)
    # 红色线条表示正态分布,蓝色线条表示样本数据,蓝色越接近红色参考线,说明越符合预期分布(这是是正态分布)
    res = stats.probplot(train_data[col], plot=plt)  # 默认检测是正态分布

# tight_layout()会调整子图之间的间隔来减少堆叠,
# 解决子图重叠问题
plt.tight_layout()
plt.savefig('./images/所有变量的直方图和Q-Q图.png')
plt.show()

###########################################################################
# KDE(Kernel Density Estimation,核密度估计)可以理解为对直方图的加窗平滑。
# 通过绘制KDE分布图,可以查看并对比训练集和测试集中特征变量的分布情况,发现两个数据集中分布不一致的特征变量。

# 对比所有变量在训练集和测试集中的KDE分布
dist_cols = 6
dist_rows = len(test_data.columns)
plt.figure(figsize=(4 * dist_cols, 4 * dist_rows))

i = 1
for col in test_data.columns:
    ax = plt.subplot(dist_rows, dist_cols, i)
    ax = sns.kdeplot(train_data[col], color='Red', shade=True)
    ax = sns.kdeplot(test_data[col], color='Blue', shade=True)
    ax.set_xlabel(col)
    ax.set_ylabel('Frequency')
    ax = ax.legend(['train', 'test'])
    i += 1

plt.savefig('./images/KDE分布图.png')
plt.show()

################################################################################
# 线性回归关系图主要用于分析变量之间的线性回归关系。

# 查看所有特征变量与target变量的线性回归关系
fcols = 6
frows = len(test_data.columns)
plt.figure(figsize=(5 * fcols, 4 * frows))

i = 0
for col in test_data.columns:
    i += 1
    ax = plt.subplot(frows, fcols, i)
    sns.regplot(x=col, y='target', data=train_data, ax=ax,
                scatter_kws={'marker': '.', 's': 3, 'alpha': 0.3},
                line_kws={'color': 'k'})
    plt.xlabel(col)
    plt.ylabel('target')

    i += 1
    ax = plt.subplot(frows, fcols, i)
    sns.distplot(train_data[col].dropna())
    plt.xlabel(col)

plt.savefig('./images/线性回归关系图.png')
plt.show()

#################################################################################
# 对特征变量的相关性进行分析,可以发现特征变量和目标变量及特征变量之间的关系,为在特征工程中提取特征做准备。
# 但是相关性选择主要用于判别线性相关,对于target变量如果存在更复杂的函数形式的影响,则建议使用树模型的特征重要性去选择。

# 计算变量与target之间的相关性系数
pd.set_option('display.max_columns', 10)
pd.set_option('display.max_rows', 10)
# 特征变量V5,V9,V11,V17,V22,V28在训练集和测试集中的分布不一致,这会导致模型的泛化能力变差,需要删除此类特征。
data_train1 = train_data.drop(['V5', 'V9', 'V11', 'V17', 'V22', 'V28'], axis=1)
train_corr = data_train1.corr()  # 返回df格式
print("相关性系数:\n", train_corr)

# 将相关系数的结果以热力图的形式显示
ax = plt.subplots(figsize=(20, 16))  # 调整画布大小
# vmin, vmax:用于指定图例中最小值与最大值的显示值
# square:bool类型参数,是否使热力图的每个单元格为正方形,默认为False
# annot:指定一个bool类型的值或与data参数形状一样的数组,如果为True,就在热力图的每个单元上显示数值
ax = sns.heatmap(train_corr, vmax=.8, square=True, annot=True)  # 画热力图
plt.savefig('./images/相关系数热力图.png')
plt.show()

########################################################################
# 根据相关系数筛选特征变量
# 寻找K个与target变量最相关的特征变量(K=10)
k = 10  # number of variables for heatmap
cols = train_corr.nlargest(k, 'target')['target'].index

cm = np.corrcoef(train_data[cols].values.T)
hm = plt.subplots(figsize=(10, 10))  #调整画布大小
hm = sns.heatmap(train_data[cols].corr(), annot=True, square=True)
plt.savefig('./images/K个与target变量最相关的特征变量相关系数热力图.png')
plt.show()

# 找出与target变量的相关系数大于0.5的特征变量
threshold = 0.5

corrmat = train_data.corr()
top_corr_features = corrmat.index[abs(corrmat["target"]) > threshold]
plt.figure(figsize=(10, 10))
# cmap:matplotlib的colormap名称或颜色对象;如果没有提供,
# 默认为cubehelix map (数据集为连续数据集时) 或 RdBu_r (数据集为离散数据集时)
g = sns.heatmap(train_data[top_corr_features].corr(), annot=True, cmap="RdYlGn")
plt.savefig('./images/与target变量相关系数大于0.5的特征变量相关系数热力图.png')
plt.show()

# 用相关系数阈值移除相关特征
threshold = 0.5

# 相关系矩阵
corr_matrix = data_train1.corr().abs()
drop_col = corr_matrix[corr_matrix["target"] < threshold].index
#data_all.drop(drop_col, axis=1, inplace=True)

##################################################################################
# 由于线性回归是基于正态分布的,因此在进行统计分析时,需要将数据转换使其符合正态分布。

# Box-Cox变换是统计建模中常用的一种数据转换方法。在连续的响应变量不满足正态分布时,可以使用Box-Cox变换,
# 这一变换可以使线性回归模型在满足线性、正态性、独立性以及方差齐性的同时,又不丢失信息。
# 在对数据做Box-Cox变换之后,可以在一定程度上减少不可观测的误差和预测变量的相关性,
# 这有利于线性模型的拟合及分析出特征的相关性。

# 在做Box-Cox变换之前,需要对数据做归一化预处理。

# 1、合并训练和测试数据集(在归一化时,对数据进行合并操作可以使训练数据和测试数据一致)
drop_columns = ['V5', 'V9', 'V11', 'V17', 'V22', 'V28']
train_x = train_data.drop(['target'], axis=1)

# data_all=pd.concat([train_data,test_data],axis=0,ignore_index=True)
data_all = pd.concat([train_x, test_data])
data_all.drop(drop_columns, axis=1, inplace=True)
print(data_all.head())

# 对合并后的每列数据进行归一化
cols_numeric = list(data_all.columns)

# 离差标准化,将数据调整到[0,1]
def scale_minmax(col):
    return (col-col.min())/(col.max()-col.min())

data_all[cols_numeric] = data_all[cols_numeric].apply(scale_minmax, axis=0)
print(data_all[cols_numeric].describe())

# 2、也可以分开对训练数据和测试数据进行归一化处理,
# 不过这种方式需要建立在训练数据和测试数据分布一致的前提下,建议在数据量大的情况下使用(数据量大,一般分布比较一致),
# 能加快归一化的速度。而数据量较小会存在分布差异较大的情况,此时,在数据分析和线下建模中应该将数据统一归一化。
train_data_process = train_data[cols_numeric]
train_data_process = train_data_process[cols_numeric].apply(scale_minmax, axis=0)

test_data_process = test_data[cols_numeric]
test_data_process = test_data_process[cols_numeric].apply(scale_minmax, axis=0)

# 对特征变量做Box-Cox变换
# 然后,计算分位数并画图展示(基于正态分布),显示特征变量与target变量的线性关系

# 意义解释:Box-Cox变换
# 做回归分析时,通常假设回归方程的残差具有有齐性,即等方差。如果残差不满足齐性,即出现异方差(残差发散),此时就可以通过做Box-Cox变换实现回归方程残差齐性。
# Box-Cox变换实现两项工作:做变换、确定lambda的值(也就是boxcox()函数返回的前两个参数)
# boxcox是将数据分布正态化,使其更加符合后续对数据分布的假设。
# boxcox可以降低skewness值,达到接近正态分布的目标。


cols_numeric_left = cols_numeric[0:13]
cols_numeric_right = cols_numeric[13:]
train_data_process = pd.concat([train_data_process, train_data['target']], axis=1)

fcols = 6
frows = len(cols_numeric_left)
plt.figure(figsize=(4 * fcols, 4 * frows))
i = 0
# 每行展示一个变量,6个图
for var in cols_numeric_left:
    dat = train_data_process[[var, 'target']].dropna()

    # 第1个子图
    i += 1
    plt.subplot(frows, fcols, i)
    sns.distplot(dat[var], fit=stats.norm)
    plt.title(var + ' Original')
    plt.xlabel('')

    # 第2个子图
    i += 1
    plt.subplot(frows, fcols, i)
    _ = stats.probplot(dat[var], plot=plt)
    # stats.skew()计算数据集的样本偏度。
    # 对于正态分布的数据,偏度应该大约为零。
    plt.title('skew=' + '{:.4f}'.format(stats.skew(dat[var])))
    plt.xlabel('')
    plt.ylabel('')

    # 第3个子图
    i += 1
    plt.subplot(frows, fcols, i)
    plt.plot(dat[var], dat['target'], '.', alpha=0.5)
    plt.title('corr=' +
              '{:.2f}'.format(np.corrcoef(dat[var], dat['target'])[0][1]))

    # 第4个子图
    i += 1
    plt.subplot(frows, fcols, i)
    # scipy.stats.boxcox(x, lmbda=None, alpha=None)
    # 参数:
    # lmbda:
    #   如果lmbda不是None,对这个值进行转换
    #   如果lmbda是None,找到最大化对数似然函数的lambda并将其返回为第二个输出参数
    # alpha:
    #   如果alpha不是None,返回lmbda的置信区间作为第三个输出参数。必须是0.0~1.0
    # 返回值:
    # boxcox:nndarray Box-Cox次方转换数组
    # maxlog:float, optional 如果lmbda参数是None,返回的第二个参数是最大化对数似然函数的lambda
    # (min_ci, max_ci):tuple of float, optional 如果lmbda参数是None并且alpha不是None,这个返回的浮点数元组表示在给定alpha下最小和最大置信限制

    #  boxcox要求输入的数据是正值,这里的输入值一般是经过预处理后的,有一个固定的范围,
    #  如果存在非正值,则需要加上一个常数,保证输入值为正值。
    trans_var, lambda_var = stats.boxcox(dat[var].dropna() + 1)
    # print("trans_var, lambda_var", trans_var, lambda_var)
    trans_var = scale_minmax(trans_var)
    # print(trans_var)
    sns.distplot(trans_var, fit=stats.norm)
    plt.title(var + ' Tramsformed')
    plt.xlabel('')

    # 第5个子图
    i += 1
    plt.subplot(frows, fcols, i)
    _ = stats.probplot(trans_var, plot=plt)
    plt.title('skew=' + '{:.4f}'.format(stats.skew(trans_var)))
    plt.xlabel('')
    plt.ylabel('')

    # 第6个子图
    i += 1
    plt.subplot(frows, fcols, i)
    plt.plot(trans_var, dat['target'], '.', alpha=0.5)
    plt.title('corr=' + '{:.2f}'.format(np.corrcoef(trans_var, dat['target'])[0][1]))

# 经过变换后,变量分布更接近正态分布,从图中可以更加直观地看出特征变量与target变量的线性相关性(标题中)
plt.savefig('./images/Box-Cox变换前后直方图和Q-Q图和散点图.png')
plt.show()

数据获取:
https://download.csdn.net/download/fgg1234567890/87832053

阿里云天池大赛工业蒸汽预测是一个时间序列预测问题,可以使用R语言中的时间序列分析和建模工具进行解决。以下是一个简单的R语言代码示例,用于预测未来的蒸汽。 首先,我们需要读入数据并对其进行预处理。这个数据集包含了两个变:日期和蒸汽。 ```r # 读入数据 data <- read.csv("data.csv") # 转换日期格式 data$DATE <- as.Date(data$DATE, format = "%Y/%m/%d") # 将日期设置为数据框的行名 rownames(data) <- data$DATE # 移除日期变 data$DATE <- NULL ``` 接下来,我们可以绘制数据的时间序列图,以便更好地了解数据的性质。 ```r # 绘制时间序列图 plot(data$V1, type = "l", xlab = "日期", ylab = "蒸汽") ``` 然后,我们可以使用时间序列分解方法,将时间序列分解为趋势、季节性和随机性三个部分,并对其进行可视化。 ```r # 时间序列分解 ts.decomp <- decompose(data$V1) # 可视化分解结果 plot(ts.decomp) ``` 分解结果表明,该时间序列具有明显的季节性和趋势,但是随机性较小。 接下来,我们可以使用ARIMA模型进行时间序列预测。ARIMA模型是一种常用的时间序列建模方法,可以用于预测未来的蒸汽。 ```r # 拟合ARIMA模型 arima.model <- arima(data$V1, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 7)) # 预测未来7天的蒸汽 forecast <- predict(arima.model, n.ahead = 7) # 输出预测结果 print(forecast$pred) ``` 以上代码中,我们使用ARIMA(1,1,1)模型,并将季节性设置为7,以便对一周内的季节性进行建模。最后,我们使用predict函数预测未来7天的蒸汽,并输出预测结果。 这是一个简单的R语言示例,用于预测未来的蒸汽。您可以根据实际情况进行修改和扩展,以获得更好的预测结果。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值