\documentclass{ctexart}%中文
\pagestyle{empty}%去掉页眉
\usepackage{graphicx}%插图
\begin{document}
%封面
\begin{flushright}
	{\Large 分\hspace{2cm}数}：\underline{\quad\qquad\qquad}
	\vskip 0.5cm
	{\Large 任课教师签字}：\underline{\quad\qquad\qquad}
\end{flushright} 
\begin{center}
	\quad \\
	\quad \\
	\heiti \fontsize{30}{17} 华\quad 北\quad 电\quad 力\quad 大\quad 学\quad 研\quad 究\quad 生\quad 结\quad 课\quad 作\quad 业
	\vskip 6cm
	%\heiti \zihao{2} 在此打印论文题目，二号黑体	
\end{center}
\vskip 3cm
\begin{quotation}
	\songti \fontsize{15}{15}
	\par\setlength\parindent{8em}
	\quad 
	
	学\hspace{0.2cm} 年\hspace{0.2cm}学\hspace{0.2cm}期：\underline{{\Large 2019-2020}学年第二学期}
	\vskip 0.5cm
	课\hspace{0.2cm} 程\hspace{0.2cm}名\hspace{0.2cm}称：\underline{数据仓库与数据挖掘\qquad}
	\vskip 0.5cm
	学\hspace{0.2cm} 生\hspace{0.1cm} 姓\hspace{0.1cm} 名：\underline{\qquad\qquad\qquad 肖雄\qquad\qquad\qquad }
	\vskip 0.5cm
	学\hspace{1.7cm} 号：\underline{\quad\qquad {\Large 2192221067}\qquad\qquad\quad}
	\vskip 0.5cm
	提\hspace{0.3cm}交\hspace{0.3cm}时\hspace{0.2cm}间：\underline{\qquad {\Large 2020}年{\Large 06}月{\Large 15}日\qquad}
	\vskip 2cm
	\centering
\end{quotation}
\title{逻辑回归分类预测的分析与应用}
\date{}
\maketitle

%中文摘要
\begin{abstract}
	基于逻辑回归模型,对二分类问题进行分类预测;通过sklearn逻辑回归库与梯度下降算法实现做对比;并分别采用留一法与十折交叉验证法对Iris数据集和Blood Transfusion Service Center数据集进行分类;实验结果表明，留一法与十折交叉验证法精度相差不大，但十折交叉验证更加高效。
	\newline%另起一行
	
	\centering%使得关键字居中
	\textbf{关键字：}逻辑回归，梯度下降，留一法，十折交叉验证法
\end{abstract}
%英文摘要
\newcommand{\enabstractname}{Abstract}
\newenvironment{enabstract}{%
	\par\small
	\noindent\mbox{}\hfill{\bfseries \enabstractname}\hfill\mbox{}\par
	\vskip 2.5ex}{\par\vskip 2.5ex}  
\begin{enabstract}
	Based on logistic regression model, the classification prediction of dichotomy problem is carried out.The sklearn logistic regression library is compared with the implementation of gradient descent algorithm.And one method and one thousand one hundred percent cross validation method were used respectively to Iris data set and Blood Transfusion Service Center data set classification;The experimental results show that the accuracy of the retention method is not different from that of the ten fold cross validation method, but the ten fold cross validation method is more efficient.
	
	\centering
	\textbf{Keywords:} Logistic regression, gradient descent, retention method, ten fold cross validation
\end{enabstract}

\section{引言}
	随着信息化社会的高速发展，信息多元化成为主要发展模式。人们使用更多的特征属性描述数据信息，通常某一数据记录使用成千上万的特征描述。在数据挖掘领域中，有众多算法模型对数据进行特征提取分类,针对不同的评估方法，对模型的要求和测试的结果也不一样，逻辑回归模型作为一种高效、易实现的模型应用十分广泛。本文基于逻辑回归模型分别用sklearn逻辑回归库和梯度下降实现二分类问题；对留一法和十折交叉验证法分别评估做比较。实验结果表明，十折交叉验证法在满足精度的同一条件下，耗时更小。
\section{逻辑回归模型}
	逻辑回归是比较常用的机器学习方法，用于估计某种事物的可能性。比如某用户购买某商品的可能性，某病人患有某种疾病的可能性，以及某广告被用户点击的可能性等\cite{1}。逻辑回归延伸了多元线性回归思想，即因变量是二值的情形，自变量为$x_{1}$,$x_{2}$，$x_{3}$,…,$x_{k}$。逻辑回归是用来测量分类结果与因变量之间的关系。逻辑回归模型的最终结果为0,1分类结果。其中1表示属于该类，0表示不属于该类别。
	\par{一般线性模型：}$$f(x)=\omega_{1}x_{1}+\omega_{2}x_{2}…+\omega_{d}x_{d}+b$$
	\par{其中$x_{1}...x_{d}$表示d个特征属性值，$\omega_{1}...\omega_{d},b$表示特征属性参数值。用向量简化为：$$f(x)={\omega}x^T+b$$}
	\par{也可表示为：$$f(x)={\beta}\hat{x}^T$$}
	\par{其中$\beta=({\omega};b)$,$\hat{x}=(x;1)$,线性回归只能预测连续的值，对于离散的二分类问题需要转化为逻辑回归，即将线性结果映射到\{0,1\}上:$$\ln\frac{f(x)}{1-f(x)}={\beta}\hat{x}^T$$}
	\par{因此，逻辑回归又称为对数几率回归：$$f(x)=\frac{1}{1+e^{-\beta\hat{x}^T}}$$}
	\par{对于逻辑回归模型，关键在于如何求得$\beta$的值，最常用的方法就是梯度下降法。}
	%文献索引
	\cite{2}
\section{梯度下降}
	梯度下降法是求解无约束优化问题的方法之一，有计算过程简单、初始收敛较快等优点，因此也常作为其他算法的核心算法，例如人工神经网络和逻辑回归，广泛应用于数据挖掘、模式识别等领域\cite{3}。
	\par{对于一阶无约束优化问题$min_{x}f(x)$,若能找到$x^{0},x^{1},x^{2}...$满足：$$f(x^{t+1})<f(x^{t}),t=0,1,2...$$}
	\par{不断执行此过程可收敛到局部极小点，根据泰勒展开式：$$f(x+\Delta x)\approx f(x)+\Delta x\nabla f(x)$$}
	\par{于是，欲满足$f(x+\Delta x)<f(x)$，可选择$$\Delta x=-\gamma\nabla f(x)$$}
	\par{其中$\gamma$是小常数。这就是梯度下降}\cite{4}。
	\par{在逻辑回归模型中，可通过极大似然法来估计$\beta$的值\cite{5}：$$\psi(\beta)=\sum_{i=1}^{m}(y_{i}p(y=1|\beta\hat{x}^T_{i})+(1-y_{i})p(y=0|\beta\hat{x}^T_{i}))=\sum_{i=1}^{m}(y_{i}\beta\hat{x}^T_{i}-ln(1+e^{\beta\hat{x}^T_{i}}))$$}
	\par{该函数为连续可导凸函数，因此可采用梯度下降来求解，因此将上式转化为最小化：$$\psi(\beta)=\sum_{i=1}^{m}(-y_{i}\beta\hat{x}^T_{i}+ln(1+e^{\beta\hat{x}^T_{i}}))$$}
	\par{迭代过程：$$\beta^{t+1}=\beta^{t}-\gamma\nabla \psi(\beta)$$}
\section{应用}
	本文使用文献2中西瓜数据集$3.0\alpha$，分布情况如图\ref{fig:xigua}
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.7\linewidth]{1}
		\caption{西瓜数据集$3.0\alpha$散点图\label{fig:xigua}}
		\label{fig:1}
	\end{figure}
	\par{特征属性为密度和含糖量，样本标签“1”表示好瓜，“0”表示坏瓜；使用留出法选择相同的训练集和测试集，通过sklearn逻辑回归库和批量梯度下降法分别进行分类训练和测试，其中梯度下降设置固定步长为0.1，在迭代15000次后趋于稳定，如图\ref{fig:diedai}}
	\begin{figure}
		\centering
		\includegraphics[width=0.7\linewidth]{2}
		\caption{批量梯度下降迭代曲线\label{fig:diedai}}
		\label{fig:2}
	\end{figure}。
	\par{经测试两者精度相差不大，测试结果如下表:}
	\par{}
	%表格
	\begin{tabular}{|c|c|}
		\hline 
		方法&精度  \\ 
		\hline 
		sklearn逻辑回归库&67\%  \\ 
		\hline 
		批量梯度下降法&66.67\%  \\ 
		\hline 
	\end{tabular} 
	\par{考虑到该数据样本过小，使用留出法拟合效果一般，因此在UCI选择Iris数据集和Blood Transfusion Service Center数据集，同样以逻辑回归模型分别对留一法和十折交叉验证法评估做比较，结果如下：}
	\par{}
	\begin{tabular}{|c|c|c|}
		\hline 
		数据集&留一法&十折交叉验证法  \\ 
		\hline 
		Iris&96.66\%&97.33\%\\ 
		\hline 
		Blood Transfusion Service Center&76.87\%&77.01\%\\ 
		\hline 
	\end{tabular} 
	\par{Iris数据集因为数据类间分散情况比较好，广泛被引用，因此拟合效果比Blood Transfusion Service Center数据集要好，从图\ref{fig:i},\ref{fig:t}就能看出：}
	\begin{figure}
		\centering
		\includegraphics[width=0.7\linewidth]{i}
		\caption{Iris数据集散点图\label{fig:i}}
		\label{fig:i}
	\end{figure}
	\begin{figure}
		\centering
		\includegraphics[width=0.7\linewidth]{t}
		\caption{Blood Transfusion Service Center数据集散点图\label{fig:t}}
		\label{fig:t}
	\end{figure}
	
	
	\par{实验结果表明，留一法和十折交叉验证法的精度相差不大。值得注意的是，十折交叉验证更加高效，耗时更少，对于数据量越大，这种现象越明显。因此选择十折交叉验证即可满足精度要求，又减少运行成本。}
\section{总结}
	逻辑回归对于二分类问题，不仅将预测值映射到$\{0,1\}$之间的值，而且还能评估出概率值，这使得在许多领域都具有广泛的应用。一个好的模型通常是由测试结果来判定，因此对训练集、测试集进行划分的评估方法起着决定性作用，本文对留一法和十折交叉验证法评估结果比较，结果显示在数据量足够的情况，选择十折交叉验证法更加高效。
	
	
\begin{thebibliography}{}
	\bibitem{1}毛林,陆全华,程涛.\emph{基于高维数据的集成逻辑回归分类算法的研究与应用[J]},
	\texttt{科技通报,2013,29(12):64-66}
	\bibitem{2}周志华.\emph{机器学习[M]},
	\texttt{2016:53-59}
	\bibitem{3}郭跃东,宋旭东.\emph{梯度下降法的分析和改进[J]},
	\texttt{科技展望,2016,26(15):115+117}
	\bibitem{4}周志华.\emph{机器学习[M]},
	\texttt{2016:407-408}
	\bibitem{5}周志华.\emph{机器学习[M]},
	\texttt{2016:59-60}
\end{thebibliography}

\end{document}

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

iris = sns.load_dataset('iris') # 在线载入自带的 iris 数据集
X = iris.values[:, 0 : 4]
y = iris.values[:, 4]

sns.set(style='white') # 风格设置
g = sns.pairplot(iris, hue='species', markers=['o', 's', 'D']) # 变量关系组图
plt.show()

from sklearn.linear_model import LogisticRegression
from sklearn import metrics
from sklearn import model_selection

log_model = LogisticRegression(max_iter=1000) # 增加最大迭代次数,也可以减少数据量
m, n = np.shape(X)

y_pred_10_fold = model_selection.cross_val_predict(log_model, X, y, cv=10)

# 打印精度
accuracy_10_fold = metrics.accuracy_score(y, y_pred_10_fold)
print('The accuracy of 10-fold cross-validation:', accuracy_10_fold)

accuracy_LOO = 0
# 计算 m 次测试的结果
for train_index, test_index in model_selection.LeaveOneOut().split(X):
    X_train, X_test = X[train_index], X[test_index] # 训练集样本,测试集样本
    y_train, y_test = y[train_index], y[test_index] # 训练集标签, 测试集标签
    log_model.fit(X_train, y_train) # 训练模型
    y_pred_LOO = log_model.predict(X_test) # 测试
    if y_pred_LOO == y_test:
        accuracy_LOO += 1
print('The accuracy of Leave-One-Out:', accuracy_LOO / m)

import numpy as np
import matplotlib.pyplot as plt


dataset = np.loadtxt('/home/data/watermelon_3a.csv', delimiter=',')

X = dataset[:, 1 : 3]
y = dataset[:, 3]
print(np.shape(X))

f1 = plt.figure(1)
plt.title('watermelon_3a')
plt.xlabel('density')
plt.ylabel('rate_sugar')
plt.scatter(X[y == 0, 0], X[y == 0, 1], marker = 'o', color = 'k', s = 100, label= 'bad')
plt.scatter(X[y == 1, 0], X[y ==1, 1], marker= 'o', color = 'g', s = 100, label = 'good')
plt.legend(loc = 'upper right')
plt.show()

from sklearn import model_selection
from sklearn.linear_model import LogisticRegression
from sklearn import metrics
import matplotlib.pylab as pl

# 切分数据集:留出法 返回 划分好的训练集测试集样本和训练集测试集标签
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.5, random_state=0)

# 训练模型
log_model = LogisticRegression()
log_model.fit(X_train, y_train)

# 模型测试
y_pred = log_model.predict(X_test)

print(metrics.confusion_matrix(y_test, y_pred))
print(metrics.classification_report(y_test, y_pred))
precision, recall, thresholds = metrics.precision_recall_curve(y_test, y_pred)

f2 = plt.figure(2)
h = 0.01
x0_min, x0_max = X[:, 0].min() - 0.1, X[:, 0].max() + 0.1
x1_min, x1_max = X[:, 1].min() - 0.1, X[:, 1].max() +0.1
x0, x1= np.meshgrid(np.arange(x0_min, x0_max, h), np.arange(x1_min, x1_max, h)) # 生成笛卡尔积坐标矩阵

z = log_model.predict(np.c_[x0.ravel(), x1.ravel()]) # c_ 按列合并, ravel 降成一维

z = z.reshape(x0.shape)
plt.contourf(x0, x1, z, cmap = pl.cm.Paired)# 等高线

plt.title('watermelon_3a')
plt.xlabel('density')
plt.ylabel('rate_sugar')
plt.scatter(X[y==0, 0], X[y==0, 1], marker='o', color='k', s=100, label='bad')
plt.scatter(X[y==1, 0], X[y==1, 1], marker='o', color='g', s=100, label='good')
plt.legend(loc='upper right')
plt.show()

# 1)实现 P59 公式3.27极大似然法

def likelihood_sub(x, y, beta):
    """
    :param x: 一个示例变量(行向量)
    :param y:一个样品标签(行向量)
    :param beta:3.27中矢量参数(行向量)
    :return: 单个对数似然 3.27
    """
    return -y * np.dot(beta, x.T) + np.math.log(1 + np.math.exp(np.dot(beta, x.T)))

def likelihood(X, y, beta):
    """
    公式 3.27 :对数似然函数(交叉熵损失函数)
    :param X: 示例变量矩阵
    :param y:样本标签矩阵
    :param beta:3.27中的矢量参数
    :return: beta 的似然值
    """
    sum = 0
    m, n = np.shape(X)

    for i in range(m):
        sum += likelihood_sub(X[i], y[i], beta)
    return sum

# 2)实现似然公式一阶偏导
def sigmoid(x, beta):
    """
    基础模型 S 形函数
    P59 对数几率回归(逻辑回归)公式 3.23
    :param x: 预测变量
    :param beta: beta 变量
    :return:S 形函数
    """
    return 1 / (1 + np.math.exp(- np.dot(beta, x.T)))
def partial_derivative(X, y, beta):
    """
    P60 似然公式一阶偏导3.30
    :param X:示例变量矩阵
    :param y:样本标签矩阵
    :param beta:3.27 中矢量参数
    :return: beta 的偏导数,梯度
    """
    m, n = np.shape(X)
    pd = np.zeros(n)

    for i in range(m):
        tmp = -y[i] + sigmoid(X[i], beta)
        for j in range(n):
            pd[j] += X[i][j] * tmp
    return pd

def gradDscent(X, y, alpha, iterations, n):
    """
    :param X:变量矩阵
    :param y:样本标签数组
    :return:3.27中beta参数最优解
    """
    cost = np.zeros(iterations) # 构建 max_times 个 0 的数组
    beta = np.mat(np.zeros(n)) # 初始化 beta

    for i in range(iterations):
        # 梯度下降
        output = partial_derivative(X, y, beta)
        beta = beta - alpha * output
        cost[i] = likelihood(X, y, beta)

    return beta, cost

def showConvergCurve(Iterations, Cost):
    """
    :param Iterations: 迭代次数
    :param Cost: 损失值数组
    """
    f1 = plt.figure(1)
    t = np.arange(Iterations)
    p1 = plt.subplot(1,1,1)
    p1.plot(t, Cost, 'r')
    p1.set_xlabel('Iterations')
    p1.set_ylabel('cost')
    p1.set_title('The Gradient Descent Convergence Curve')

    plt.show()

def showLogRegression(X, y, Beta, N):
    f2 = plt.figure(2)

    plt.title('The Logistic Regression Fitted Curve')
    plt.xlabel('density')
    plt.ylabel('rate_sugar')
    # f = Beta * X.transpose()
    # plt.plot(X[:, 2], f.tolist()[0], 'r', label = 'Prediction')
    min_x = min(X[:, 0])
    max_x = max(X[:, 0])
    y_min_x = (- Beta.tolist()[0][2] - Beta.tolist()[0][0] * min_x) / Beta.tolist()[0][1] # 由线性模型 y =  w1 * x1 + w2 * x2 +b
    y_max_x = (- Beta.tolist()[0][2] - Beta.tolist()[0][0] * max_x) / Beta.tolist()[0][1]
    plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')
    plt.scatter(X[y == 0, 0], X[y == 0, 1], marker='o', color='k', s=100, label='bad')
    plt.scatter(X[y == 1, 0], X[y == 1, 1], marker='o', color='g', s=100, label='good')
    plt.legend(loc='upper right')
    plt.show()

def testLogRegres(Beta, test_x, test_y):
    m, n = np.shape(test_x)
    matchCount = 0
    for i in range(m):
        predict = sigmoid(test_x[i], Beta) > 0.5
        if predict == bool(test_y[i]):
            matchCount += 1
    accuracy = float(matchCount) / m
    return accuracy

def loadData():
    dataset = np.loadtxt('/home/data/watermelon_3a.csv', delimiter=',')

    X = dataset[:, 1: 3]
    tmp = np.ones(X.shape[0])
    X = np.insert(X, 2, values=tmp, axis=1) # 在最后一列插入全是 1 的列
    y = dataset[:, 3]
    return X, y

def main():
    alpha = 0.1  # 迭代步长
    iterations = 1500  # 迭代次数上限
    X, y = loadData()
    test_x = X
    test_y = y
    m, n = np.shape(X)
    beta, cost = gradDscent(X, y, alpha, iterations, n)
    print(beta)
    showConvergCurve(iterations, cost)
    showLogRegression(X, y, beta, n)
    accuracy = testLogRegres(beta, test_x, test_y)
    print('The classify accuracy is: %.3f%%' %(accuracy * 100))

import numpy as np
import pandas as pd

## 定义结果的加权平均函数
def Weighted_method(test_pre1,test_pre2,test_pre3,w=[1/3,1/3,1/3]):
    Weighted_result = w[0]*pd.Series(test_pre1)+w[1]*pd.Series(test_pre2)+w[2]*pd.Series(test_pre3)
    return Weighted_result

from sklearn import metrics
# 各模型的预测结果计算MAE
print('Pred1 MAE:',metrics.mean_absolute_error(y_test_true, test_pre1))
print('Pred2 MAE:',metrics.mean_absolute_error(y_test_true, test_pre2))
print('Pred3 MAE:',metrics.mean_absolute_error(y_test_true, test_pre3))

## 根据加权计算MAE
w = [0.3,0.4,0.3] # 定义比重权值
Weighted_pre = Weighted_method(test_pre1,test_pre2,test_pre3,w)
print('Weighted_pre MAE:',metrics.mean_absolute_error(y_test_true, Weighted_pre))

## 定义结果的加权平均函数
def Mean_method(test_pre1,test_pre2,test_pre3):
    Mean_result = pd.concat([pd.Series(test_pre1),pd.Series(test_pre2),pd.Series(test_pre3)],axis=1).mean(axis=1)
    return Mean_result

Mean_pre = Mean_method(test_pre1,test_pre2,test_pre3)
print('Mean_pre MAE:',metrics.mean_absolute_error(y_test_true, Mean_pre))

## 定义结果的加权平均函数
def Median_method(test_pre1,test_pre2,test_pre3):
    Median_result = pd.concat([pd.Series(test_pre1),pd.Series(test_pre2),pd.Series(test_pre3)],axis=1).median(axis=1)
    return Median_result

Median_pre = Median_method(test_pre1,test_pre2,test_pre3)
print('Median_pre MAE:',metrics.mean_absolute_error(y_test_true, Median_pre))

from sklearn import linear_model

def Stacking_method(train_reg1,train_reg2,train_reg3,y_train_true,test_pre1,test_pre2,test_pre3,model_L2= linear_model.LinearRegression()):
    model_L2.fit(pd.concat([pd.Series(train_reg1),pd.Series(train_reg2),pd.Series(train_reg3)],axis=1).values,y_train_true)
    Stacking_result = model_L2.predict(pd.concat([pd.Series(test_pre1),pd.Series(test_pre2),pd.Series(test_pre3)],axis=1).values)
    return Stacking_result

## 生成一些简单的样本数据，test_prei 代表第i个模型的预测值
train_reg1 = [3.2, 8.2, 9.1, 5.2]
train_reg2 = [2.9, 8.1, 9.0, 4.9]
train_reg3 = [3.1, 7.9, 9.2, 5.0]
# y_test_true 代表第模型的真实值
y_train_true = [3, 8, 9, 5] 

test_pre1 = [1.2, 3.2, 2.1, 6.2]
test_pre2 = [0.9, 3.1, 2.0, 5.9]
test_pre3 = [1.1, 2.9, 2.2, 6.0]

# y_test_true 代表第模型的真实值
y_test_true = [1, 3, 2, 6] 

model_L2= linear_model.LinearRegression()
Stacking_pre = Stacking_method(train_reg1,train_reg2,train_reg3,y_train_true,
                               test_pre1,test_pre2,test_pre3,model_L2)
print('Stacking_pre MAE:',metrics.mean_absolute_error(y_test_true, Stacking_pre))

# reduce_mem_usage 函数通过调整数据类型，帮助我们减少数据在内存中占用的空间
def reduce_mem_usage(df):
    """ iterate through all the columns of a dataframe and modify the data type
        to reduce memory usage.
    """
    start_mem = df.memory_usage().sum()
    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))

    for col in df.columns:
        col_type = df[col].dtype

        if col_type != object:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
        else:
            df[col] = df[col].astype('category')

    end_mem = df.memory_usage().sum()
    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))
    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))
    return df

sample_feature = reduce_mem_usage(pd.read_csv("data_for_tree.csv"))

# 返回 x in sample_feature.columns not include ['price','brand','model','brand'] 的列表
continuous_feature_names = [x for x in sample_feature.columns if x not in ['price','brand','model','brand']]

sample_feature = sample_feature.dropna().replace('-', 0).reset_index(drop=True)
sample_feature['notRepairedDamage'] = sample_feature['notRepairedDamage'].astype(np.float32)
train = sample_feature[continuous_feature_names + ['price']]

train_X = train[continuous_feature_names]
train_y = train['price']
`

from sklearn.linear_model import LinearRegression

model = LinearRegression(normalize=True)

model = model.fit(train_X, train_y)

# 查看训练的线性回归模型的截距（intercept）与权重(coef)
print('intercept:'+ str(model.intercept_))

print(sorted(dict(zip(continuous_feature_names, model.coef_)).items(), key=lambda x:x[1], reverse=True))

from matplotlib import pyplot as plt
subsample_index = np.random.randint(low=0, high=len(train_y), size=50)
#绘制特征v_9的值与标签的散点图，图片发现模型的预测结果（蓝色点）与真实标签（黑色点）的分布差异较大，
# 且部分预测值出现了小于0的情况，说明我们的模型存在一些问题
plt.scatter(train_X['v_9'][subsample_index], train_y[subsample_index], color='black')
plt.scatter(train_X['v_9'][subsample_index], model.predict(train_X.loc[subsample_index]), color='blue')
plt.xlabel('v_9')
plt.ylabel('price')
plt.legend(['True Price','Predicted Price'],loc='upper right')
print('The predicted price is obvious different from true price')
plt.show()

import seaborn as sns
print('It is clear to see the price shows a typical exponential distribution')
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
sns.distplot(train_y)
plt.subplot(1,2,2)
sns.distplot(train_y[train_y < np.quantile(train_y, 0.9)])
plt.show()

# 在这里我们对标签进行了 log(x+1) 变换，使标签贴近于正态分布
train_y_ln = np.log(train_y + 1)
print('The transformed price seems like normal distribution')
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
sns.distplot(train_y_ln)
plt.subplot(1,2,2)
sns.distplot(train_y_ln[train_y_ln < np.quantile(train_y_ln, 0.9)])
plt.show()

model = model.fit(train_X, train_y_ln)

print('intercept:'+ str(model.intercept_))
sorted(dict(zip(continuous_feature_names, model.coef_)).items(), key=lambda x:x[1], reverse=True)

#再次进行可视化，发现预测结果与真实值较为接近，且未出现异常状况
plt.scatter(train_X['v_9'][subsample_index], train_y[subsample_index], color='black')
plt.scatter(train_X['v_9'][subsample_index], np.exp(model.predict(train_X.loc[subsample_index])), color='blue')
plt.xlabel('v_9')
plt.ylabel('price')
plt.legend(['True Price','Predicted Price'],loc='upper right')
print('The predicted price seems normal after np.log transforming')
plt.show()

##使用线性回归模型，对未处理标签的特征数据进行五折交叉验证
from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_absolute_error,  make_scorer
def log_transfer(func):
    def wrapper(y, yhat):
        result = func(np.log(y), np.nan_to_num(np.log(yhat)))
        return result
    return wrapper

scores = cross_val_score(model, X=train_X, y=train_y, verbose=1, cv = 5, scoring=make_scorer(log_transfer(mean_absolute_error)))

print('AVG:', np.mean(scores))

#使用线性回归模型，对处理过标签的特征数据进行五折交叉验证（
scores = cross_val_score(model, X=train_X, y=train_y_ln, verbose=1, cv = 5, scoring=make_scorer(mean_absolute_error))

print('AVG:', np.mean(scores))

scores = pd.DataFrame(scores.reshape(1,-1))
scores.columns = ['cv' + str(x) for x in range(1, 6)]
scores.index = ['MAE']
print(scores)

# 采用时间顺序对数据集进行分隔 选用靠前时间的4/5样本当作训练集，靠后时间的1/5当作验证集
import datetime
sample_feature = sample_feature.reset_index(drop=True)
split_point = len(sample_feature) // 5 * 4 # 取整除 - 返回商的整数部分（向下取整）

train = sample_feature.loc[:split_point].dropna()
val = sample_feature.loc[split_point:].dropna()

train_X = train[continuous_feature_names]
train_y_ln = np.log(train['price'] + 1)
val_X = val[continuous_feature_names]
val_y_ln = np.log(val['price'] + 1)

model = model.fit(train_X, train_y_ln)
print(mean_absolute_error(val_y_ln, model.predict(val_X)))

# 绘制学习率曲线与验证曲线
from sklearn.model_selection import learning_curve, validation_curve
def plot_learning_curve(estimator, title, X, y, ylim=None, cv=None,n_jobs=1, train_size=np.linspace(.1, 1.0, 5 )):
    plt.figure()
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel('Training example')
    plt.ylabel('score')
    train_sizes, train_scores, test_scores = learning_curve(estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_size, scoring = make_scorer(mean_absolute_error))
    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)
    plt.grid()#区域

    # x：第一个参数表示覆盖的区域，我直接复制为x，表示整个x都覆盖
    # 0：表示覆盖的下限
    # y：表示覆盖的上限是y这个曲线
    # facecolor：覆盖区域的颜色
    # alpha：覆盖区域的透明度[0,1],其值越大，表示越不透明
    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std, alpha=0.1,
                     color="g")
    plt.plot(train_sizes, train_scores_mean, 'o-', color='r',
             label="Training score")
    plt.plot(train_sizes, test_scores_mean,'o-',color="g",
             label="Cross-validation score")
    plt.legend(loc="best")
    return plt

plot_learning_curve(LinearRegression(), 'Liner_model', train_X[:1000], train_y_ln[:1000], ylim=(0.0, 0.5), cv=5, n_jobs=1)
plt.show()

train = sample_feature[continuous_feature_names + ['price']].dropna()

train_X = train[continuous_feature_names]
train_y = train['price']
train_y_ln = np.log(train_y + 1)

# 线性模型 & 嵌入式特征选择
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso

models = [LinearRegression(),
          Ridge(),
          Lasso()]

result = dict()
for model in models:
    model_name = str(model).split('(')[0]
    scores = cross_val_score(model, X=train_X, y=train_y_ln, verbose=0, cv = 5, scoring=make_scorer(mean_absolute_error))
    result[model_name] = scores
    print(model_name + ' is finished')

# 对三种方法的效果对比
result = pd.DataFrame(result)
result.index = ['cv' + str(x) for x in range(1, 6)]
print(result)

model = LinearRegression().fit(train_X, train_y_ln)
print('intercept:'+ str(model.intercept_))
sns.barplot(abs(model.coef_), continuous_feature_names)
plt.show()

model = Ridge().fit(train_X, train_y_ln)
print('intercept:'+ str(model.intercept_))
sns.barplot(abs(model.coef_), continuous_feature_names)
plt.show()

model = Lasso().fit(train_X, train_y_ln)
print('intercept:'+ str(model.intercept_))
sns.barplot(abs(model.coef_), continuous_feature_names)
plt.show()

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from operator import itemgetter

Train_data = pd.read_csv("./datalab/used_car_train_20200313.csv", sep=" ")
Test_data = pd.read_csv("./datalab/used_car_testA_20200313.csv", sep=" ")

print(Train_data.shape)

print(Train_data.head())

print(Train_data.columns)

## 删除异常值
# 这里我包装了一个异常值处理的代码，可以随便调用
def outliers_proc(data, col_name, scale=3):
    """
    用于清洗异常值，默认用box_plot(scale=3)进行清洗
    :param data:接受 pandas 数据格式
    :param col_name:pandas 列名
    :param scale:尺度
    """
    def box_plot_outliers(data_ser, box_scale):
        """
        利用箱线图去除异常值
        :param data_ser:接收 pandas.Series 数据格式
        :param box_scale: 箱线图尺度 （规定大于上四分位数1.5倍四分位数差 的值，或者小于下四分位数1.5倍四分位数差的值，划为异常值）
        """
        iqr = box_scale * (data_ser.quantile(0.75) - data_ser.quantile(0.25)) # 3倍四分位数差
        val_low = data_ser.quantile(0.25) - iqr # 下限=Q1-3IQR
        val_up = data_ser.quantile(0.75) + iqr # 上限=Q3+3IQR
        rule_low = (data_ser < val_low)
        rule_up = (data_ser > val_up) # 返回 pandas.Series 中对应值的bool
        return (rule_low, rule_up), (val_low, val_up)

    data_n = data.copy() # copy 数据
    data_series = data_n[col_name] # 返回指定 col_name 数据
    rule, value = box_plot_outliers(data_series, box_scale=scale)
    index = np.arange(data_series.shape[0])[rule[0]|rule[1]] # 返回rule_low, rule_up中为True的下标的列表
    print("Delete number is:{}".format(len(index))) # 打印下标列表中个数
    data_n = data_n.drop(index) # 删除(删除后下标没变)
    data_n.reset_index(drop=True, inplace=True) # 重置索引（drop=True删除原来的索引;inplace=True当前修改状态应用到原来Series中）
    print("Now column number is:{}".format(data_n.shape[0])) # 查看删除后的数据个数
    index_low = np.arange(data_series.shape[0])[rule[0]]
    outliers = data_series.iloc[index_low] # ilco-按下标进行索引
    print("Description of data larger than the lower bound is:")
    print(pd.Series(outliers).describe())
    index_up = np.arange(data_series.shape[0])[rule[1]]
    outliers = data_series.iloc[index_up]
    print("Description of data larger than the upper bound is:")
    print(pd.Series(outliers).describe())

    fig, ax = plt.subplots(1, 2, figsize=(10, 7)) # 创建子图:1行2列
    sns.boxplot(y=data[col_name], data=data, palette="Set1", ax=ax[0]) # 箱线图
    sns.boxplot(y=data_n[col_name], data=data_n, palette="Set1", ax=ax[1])
    plt.show()
    return data_n


## 我们可以删掉一些异常数据,以 power 为例
## 这里删不删可以自行判断
## 但是注意 Test 的数据不能删
Train_data = outliers_proc(Train_data, "power", scale=3)

# 训练集和测试集放在一起，方便构造特征
Train_data["train"] = 1 # 添加新字段，并设置值为1
Test_data["train"] = 1
data = pd.concat([Train_data,Test_data],ignore_index=True) # 连接函数 ignore_index=True重置索引

# 使用时间：data["createDate"] - data["regDate"], 反应汽车使用时间，一般来说价格与使用时间成反比
# 不过要注意, 数据里有时间出错的格式, 所以我们需要 errors = "coerce"
data["used_time"] = (pd.to_datetime(data["creatDate"], format="%Y%m%d", errors="coerce") -
pd.to_datetime(data["regDate"],format="%Y%m%d",errors="coerce")).dt.days # to_datetime将参数转换为日期 dt.days每个元素的天数

# 看一下空数据, 有 15k 个样本的时间有问题的, 我们可以选择删除, 也可以选择放着
# 但是这里不建议删除, 因为删除缺失数据占总样本量过大, 7.5%
# 我们可以先放着, 因为如果我们 XGBoost 之类的决策树, 其本身就能处理缺失值, 所以可以不用管
print(data["used_time"].isnull().sum())

# 从邮编中提取城市信息, 相当于加入了先验知识
#print(data["regionCode"])
# 增加city 字段, 并从 regionCode 值的倒数第三位切片(apply 对regionCode每个元素运行指定运算 lambda 匿名函数)
data["city"] = data["regionCode"].apply(lambda x : str(x)[:-3])

# 计算某品牌的销售统计量, 还可以计算其他特征的统计量
# 这里以 train 的数据计算统计量
Train_gb = Train_data.groupby("brand") # 分组
all_info = {}
for kind, kind_data in Train_gb:
    info = {}
    kind_data = kind_data[kind_data["price"] > 0]
    info["brand_amount"] = len(kind_data)
    info["brand_price_max"] = kind_data.price.max()
    info["brand_price_median"] = kind_data.price.median()
    info["brand_price_min"] = kind_data.price.min()
    info["brand_price_sum"] = kind_data.price.sum()
    info["brand_price_std"] = kind_data.price.std() # 样本方差
    info["brand_price_average"] = round(kind_data.price.sum() / (len(kind_data)+1), 2) # round(2)取近似值保留两位数
    all_info[kind] = info

brand_fe = pd.DataFrame(all_info).T.reset_index().rename(columns={"index":"brand"}) # T转置   reset_index还原索引  rename并修改列名
data = data.merge(brand_fe, how="left", on="brand") # 合并数据 data链接在brand_fe "brand"字段左边

# 数据分桶 以 power 为例
# 这时候缺失值也进桶了
# 为什么要分桶：
# 1. 离散后稀疏向量内积乘法运算速度更快, 计算结果也方便存储, 容易扩展
# 2. 离散后的特征对异常值更具鲁棒性, 如 age>30 为 1 否则为 0 , 对于年龄为 200 的也不会对模型造成很大的干扰
# 3. LR 属于广义线性模型, 表达能力有限, 经过离散化后, 每个变量有单独的权重, 这相当于引入了非线性, 能够提升模型的表达能力, 加大拟合
# 4. 离散后特征可以进行特征交叉, 提升表达能力, 由 M+N 个变量变成　Ｍ*N 个变量, 进一步引入非线性, 提升了表达能力
# 5. 特征离散后模型更稳定, 如用户年龄区间, 不会因为用户年龄长了一岁就变化
# 当然还有很多原因,　LightGBM 在改进 XGBoost 时就增加了数据分桶, 增强了模型的泛化性

bin = [i*10 for i in range(31)]
data["power_bin"] = pd.cut(data["power"], bin, labels=False) # 分桶 cut切分数据(必须是一维的) bin定义区间 labels=False返回第几个bin（从0开始）
print(data[["power_bin", "power"]].head())

# 删除不需要的数据
data = data.drop(["creatDate", "regDate", "regionCode"], axis=1) # drop函数默认删除行，列需要加axis = 1
print(data.shape)

print(data.columns)

# 目前的数据其实已经可以给树模型使用了, 所以我们导出一下
data.to_csv("data_for_tree.csv", index=0) # index=0不保存行索引

# 我们可以再构造一份特征给 LR NN 之类的模型用
# 之所以分开构造是因为, 不同模型对数据的要求不同
# 先看下数据分布：
data["power"].plot.hist()
plt.show()

# 我们刚刚已经对 train 进行异常值处理了，但是现在还有这么奇怪的分布是因为 test 中的 power 异常值，
# 所以我们其实刚刚 train 中的 power 异常值不删为好，可以用长尾分布截断来代替
Train_data["power"].plot.hist()
plt.show()

# 我们对其取 log, 再做归一化
from sklearn import preprocessing
# 将数据的每一个特征缩放到给定的范围，将数据的每一个属性值减去其最小值，然后除以其极差（最大值 - 最小值）
min_max_scaler = preprocessing.MinMaxScaler()
data["power"] = np.log(data["power"] + 1)
# 归一化：(0,1)标准化
data["power"] = ((data["power"] - np.min(data["power"])) / (np.max(data["power"]) - np.min(data["power"])))
data["power"].plot.hist()
plt.show()

# km 的比较正常, 应该已经做过分桶了
data["kilometer"].plot.hist()
plt.show()

# 所以可以直接作归一化
data["kilometer"] = ((data["kilometer"] - np.min(data["kilometer"])) /
                     (np.max(data["kilometer"]) - np.min(data["kilometer"])))
data["kilometer"].plot.hist()
plt.show()

# 除此之外 还有我们刚刚构造的统计量特征：
# 'brand_amount', 'brand_price_average', 'brand_price_max',
# 'brand_price_median', 'brand_price_min', 'brand_price_std',
# 'brand_price_sum'
# 这里不再一一举例分析了，直接做变换，
def max_min(x):
    return (x - np.min(x)) / (np.max(x) - np.min(x))

data['brand_amount'] = ((data['brand_amount'] - np.min(data['brand_amount'])) /
                        (np.max(data['brand_amount']) - np.min(data['brand_amount'])))
data['brand_price_average'] = ((data['brand_price_average'] - np.min(data['brand_price_average'])) /
                               (np.max(data['brand_price_average']) - np.min(data['brand_price_average'])))
data['brand_price_max'] = ((data['brand_price_max'] - np.min(data['brand_price_max'])) /
                           (np.max(data['brand_price_max']) - np.min(data['brand_price_max'])))
data['brand_price_median'] = ((data['brand_price_median'] - np.min(data['brand_price_median'])) /
                              (np.max(data['brand_price_median']) - np.min(data['brand_price_median'])))
data['brand_price_min'] = ((data['brand_price_min'] - np.min(data['brand_price_min'])) /
                           (np.max(data['brand_price_min']) - np.min(data['brand_price_min'])))
data['brand_price_std'] = ((data['brand_price_std'] - np.min(data['brand_price_std'])) /
                           (np.max(data['brand_price_std']) - np.min(data['brand_price_std'])))
data['brand_price_sum'] = ((data['brand_price_sum'] - np.min(data['brand_price_sum'])) /
                           (np.max(data['brand_price_sum']) - np.min(data['brand_price_sum'])))



# 对类别特征进行 OneEncoder

data = pd.get_dummies(data, columns=['model', 'brand', 'bodyType', 'fuelType',
                                     'gearbox', 'notRepairedDamage', 'power_bin']) # 装换虚伪变量

print(data.shape)

print(data.columns)

# 这份数据可以给 LR 用
data.to_csv("data_for_lr.csv", index=0)

# 1）过滤式
# 相关性分析
print(data['power'].corr(data['price'], method='spearman')) #spearman：非线性的，非正太分析的数据的相关系数
print(data['kilometer'].corr(data['price'], method='spearman'))
print(data['brand_amount'].corr(data['price'], method='spearman'))
print(data['brand_price_average'].corr(data['price'], method='spearman'))
print(data['brand_price_max'].corr(data['price'], method='spearman'))
print(data['brand_price_median'].corr(data['price'], method='spearman'))

# 当然也可以直接看图
data_numeric = data[['power', 'kilometer', 'brand_amount', 'brand_price_average',
                     'brand_price_max', 'brand_price_median']]
correlation = data_numeric.corr() #返回data_numeric 相关性矩阵

f, ax = plt.subplots(figsize=(7,7))
plt.title("Correlation of Numeric Features with Price", y=1, size=16)
square=True # 将坐标轴方向设置为“equal”，以使每个单元格为方形 , vmax:色彩映射的值
sns.heatmap(correlation, square=True, vmax=0.8)
plt.show()

# 下一章介绍，Lasso 回归和决策树可以完成嵌入式特征选择
# 大部分情况下都是用嵌入式做特征筛选

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from operator import itemgetter

# %matplotlib inline 在终端中可以代替plt。show() 直接生成图

Train_data = pd.read_csv("./datalab/used_car_train_20200313.csv", sep=" ")
Test_data = pd.read_csv("./datalab/used_car_testA_20200313.csv", sep=" ")

# print(Train_data.shape)
# print(Train_data.head())

#print(Train_data.columns)

## 删除异常值
# 这里我包装了一个异常值处理的代码，可以随便调用
def outliers_proc(data, col_name, scale=3):
    """
    用于清洗异常值，默认用box_plot(scale=3)进行清洗
    :param data:接受 pandas 数据格式
    :param col_name:pandas 列名
    :param scale:尺度
    """
    def box_plot_outliers(data_ser, box_scale):
        """
        利用箱线图去除异常值
        :param data_ser:接收 pandas.Series 数据格式
        :param box_scale: 箱线图尺度 （规定大于上四分位数1.5倍四分位数差 的值，或者小于下四分位数1.5倍四分位数差的值，划为异常值）
        """
        iqr = box_scale * (data_ser.quantile(0.75) - data_ser.quantile(0.25)) # 3倍四分位数差
        val_low = data_ser.quantile(0.25) - iqr # 下限=Q1-3IQR
        val_up = data_ser.quantile(0.75) + iqr # 上限=Q3+3IQR
        rule_low = (data_ser < val_low)
        rule_up = (data_ser > val_up) # 返回 pandas.Series 中对应值的bool
        return (rule_low, rule_up), (val_low, val_up)

    data_n = data.copy() # copy 数据
    data_series = data_n[col_name] # 返回指定 col_name 数据
    rule, value = box_plot_outliers(data_series, box_scale=scale)
    index = np.arange(data_series.shape[0])[rule[0]|rule[1]] # 返回rule_low, rule_up中为True的下标的列表
    #print("Delete number is:{}".format(len(index))) # 打印下标列表中个数
    data_n = data_n.drop(index) # 删除(删除后下标没变)
    data_n.reset_index(drop=True, inplace=True) # 重置索引（drop=True删除原来的索引;inplace=True当前修改状态应用到原来Series中）
    #print("Now column number is:{}".format(data_n.shape[0])) # 查看删除后的数据个数
    index_low = np.arange(data_series.shape[0])[rule[0]]
    outliers = data_series.iloc[index_low] # ilco-按下标进行索引
    #print("Description of data larger than the lower bound is:")
    #print(pd.Series(outliers).describe())
    index_up = np.arange(data_series.shape[0])[rule[1]]
    outliers = data_series.iloc[index_up]
    #print("Description of data larger than the upper bound is:")
    #print(pd.Series(outliers).describe())

    #fig, ax = plt.subplots(1, 2, figsize=(10, 7)) # 创建子图:1行2列
    #sns.boxplot(y=data[col_name], data=data, palette="Set1", ax=ax[0]) # 箱线图
    #sns.boxplot(y=data_n[col_name], data=data_n, palette="Set1", ax=ax[1])
    #plt.show()
    return data_n


## 我们可以删掉一些异常数据,以 power 为例
## 这里删不删可以自行判断
## 但是注意 Test 的数据不能删
Train_data = outliers_proc(Train_data, "power", scale=3)



## 特征构造
# 训练集和测试集放在一起，方便构造特征
Train_data["train"] = 1 # 添加新字段，并设置值为1
Test_data["train"] = 1
data = pd.concat([Train_data,Test_data],ignore_index=True) # 连接函数 ignore_index=True重置索引

# 使用时间：data["createDate"] - data["regDate"], 反应汽车使用时间，一般来说价格与使用时间成反比
# 不过要注意, 数据里有时间出错的格式, 所以我们需要 errors = "coerce"
data["used_time"] = (pd.to_datetime(data["creatDate"], format="%Y%m%d", errors="coerce") -
pd.to_datetime(data["regDate"],format="%Y%m%d",errors="coerce")).dt.days # to_datetime将参数转换为日期 dt.days每个元素的天数

# 看一下空数据, 有 15k 个样本的时间有问题的, 我们可以选择删除, 也可以选择放着
# 但是这里不建议删除, 因为删除缺失数据占总样本量过大, 7.5%
# 我们可以先放着, 因为如果我们 XGBoost 之类的决策树, 其本身就能处理缺失值, 所以可以不用管
#print(data["used_time"].isnull().sum())

# 从邮编中提取城市信息, 相当于加入了先验知识
#print(data["regionCode"])
# 增加city 字段, 并从 regionCode 值的倒数第三位切片(apply 对regionCode每个元素运行指定运算 lambda 匿名函数)
data["city"] = data["regionCode"].apply(lambda x : str(x)[:-3])

# 计算某品牌的销售统计量, 还可以计算其他特征的统计量
# 这里以 train 的数据计算统计量
Train_gb = Train_data.groupby("brand") # 分组
all_info = {}
for kind, kind_data in Train_gb:
    info = {}
    kind_data = kind_data[kind_data["price"] > 0]
    info["brand_amount"] = len(kind_data)
    info["brand_price_max"] = kind_data.price.max()
    info["brand_price_median"] = kind_data.price.median()
    info["brand_price_min"] = kind_data.price.min()
    info["brand_price_sum"] = kind_data.price.sum()
    info["brand_price_std"] = kind_data.price.std() # 样本方差
    info["brand_price_average"] = round(kind_data.price.sum() / (len(kind_data)+1), 2) # round(2)取近似值保留两位数
    all_info[kind] = info

brand_fe = pd.DataFrame(all_info).T.reset_index().rename(columns={"index":"brand"}) # T转置   reset_index还原索引  rename并修改列名
data = data.merge(brand_fe, how="left", on="brand") # 合并数据 data链接在brand_fe "brand"字段左边




# 数据分桶 以 power 为例
# 这时候缺失值也进桶了
# 为什么要分桶：
# 1. 离散后稀疏向量内积乘法运算速度更快, 计算结果也方便存储, 容易扩展
# 2. 离散后的特征对异常值更具鲁棒性, 如 age>30 为 1 否则为 0 , 对于年龄为 200 的也不会对模型造成很大的干扰
# 3. LR 属于广义线性模型, 表达能力有限, 经过离散化后, 每个变量有单独的权重, 这相当于引入了非线性, 能够提升模型的表达能力, 加大拟合
# 4. 离散后特征可以进行特征交叉, 提升表达能力, 由 M+N 个变量变成　Ｍ*N 个变量, 进一步引入非线性, 提升了表达能力
# 5. 特征离散后模型更稳定, 如用户年龄区间, 不会因为用户年龄长了一岁就变化
# 当然还有很多原因,　LightGBM 在改进 XGBoost 时就增加了数据分桶, 增强了模型的泛化性

bin = [i*10 for i in range(31)]
data["power_bin"] = pd.cut(data["power"], bin, labels=False) # 分桶 cut切分数据(必须是一维的) bin定义区间 labels=False返回第几个bin（从0开始）
#print(data[["power_bin", "power"]].head())

# 删除不需要的数据
data = data.drop(["creatDate", "regDate", "regionCode"], axis=1) # drop函数默认删除行，列需要加axis = 1
#print(data.shape)
#print(data.columns)

# 目前的数据其实已经可以给树模型使用了, 所以我们导出一下
#data.to_csv("data_for_tree.csv", index=0) # index=0不保存行索引


# 我们可以再构造一份特征给 LR NN 之类的模型用
# 之所以分开构造是因为, 不同模型对数据的要求不同
# 先看下数据分布：
#data["power"].plot.hist()
#plt.show()

# 我们刚刚已经对 train 进行异常值处理了，但是现在还有这么奇怪的分布是因为 test 中的 power 异常值，
# 所以我们其实刚刚 train 中的 power 异常值不删为好，可以用长尾分布截断来代替
#Train_data["power"].plot.hist()
#plt.show()

# 我们对其取 log, 再做归一化
from sklearn import preprocessing
# 将数据的每一个特征缩放到给定的范围，将数据的每一个属性值减去其最小值，然后除以其极差（最大值 - 最小值）
min_max_scaler = preprocessing.MinMaxScaler()
data["power"] = np.log(data["power"] + 1)
# 归一化：(0,1)标准化
data["power"] = ((data["power"] - np.min(data["power"])) / (np.max(data["power"]) - np.min(data["power"])))
#data["power"].plot.hist()
#plt.show()

# km 的比较正常, 应该已经做过分桶了
# data["kilometer"].plot.hist()
# plt.show()

# 所以可以直接作归一化
data["kilometer"] = ((data["kilometer"] - np.min(data["kilometer"])) /
                     (np.max(data["kilometer"]) - np.min(data["kilometer"])))
#data["kilometer"].plot.hist()
#plt.show()

# 除此之外 还有我们刚刚构造的统计量特征：
# 'brand_amount', 'brand_price_average', 'brand_price_max',
# 'brand_price_median', 'brand_price_min', 'brand_price_std',
# 'brand_price_sum'
# 这里不再一一举例分析了，直接做变换，
def max_min(x):
    return (x - np.min(x)) / (np.max(x) - np.min(x))

data['brand_amount'] = ((data['brand_amount'] - np.min(data['brand_amount'])) /
                        (np.max(data['brand_amount']) - np.min(data['brand_amount'])))
data['brand_price_average'] = ((data['brand_price_average'] - np.min(data['brand_price_average'])) /
                               (np.max(data['brand_price_average']) - np.min(data['brand_price_average'])))
data['brand_price_max'] = ((data['brand_price_max'] - np.min(data['brand_price_max'])) /
                           (np.max(data['brand_price_max']) - np.min(data['brand_price_max'])))
data['brand_price_median'] = ((data['brand_price_median'] - np.min(data['brand_price_median'])) /
                              (np.max(data['brand_price_median']) - np.min(data['brand_price_median'])))
data['brand_price_min'] = ((data['brand_price_min'] - np.min(data['brand_price_min'])) /
                           (np.max(data['brand_price_min']) - np.min(data['brand_price_min'])))
data['brand_price_std'] = ((data['brand_price_std'] - np.min(data['brand_price_std'])) /
                           (np.max(data['brand_price_std']) - np.min(data['brand_price_std'])))
data['brand_price_sum'] = ((data['brand_price_sum'] - np.min(data['brand_price_sum'])) /
                           (np.max(data['brand_price_sum']) - np.min(data['brand_price_sum'])))



# 对类别特征进行 OneEncoder

data = pd.get_dummies(data, columns=['model', 'brand', 'bodyType', 'fuelType',
                                     'gearbox', 'notRepairedDamage', 'power_bin']) # 装换虚伪变量

#print(data.shape)
#print(data.columns)

# 这份数据可以给 LR 用
#data.to_csv("data_for_lr.csv", index=0)






# 特征筛选
# 1）过滤式
# 相关性分析
# print(data['power'].corr(data['price'], method='spearman')) #spearman：非线性的，非正太分析的数据的相关系数
# print(data['kilometer'].corr(data['price'], method='spearman'))
# print(data['brand_amount'].corr(data['price'], method='spearman'))
# print(data['brand_price_average'].corr(data['price'], method='spearman'))
# print(data['brand_price_max'].corr(data['price'], method='spearman'))
# print(data['brand_price_median'].corr(data['price'], method='spearman'))

# 当然也可以直接看图
# data_numeric = data[['power', 'kilometer', 'brand_amount', 'brand_price_average',
#                      'brand_price_max', 'brand_price_median']]
# correlation = data_numeric.corr() #返回data_numeric 相关性矩阵
#
# f, ax = plt.subplots(figsize=(7,7))
# plt.title("Correlation of Numeric Features with Price", y=1, size=16)
# square=True # 将坐标轴方向设置为“equal”，以使每个单元格为方形 , vmax:色彩映射的值
# sns.heatmap(correlation, square=True, vmax=0.8)
# plt.show()



# # 2)包裹式
from mlxtend.feature_selection import SequentialFeatureSelector as SFS #序列特征算法的实现——贪婪搜索算法
from sklearn.linear_model import  LinearRegression # 基于最小二乘法的线性回归
sfs = SFS(LinearRegression(), # 分类器或回归矩阵
          k_features=10, # 要选择的特征数量
          forward=True, #  如果为True，则向前选择，否则为反向选择
          floating=False, # 如果为True，则添加条件排除/包含。
          scoring="r2", # 对于sklearn回归变量使用“ r2”
          cv=0) # 如果cv为None、False或0，则不进行交叉验证
x = data.drop(["price"], axis=1)
x = x.fillna(0)
y = data["price"]

sfs.fit(x, y) # 执行特征选择并从训练数据中学习模型 x训练样本 y目标值
sfs.k_feature_names_


# 画出来，可以看到边际效益
from mlxtend.plotting import plot_sequential_feature_selection as plot_sfs
import matplotlib.pyplot as plt
fig1 = plot_sfs(sfs.get_metric_dict(), kind='std_dev')
plt.grid()

# 导入warnings包，利用过滤器来实现忽略警告语句
import warnings
warnings.filterwarnings("ignore")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import missingno as msno

## pd.set_option('display.max_columns', None)# 显示所有列
## pd.set_option('display.max_row', None)# 显示所有行
## 1)载入训练集和测试集
Train_data = pd.read_csv("./datalab/used_car_train_20200313.csv", sep = " ")
Test_data = pd.read_csv("./datalab/used_car_testA_20200313.csv", sep = " ")

print(Train_data.shape)

## 3)通过describe()来熟悉相关统计量
print(Train_data.describe())

## 4)通过info()来熟悉数据类型
print(Train_data.info())

## 5) 查看每列的存在nan情况
print(Train_data.isnull().sum())

#nan可视化
missing = Train_data.isnull().sum()
missing = missing[missing > 0]
missing.sort_values(inplace=True) # 排序
missing.plot.bar() # 绘柱状图
plt.tight_layout() # 自动调整子图参数
plt.show()

## 可视化看下缺省值
msno.matrix(Train_data.sample(250))
plt.show()

msno.bar(Train_data.sample(1000)) # 条形图
plt.show()

Train_data["notRepairedDamage"].replace("-", np.nan, inplace=True) # 将数据中‘-’替换成nan值
print(Train_data["notRepairedDamage"].value_counts())

print(Train_data.isnull().sum())

print(Train_data["seller"].value_counts())

print(Train_data["offerType"].value_counts())

# 删除严重倾斜的数据
del Train_data["seller"]
del Train_data["offerType"]
print(Train_data.info())
print(Train_data.shape)

print(Train_data["price"])
print(Train_data["price"].value_counts())

## 1)总体分布情况(无界约翰逊分布等）
import scipy.stats as st
y = Train_data["price"]
plt.figure(1); plt.title("Johnson SU") # 创建新图
sns.distplot(y, kde=False, fit=st.johnsonsu)
plt.figure(2); plt.title("Normal")
sns.distplot(y, kde=False, fit=st.norm)
plt.figure(3); plt.title("Log Normal")
sns.distplot(y, kde=False, fit=st.lognorm)
plt.show() # 最佳拟合是无界约翰逊分布

## 2)查看skewness and kurtosis
sns.distplot(Train_data["price"])
print("Skewness: %f" % Train_data["price"].skew()) # 偏度
print("Kurtosis: %f" % Train_data["price"].kurt()) # 峰度
plt.show()

print(Train_data.skew())
print(Train_data.kurt())

sns.distplot(Train_data.skew(), color="blue", axlabel="Skewness")
plt.show()

sns.distplot(Train_data.kurt(), color="orange", axlabel="Kurtness")
plt.show()

# 3)查看预测值的具体频数
plt.hist(Train_data["price"], orientation="vertical", histtype="bar", color="red")
plt.show() # 直方图

# log变换之后的分布比较均匀，可以进行log变换进行预测，这也是预测问题常用的trick
plt.hist(np.log(Train_data["price"]), orientation="vertical", histtype="bar", color="red")
plt.show()

# 分离label即预测值
Y_train = Train_data['price']

# 这个区别方式适用于没有直接label coding的数据
# 这里不适用，需要人为根据实际含义来区分
# 数字特征
# numeric_features = Train_data.select_dtypes(include=[np.number])
# numeric_features.columns
# # 类型特征
# categorical_features = Train_data.select_dtypes(include=[np.object])
# categorical_features.columns

# 数字特征
#numeric_features = ['power', 'kilometer', 'v_0', 'v_1', 'v_2', 'v_3', 'v_4', 'v_5', 'v_6', 'v_7', 'v_8', 'v_9', 'v_10', 'v_11', 'v_12', 'v_13','v_14' ]
# 类别特征
#categorical_features = ['name', 'model', 'brand', 'bodyType', 'fuelType', 'gearbox', 'notRepairedDamage', 'regionCode']

## 类别特征nunique分布——Train_data
for cat_fea in categorical_features:
    print(cat_fea+"的特征分布如下：")
    print("{}特征有{}个不同的值".format(cat_fea, Train_data[cat_fea].nunique()))
    print(Train_data[cat_fea].value_counts())

numeric_features.append("price")
print(numeric_features)

print(Train_data.head())

## 1)相关性分析
price_numeric = Train_data[numeric_features]
correlation = price_numeric.corr() # 返回一个相关系数的矩阵
print(correlation["price"].sort_values(ascending=False),"\n") # 降序排序

f , ax = plt.subplots(figsize = (7, 7))
plt.title("Correlation of Numeric Features with Price")
sns.heatmap(correlation, square=True, vmax=0.8) # 热图（显示相关系数）
plt.show()

## 2)查看几个特征的偏度和峰度
for col in numeric_features:
    print("{:15}".format(col),"Skewness:{:05.2f}".format(Train_data[col].skew()),
    "   ",
          "Kurtosis:{:06.2f}".format(Train_data[col].kurt()))

## 3)每个数字特征得分布可视化
f = pd.melt(Train_data, value_vars=numeric_features) # 转换
g = sns.FacetGrid(f,col="variable", col_wrap=2, sharex=False,sharey=False) # 以”variable“作“格子"绘图
# plt.show()
g = g.map(sns.distplot, "value") # 以”value“绘制到”格子”图中
plt.show()

## 4）数字特征相互之间的关系可视化
sns.set() # 风格设置
colunms = ["price", "v_12", "v_8", "v_0", "power", "v_5", "v_2", "v_6", "v_1", "v_14"]
sns.pairplot(Train_data[colunms],size=2, kind="scatter", diag_kind="kde") # 多变量图
plt.show()

print(Y_train)

## 5)多变量互相关系回归关系可视化
fig,((ax1, ax2), (ax3, ax4), (ax5, ax6), (ax7, ax8), (ax9, ax10)) = plt.subplots(nrows=5, ncols=2, figsize=(24, 20)) # 生成5行2列十个子图
# ['v_12', 'v_8' , 'v_0', 'power', 'v_5',  'v_2', 'v_6', 'v_1', 'v_14']
v_12_scatter_plot = pd.concat([Y_train,Train_data["v_12"]], axis=1) # 合并成一列
#print(v_12_scatter_plot)
sns.regplot(x="v_12", y="price", data=v_12_scatter_plot,scatter=True,fit_reg=True,ax=ax1) # 数据与回归模型拟合

v_8_scatter_plot = pd.concat([Y_train,Train_data['v_8']],axis = 1)
sns.regplot(x='v_8',y = 'price',data = v_8_scatter_plot,scatter= True, fit_reg=True, ax=ax2)

v_0_scatter_plot = pd.concat([Y_train,Train_data['v_0']],axis = 1)
sns.regplot(x='v_0',y = 'price',data = v_0_scatter_plot,scatter= True, fit_reg=True, ax=ax3)

power_scatter_plot = pd.concat([Y_train,Train_data['power']],axis = 1)
sns.regplot(x='power',y = 'price',data = power_scatter_plot,scatter= True, fit_reg=True, ax=ax4)

v_5_scatter_plot = pd.concat([Y_train,Train_data['v_5']],axis = 1)
sns.regplot(x='v_5',y = 'price',data = v_5_scatter_plot,scatter= True, fit_reg=True, ax=ax5)

v_2_scatter_plot = pd.concat([Y_train,Train_data['v_2']],axis = 1)
sns.regplot(x='v_2',y = 'price',data = v_2_scatter_plot,scatter= True, fit_reg=True, ax=ax6)

v_6_scatter_plot = pd.concat([Y_train,Train_data['v_6']],axis = 1)
sns.regplot(x='v_6',y = 'price',data = v_6_scatter_plot,scatter= True, fit_reg=True, ax=ax7)

v_1_scatter_plot = pd.concat([Y_train,Train_data['v_1']],axis = 1)
sns.regplot(x='v_1',y = 'price',data = v_1_scatter_plot,scatter= True, fit_reg=True, ax=ax8)

v_14_scatter_plot = pd.concat([Y_train,Train_data['v_14']],axis = 1)
sns.regplot(x='v_14',y = 'price',data = v_14_scatter_plot,scatter= True, fit_reg=True, ax=ax9)

v_13_scatter_plot = pd.concat([Y_train,Train_data['v_13']],axis = 1)
sns.regplot(x='v_13',y = 'price',data = v_13_scatter_plot,scatter= True, fit_reg=True, ax=ax10)
plt.show()

## 1）nunique分布
for fea in categorical_features:
    print(Train_data[fea].nunique())

print(categorical_features)

## 2)类别箱形图可视化
# 因为 name和 regionCode的类别太稀疏了，这里我们把不稀疏的几类画一下
categorical_features = ["model",
                        "brand",
                        "bodyType",
                        "fuelType",
                        "gearbox",
                        "notRepairedDamage"]
for c in categorical_features:
    Train_data[c] = Train_data[c].astype("category") # 强制转换数据类型
    if Train_data[c].isnull().any(): # 检查字段缺失
        Train_data[c] = Train_data[c].cat.add_categories(["MISSING"]) # 添加新类别
        Train_data[c] = Train_data[c].fillna("MISSING") # 填充为NAN的值
def boxplot(x, y, **kwargs):
    sns.boxplot(x=x, y=y) # 箱形图
    x=plt.xticks(rotation=90) #  设置坐标轴

f = pd.melt(Train_data, id_vars=["price"], value_vars=categorical_features)
g = sns.FacetGrid(f,col="variable", col_wrap=2, sharex=False,sharey=False,size=5)
g = g.map(boxplot, "value", "price")
plt.show()

print(Train_data.columns)

## 3)类别特征的小提琴图可视化
catg_list = categorical_features
target = "price"
for catg in catg_list:
    sns.violinplot(x=catg,y=target,data=Train_data)
    plt.show()

print(categorical_features)

## 4)类别特征的柱形图可视化
def bar_plot(x,y,**kwargs): # 柱形图
    sns.barplot(x=x,y=y)
    x=plt.xticks(rotation=90)
f = pd.melt(Train_data, id_vars=["price"], value_vars=categorical_features)
g = sns.FacetGrid(f, col="variable",col_wrap=2,sharex=False,sharey=False,size=5)
g = g.map(bar_plot, "value", "price")
plt.show()

## 5)类别特征的每个类别频数可视化
def count_plot(x,**kwargs): # 计数直方图
    sns.countplot(x=x)
    x=plt.xticks(rotation=90)
f = pd.melt(Train_data,value_vars=categorical_features)
g = sns.FacetGrid(f,col="variable", col_wrap=2,sharex=False,sharey=False,size=5)
g = g.map(count_plot,"value")
plt.show()

## 生成数据报告
import pandas_profiling

pfr = pandas_profiling.ProfileReport(Train_data)
pfr.to_file("./example.html")