第6章 电力窃漏电用户自动识别

数据挖掘建模目的

1.1 归纳出窃漏电用户的关键特征，构建窃漏电用户的识别模型。

1.2 利用实时监测数据，调用窃漏电用户识别模型实现实时诊断。

分析方法与过程

2.1 数据抽取

从营销系统、自动化设备及往年的窃漏电用户数据抽取数据。（实际上包含了自动抽取和人工标示两种方法）

2.2 数据探索

初步分析的过程，可以说是技术选型阶段。而本案例中采用的是分布分析和周期性分析方法。

2.2.1 分布分析

根据往年的窃漏电用户数据中窃漏电用户的行业分布进行分类，做出柱状图进行分析。（以groupby方法可以实现分类目的）

2.2.2 周期性分析

通过对正常用户和窃漏电用户的用电行为进行分析，总结两者的行为特征，为区分正常用户和窃漏电用户做为判断标准。

其中用电量趋势图可以通过plot方法获得，较为简单。

2.3 数据预处理

2.3.1 数据清洗

把不具备特征性的数据首先进行清除。

2.3.2 缺失值处理

缺失值常见的处理方法分为：删除记录，数据插补，不处理。

插补方法	描述
均值/中位数/众数	根据数据类型进行选择
适用固定值	对于明确的常量可以使用固定值插补
最近临插补	在记录中寻找最接近的样本进行插补，比如一些连续的数据，有规律的数据
回归方法	根据拟合模型来预测缺失量（实质就是做出一个函数来算出缺失值）
插值法	根据已知点建立插值函数，插补缺失值

常用的插值法有：拉格朗日插值法和牛顿插值法，Hermite插值，分段插值，样条插值法等。

拉格朗日插值法和牛顿插值法结果是相同的，所以scipy中仅提供了拉格朗日插值法函数（from scipy.interpolate import lagrange）。

对于不能忽略的缺失值应当进行处理补偿，本案例采用的是拉格朗日插值法对缺失值进行插补。（第4章P61页具体介绍拉格朗日插值法）

2.3 数据变换

2.3.1 电量趋势下降指标

问题：T值正确的计算公式是什么？

遗留问题：通过计算获得ki后，即使正常使用情况下，ki值也会出现上升和下降，那么本指标的作用呢？我认为ki>= ki-1 值应当为-1。或者最终以T占总量的百分比确定是否异常

2.3.2 线损指标

以当日为基点，（前五日的平均线损率-后五日的平均线损率）/后五日平均线损率，进行比较。（对应代码P71页）

2.4 构件专家样本

2.5 模型构建

2.5.1 构件窃漏电用户识别模型

构件分类预测模型（LM神经网络和CART决策树），利用80%做为训练样本，20%做为测试样本。

2.5.2 LM神经网络

LM神经网络模型完整代码如下（keras 2.0.8 backend tf）：

import pandas as pd
from random import shuffle

from keras.models import Sequential
from keras.layers.core import Dense, Activation

datafile = r"E:\Download\百度云\图书配套数据、代码\chapter6\demo\data\model.xls"
data = pd.read_excel(datafile)

netfile = r"E:\\net.h5"  # 模型保存文件

# 读取数据
data = data.as_matrix()  # 转换为矩阵，类型为numpy.ndarray
shuffle(data)  # 打乱矩阵
# 形成训练样本和测试样本
p = 0.8  # 训练样本比例
train = data[:int(len(data) * p), :]
test = data[int(len(data) * p):, :]
# 建立神经网络
net = Sequential()  # 建立神经网络
net.add(Dense(input_dim=3, units=10))
net.add(Activation('relu'))  # 输入层使用relu激活函数
net.add(Dense(input_dim=10, units=1))
net.add(Activation('sigmoid'))  # 隐藏层使用relu激活函数
# 编译神经网络
# net.compile(loss = "binary_crossentropy", optimizer='adam')
net.compile(loss="binary_crossentropy", optimizer='adam', metrics=['accuracy'])  # 常用方式
# 训练模型
hist = net.fit(train[:, :3], train[:, 3], epochs=1000, batch_size=1)  # 条件，结果，训练次数,批尺寸
# 保存模型
net.save_weights(netfile)  # 不能使用中文路径
# 读取模型
from keras.models import load_model

net = load_model("E:\\net.h5")

predict_result = net.predict_classes(train[:, :3]).reshape(len(train))

from ch06.cm_plot import *

cm_plot(train[:, 3], predict_result).show()

# ROC曲线评价
from  sklearn.metrics import roc_curve
from matplotlib import pyplot as plt

predict_result_test = net.predict(test[:, :3]).reshape(len(test))    # 生成预测结果
fpr, tpr, thresholds = roc_curve(test[:, 3], predict_result_test, pos_label=1)    # 实际结果，预测结果，

plt.plot(fpr, tpr, linewidth=2, label='ROC OF LM')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.ylim(0, 1.05)
plt.xlim(0, 1.05)
plt.legend(loc=4)
plt.show()

2.5.3 CART决策树模型

CART的ROC模型评价结果不稳定，原因未知。

CART决策树模型完整代码如下（scikit-learn 0.18.1）：

from sklearn.tree import DecisionTreeClassifier
import pandas as pd
from random import shuffle

datafile = r"E:\Download\百度云\图书配套数据、代码\chapter6\
4000
demo\data\model.xls"
data = pd.read_excel(datafile)

treefile = r"E:\tree.pkl"  # 模型保存文件

# 读取数据
data = data.as_matrix()  # 转换为矩阵，类型为numpy.ndarray
shuffle(data)  # 打乱矩阵
# 形成训练样本和测试样本
p = 0.8  # 训练样本比例
train = data[:int(len(data) * p), :]
test = data[int(len(data) * p):, :]

tree = DecisionTreeClassifier()  # 建立模型
tree.fit(train[:, :3], train[:, 3])  # 训练模型

# 保存模型
from sklearn.externals import joblib

joblib.dump(tree, treefile)

# 读取模型
tree = joblib.load(treefile)
predict_result = tree.predict(train[:, :3])

from ch06.cm_plot import *

cm_plot(train[:, 3], tree.predict(train[:, :3])).show()

# ROC曲线评价
from sklearn.metrics import roc_curve
from matplotlib import pyplot as plt

# predict_result_test = tree.predict(test[:, :3]).reshape(len(test))    # 预测样本的分类（是，不是）
predict_result_test = tree.predict_proba(test[:, :3])[:, 1]    # 预测样本分类的可能性（是的可能性，不是的可能性）
fpr, tpr, thresholds = roc_curve(test[:, 3], predict_result_test, pos_label=1)

plt.plot(fpr, tpr, linewidth=2, label='ROC OF CART')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.ylim(0, 1.05)
plt.xlim(0, 1.05)
plt.legend(loc=4)
plt.show()