神经网络是当前机器学习领域普遍所应用的，例如可利用神经网络进行图像识别、语音识别等，从而将其拓展应用于自动驾驶汽车等领域。神经网络的衍生变种目前有很多种，如CNN、RNN、GAN等，它们在不同应用场景有着各自针对性。但最简单且原汁原味的神经网络则是多层感知器（Muti－Layer Perception ，MLP）也叫人工神经网络(ANN,Artificial Neural Network)，除了输入输出层，它中间可以有多个隐层，最简单的MLP只含一个隐层，即三层的结构。只有理解经典的原版，才能更好的去理解功能更加强大的现代变种。

为了更好的理解，我画了手绘稿作为插图配合最后的代码来梳理人工神经网络的基本过程。

在理解神经网络前，不得不先引申出一点就是我们人类大脑中的基本单元———神经元。

虽然神经元有各种形式，但是所有的神经元大致都是将电信号从一端传输到另一端，然后沿着轴突把信号从树突传到树突。最后这些信号从一个神经元传送到另一个神经元。比方说我们的身体可以感知光、触感，声音等信号，就是来自感官神经元的信号最终传到了我们的大脑，而大脑也是由各种神经元组成的。

而人工神经网络是在理解和抽象了人脑结构及其对外界刺激的响应机制后，以生物神经系统的基本原理和结构为范本原型，以网络拓扑为理论基础，对复杂信息进行非线性关系表示和逻辑操作的一种数学模型。这里要说明的是并没有任何科学证明神经网络就是人脑结构的翻版，只能说是借鉴于人们对人脑结构的了解而启发创造的一种人工神经网络模型。

2、人工神经网络大致流程

手写数字识别被称为神经网络领域的 "Hello World"，因为识别人手写笔记的这个问题是检验人工智能基本的理想挑战，要让计算机准确区分图像中包含的内容。所以我把对MNIST的图像识别视为进入人工神经网络的大门。

在开始讲解代码实现神经网络之前，让我们先梳理一下神经网络流程的大致过程。首先，数据集是构建算法模型最重要的起点。如下图所示整个人工神经网络的过程大致可以分为这几步，导入数据集，并且对数据集进行必要的转换以此达到符合模型所需的矩阵数据格式；接着用划分的训练集对神经网络模型进行训练，当模型的准确度训练到一定程度后，最终就能通过训练好的神经网络模型进行图像识别查验了。

3、准备MNIST数据集

MNIST 数据集是著名的手写图像集合，常用来测量和比较机器学习算法的性能。该数据集包括了用于训练机器学习模型的 60,000 个图像和用于测试性能的 10,000 个图像。现成的MNIST的数据集主要是CSV文件，它的格式如下图分解图所示。每一行由一连串数字分别以逗号隔开。其中第一个值是数据标签，代表图像的实际"数字"，比如第一个值为 2那这个图像应该被识别为 2。随后的一连串数值由逗号分隔代表手写数字的组成像素值，由于像素的尺寸是28*28=784，所以除去首位标签值后，还有一共784个值。

4、神经网络的运行基本原理

如之前开头所说的那样，神经网络的运行过程大致可以分为接收导入的数据，然后进行训练，最后用输出结果误差对过程中的链接权重进行更新，最后生成模型。这中间训练的过程尤为重要，如下图所示我们的数据集是由784个数值组成的，这784个数值作为第一层输入层被分为了784个输入层节点，经过与下一层的链接权重进行一系列的组合计算并应用激活函数后完成第一遍输出值。接着通过这个输出值与目标值的比较得到误差值，最后使用误差值反向更新权重以此来调节链接权重的大小，反复几次操作后神经网络模型生成完毕。

4.1、神经网络的计算过程 --- 信号正向传播

神经网络信号的传播是由输入到最终输出的过程，输入的数值经过 输入层 到 隐藏层 之间的链接权重进行组合计算，最后应用激活函数来调节输出信号值。这样依次从输入层到最终输出层计算的过程称为信号的正向传播。

如上图所示为了更好的说明，以求得隐藏层第一个节点$J_ {1}$为例，首先在求任何层的节点数值时需要通过上一层（即输入层）的输入数值与对应的链接权重相乘来获得对应的隐藏层节点数值，这是一个组合计算的过程。可以发现隐藏层 $J_ {1}$ 节点的数值是由输入层 $I_ {1}$ 和 $I_ {2}$ 两个节点组成，所以分别将输入数值乘以与目标节点对应的链接权重：$I_ {1} \times W_ {1,1} + I_ {2} \times W_ {2,1}=J_ {1}$ 求得隐藏层 $J_ {1}$ 节点的数值，之后必须将数值应用于Sigmoid激活函数来调节输入信号。

想要求得更多的隐藏层节点输入信号就需要更多的输入层数值乘以对应链接权重，并最后应用激活函数求得，但是用这样的书写来计算所有节点过于麻烦，最方便的办法就是使用矩阵。通过上图可以发现，其实可以使用矩阵乘法表示所有组合计算，$X = {W}\cdot{I}$，这里$X$ 是输出矩阵， $W$ 是权重矩阵，$I$ 是输入矩阵。所以正向传播求得每一层节点输出值的过程就是 $Sigmoid(X) = {W}\cdot{I}$。

4.2、神经网络的计算过程 --- 误差反向向传播

正向传播的计算过程我们已经明白了，但是这样来判断一个结果来说还是远远不够的。因为所有传入的数据仅仅是通过组合计算在神经网络层上进行了一次正向传播而已，正由于初始的链接权重通常是随机值，基于这样的组合计算得到的输出值必然与目标值会有很大的误差，因此还需要通过这个误差值反向对链接权重进行优化（如误差过大，就调小权重，反之增大），这样做的目的是为了神经网络经过层层计算后得到的输出值更接近实际目标值。

为了更清晰描述反向传播过程，我将它分解成为两个细节部分，第一部分求得前一层（即隐藏层）误差的过程，第二部分为借助这个误差来更新权重的过程。

神经网络在进行反向传播前必然进行了一次正向传播，在正向传播的过程中经过输出层后会得到输出层的输出值，而通过目标值与这个输出值相减就可以得到输出层的实际误差值。而这个误差值还可以反向得到上一层（隐藏层）的误差值。你或许会问为什么还要计算出隐藏层的误差呢，因为我们需要使用隐藏的误差值来更新隐藏层到输出层之间的权重值。

但是反向计算上一层的误差值的方法可以有很多种，大致可以分为如下4种。

（1）直接将误差值以链接权重的个数N来平分，如输出层的误差值是由2个链接对应的节点组成，所以就以$\frac{1}{2}$分割误差值：

（2）采用以链接权重值的占比来计算：

这种方法的思想是不同的链接权重具有不同的权重强度，所以理应以权重值大小的比例来分割误差值。如下图$W_ {1,1}$的权重值为3.0，$W_ {2,1}$的权重值为1.0，那么$W_ {1,1}$的权重强度比重为$\frac{W_ {1,1}}{W_ {1,1}+W_ {2,1}}=\frac{3}{4}$，$W_ {2,1}$的权重强度比重为$\frac{W_ {2,1}}{W_ {1,1} + W_ {2,1}}=\frac{1}{4}$，也就是说上一层的节点$1$将获得$\frac{3}{4}$比例的误差值，而节点$2$将获得$\frac{1}{4}$比例的误差值，使用这种方式的好处是后续能更精确的优化权重数值，但是这也有个缺点，就是我们很难通过这种方式去进行大量的计算，即使使用了矩阵计算。

（3）直接用误差与链接权重相乘 $E_ {hidden}=W_ {i,j}^T⋅E_ {output}$（这也是本文使用的方法）。

以上图为例，反向计算 $(E)J_ {1} 节点的误差值是由 W_ {1,1} 和 W_ {1,2}$ 链接权重组成，所以计算该节点的过程就是$(E)J_ {1}=(E)K_ {1} \times W_ {1,1} + (E)K_ {2} \times W_ {1,2}$ 这种方法的好处是计算非常简单，即使反馈的误差值过大或过小，但在下一轮的迭代中，神经网路也可以自行纠正。

注意，由于是反向计算，链接权重排列的顺序与之前正向传播时的顺序不同，所以这里需要对链接权重的矩阵进行转置。

（4）与第3种方式相近，但是表达式额外除以了对应的链接权重个数N：

$E_ {hidden}=\frac{W_ {i,j}^T}{N}.E_ {output}$

这样做的方式有点类似某种形式的归一化，对于某些过大或者过小误差来说，可以缩小一定误差的比例范围。

为了更好的说明这四种方法实际计算的效果如何，如下统计图是对不同方式的性能计算展示。蓝色代表（第3种）使用权重值与输出层误差值直接相乘，红色代表（第4种）仅仅增加归一化的形式，黄色代表（第1种）直接将误差值以链接权重的个数N来均匀分割，绿色代表（第2种）以链接权重值的实际占比计算。

最后可以发现蓝色和绿色线形图在多次循环训练后都得到了不错的结果。而出于易于计算的考虑，最终选择直接将输出层的误差值与链接权重值相乘的这种方式来计算隐藏层的误差值。

这里不得不提到的一点是，当神经网络多于2层时，那么只需从最终输出层往回工作，重复应用相同的思路计算出前一层的误差值就可以了。

在获得了上一层（隐藏层）的误差值后，就可以真正的通过误差来调整对应的链接权重了。为了希望得到最小化的误差函数，我们试图要知道误差对链接权重的改变有多敏感，换个方式说就是当权重$W_ {ij}$ 这个自变量变化时，误差$E$这个因变量是如何改变的，所以我们需要做的是对链接权重求导，就是$\frac{\partial{E}}{\partial{W_ {jk}}}$ 。

如上求导的过程中使用了链式法则，而在链式法则的计算中需要代入前一层（隐藏层）的误差值，这也再次证明了我们为什么之前反复提到需要求得前一层的误差值是必须的了。通过推导可以得到如下两个链接权重的更新公式：

隐藏层与输出层的权重变化率： $\Delta{W_ {jk}}=-E_ {k}\times{Sigmoid(O_ {k})\times(1-Sigmoid(O_ {k}))\cdot{O^T_ {j}}}$

输入层与隐藏层的权重变化率： $\Delta{W_ {ij}}=-E_ {j}\times{Sigmoid(O_ {j})\times(1-Sigmoid(O_ {j}))\cdot{O^T_ {i}}}$

最后我们只需要把旧权重与新权重进行相减来更新权重：

$newW_ {jk}=oldW_ {jk}-\Delta{W_ {jk}}$

由于直接选择合适的权重太难了，另一种方法是通过误差函数的梯度下降来采取小步长，迭代的改进权重，也就是所谓的梯度下降法，为了使每次调整斜率以合适的步长幅度调整，可以引入learning rate学习率的概念，所以最终的公式可以演变为：

$newW_ {jk}=oldW_ {jk}+\alpha\times\Delta{W_ {jk}}$

至此这就是对于神经网络如何训练并自我更新权重的所有计算过程。

Python代码实现

接下来就是纯Python实现人工神经网络的所有代码，先导入最基本的模块库，这里除了最基本的numpy和pandas外，还因为后面需要更好的呈现识别图像的结果，所以还用到了matplotlib。另外由于要进行大量的数学运算所以不得不借助scipy这个科学计算包。由于神经网络训练是一个随机的过程，有时候工作的不错，有时候就很糟糕，所以为了更好的每次测试比对，这里设置了随机种子以保证每次的随机影响能保持一致。

import numpy as np  
import pandas as pd  
import matplotlib.pyplot as plt  
import scipy.special

from scipy import stats

np.random.seed(80)  

神经网络的整个结构其实很简单，如下代码由一个名为NeuralNetwork的类以及内部几个主要方法构成。主要包含3大方法：

_init_()方法：定义神经网络的基本参数，其中分别设置了输入层、隐藏层、输出层的节点数参数，以及输入层与隐藏层的权重值、隐藏层与输出层的权重值、学习率、激活函数。

train()方法：由于神经网络的训练模式组成是先通过正向传播得到基本误差，然后反向传播使用误差更新权重。所以在该方法里分别要实现正反向传播操作的过程。

query()方法：经过之前的一系列训练后，神经网络的连接权重值都已经优化到最佳数值，所以对于查询方法，只需要传入数据集进行一遍正向传播就可以输出结果。

### neural network class definition ###
class NeuralNetwork:  
    ### initialise the neural network elements
    def __init__(self, input_nodes, hidden_nodes, output_nodes, learning_rate):
        self.i_nodes = input_nodes
        self.h_nodes = hidden_nodes
        self.o_nodes = output_nodes
        self.lr = learning_rate

        self.wih = np.random.normal(0.0, pow(self.h_nodes, -0.5), (self.h_nodes, self.i_nodes))
        self.who = np.random.normal(0.0, pow(self.o_nodes, -0.5), (self.o_nodes, self.h_nodes))

        # activation function is the sigmoid function
        self.activation_function = lambda x: scipy.special.expit(x)


    ### training the neural network
    def train(self, inputs_list, targets_list):
        # convert the list to 2d array, because the each layer of data format must be 2D
        inputs = np.array(inputs_list, ndmin = 2).T
        targets = np.array(targets_list, ndmin = 2).T

        # calculate the signals into hidden layer
        hidden_inputs = np.dot(self.wih, inputs)
        # calculate the signals emerging from hidden layer
        hidden_outputs = self.activation_function(hidden_inputs)

        # calculate the signals into output layer
        final_inputs = np.dot(self.who, hidden_outputs)
        # calculate the signals emerging from final layer
        final_outputs = self.activation_function(final_inputs)

        # error is the (target - actual): 
        output_errors = targets - final_outputs

        # hidden layer error is the output_errors, split by weights, recombined at hidden nodes
        hidden_errors = np.dot(self.who.T, output_errors)

        # update the weights for the links between the hidden and output layers
        self.who += self.lr * np.dot((output_errors * final_outputs * (1.0 - final_outputs)), np.transpose(hidden_outputs))
        # update the weights for the links between the hidden and input layers
        self.wih += self.lr * np.dot((hidden_errors * hidden_outputs * (1.0 - hidden_outputs)), np.transpose(inputs))


    ### query the neural network
    def query(self, inputs_list):
        # covert inputs list to 2d array
        inputs = np.array(inputs_list, ndmin = 2).T

        # calculate signals into hidden layer
        hidden_inputs = np.dot(self.wih, inputs)
        # calculate signals emerging from hidden layer
        hidden_outputs = self.activation_function(hidden_inputs)

        # calculate signals into output layer
        final_inputs = np.dot(self.who, hidden_outputs)
        # calculate signals emerging from output layer
        final_outputs = self.activation_function(final_inputs)

        return final_outputs

如下代码是先创建神经网络对象，并这个对象设置基本参数，由于我们之前已经提到一个数字有784个像素数值构成，所以初始的输入节点自然为784个，而隐藏层的节点数通过一系列调参观察后，这里设置为200比较合适，而最后的输出数值设置为10的原因是因为我们要识别的阿拉伯数字为0~9十个数字，结果必然是10个节点中的某一个，最后的学习率设置为0.1也是通过调参观察而定的。最后设置了5个世代作为循环次数，是因为经过几次实验后发现使用5个世代循环训练在当前的过程中效果更佳，如果再增加下去就会出现过拟合的情况。

不得不提到的一点是由于数据的有限性，所以在已有数据的基础上代码最后还使用了rotate逆时针顺时针旋转的方式，以使现有数据产生变种来增加数据量，这样可以为模型增加训练强度。

### training the neural network ###
# loading the mnist traning data CSV file into list
training_data_file = open("neuralnetwork/mnist_dataset/mnist_train.csv", 'r')  
training_data_list = training_data_file.readlines()  
training_data_file.close()

# define the number of input, hidden, output nodes and learning rate
input_nodes = 784  
hidden_nodes = 200  
output_nodes = 10  
learning_rate = 0.1

# create the instance of neural network
n = NeuralNetwork(input_nodes, hidden_nodes, output_nodes, learning_rate)

# epochs is the number of times the training data set is used for training
epochs = 5

for e in range(epochs):  
    # go through all records in the traning data set
    for record in training_data_list:
        # split the record by the ',' commas
        all_values = record.split(',')

        # scale and shift the inputs
        inputs = (np.array(all_values[1:], dtype = np.float) / 255.0 * 0.99) + 0.01

        # create the target output values(because we use the sigmoid function, 
        # its range is 0<s<1, all 0.01, expect the desired label which is 0.99)
        targets = np.zeros(output_nodes) + 0.01

        # all_values[0] is the target label for this record
        targets[int(all_values[0])] = 0.99
        # training the model
        n.train(inputs, targets)

        ## create rotated variations in order to increase the sample capacity
        # rotated anticlockwise by x degrees
        inputs_plusx_img = scipy.ndimage.interpolation.rotate(inputs.reshape(28,28), 10, cval=0.01, order=1, reshape=False)
        # training the model
        n.train(inputs_plusx_img.reshape(784), targets)

        # rotated clockwise by x degrees
        inputs_minusx_img = scipy.ndimage.interpolation.rotate(inputs.reshape(28,28), -10, cval=0.01, order=1, reshape=False)
        # training the model
        n.train(inputs_minusx_img.reshape(784), targets)      

接下来使用测试集对神经网络模型进行测试并评分。

### test the neural network ###
# load the test data
test_data_file = open("neuralnetwork/mnist_dataset/mnist_test.csv", "r")  
test_data_list = test_data_file.readlines()  
test_data_file.close()


# scorecard for how well the network performs, initially empty
scorecard = []

# go through all the records in the test data set
for record in test_data_list:  
    # split the record by the ',' commas
    all_values = record.split(',')

    # convert answer is first label
    correct_label = int(all_values[0])

    inputs = (np.array(all_values[1:], dtype = np.float) / 255.0 * 0.99) + 0.01

    ## query the neural network
    outputs = n.query(inputs)

    # the index of the highest value corrsponds to the label
    predict_label = np.argmax(outputs)

    if (predict_label == correct_label):
        scorecard.append(1)
    else:
        scorecard.append(0)

# calculate the performance score, the fraction of correct answers
scorecard_array = np.array(scorecard)  
print(f'performance: {round(scorecard_array.sum() / scorecard_array.size * 100, 2)}%')  

测试评分结果为96.65%的正确率，这个结果还是十分不错的。

performance: 96.65%  

最后为了更好的展现，我对模型进行批量输入并进行可视化展示，一起来看看实际效果如何。

### test the custom image ###
from pathlib2 import Path  
import imageio  
import matplotlib.pyplot as plt


img_array_list = []

files_path = Path.cwd().joinpath('neuralnetwork/my_own_images')  
if files_path.exists():  
    for file in files_path.iterdir():
        if file.match('*.png'):
            img_array = imageio.imread(str(file), as_gray = True)
            img_array_list.append(255.0 - img_array)

fig, axes = plt.subplots(figsize=(6, 6), nrows = 3, ncols = 3)  
i = 0

for row in range(3):  
    for col in range(3):
        # reshape the 784 based on white ground
        img_data = img_array_list[i].reshape(784)
        # scale and shift the inputs
        inputs = (img_data / 255.0 * 0.99) + 0.01
        # query the neural network
        outputs = n.query(inputs)

        # the index of the highest value corresponds to the target label
        label = np.argmax(outputs)

        axes[row][col].imshow(img_array_list[i], cmap = 'Greys', interpolation = 'None')
        axes[row][col].set_title(f'predict:{label}')
        axes[row][col].set_xticks([])
        axes[row][col].set_yticks([])

        i += 1

如下图所示，从图示结果看预测成功率还不错，除了把6错误的预测为9外，其它均预测正确。我猜可能是由于背景有色差噪点的关系导致没有正确识别，加上先前训练模型的样本中可能也不存在这样的情况。

之前提到，带有噪声的数字6没有被正确识别，可能是由于训练模型时不存在类似的样本作为训练数据，所以导致判断错误。为了验证这点，我干脆把这个未识别的图像作为训练集单独拎出来再次训练下模型，看看模型最终能不能识别正确。

将未识别正确的数据作为样本训练了10次后，果真最后可以将带有噪声的数字6识别正确了。当然这种方式属于拿测试集当训练集的作弊行为了。

如下我先模拟一个数据集来说明如何解决这个问题，以DataFrame类型为数据集输出10行6列数据。

既然这篇文章要说明的是SettingWithCopyWarning的解决方案，那首先要知道SettingWithCopyWarning什么时候才会遇到，它是如何产生的。

在演示如何产生SettingWithCopyWarning警告前，我打算先说明一个知识点，那就是当对DataFrame数据进行操作的时候会出现浅拷贝（Shallow Copy）与深拷贝（Deep Copy）两种模式。下图描述了浅拷贝和深拷贝的区别，假设B复制了A，当修改B时，如果A也跟着变了，说明这是浅拷贝。如果A没变，那就是深拷贝。

说的简单点就是浅拷贝只是创建了对象的一个引用，而深拷贝则是创建了对象的一个独立的实体副本。

明白了深浅拷贝的区别后，来看如下例子，下图代码的意图是从数据中筛选出所有性别为男性的数据，后面会打算把这些数据设置为数字0。

注意，下面的代码打算把这些数据设置为数字0时因为使用多个中括号来索引操作，也就是显式链式赋值的方式。第一次是访问操作（get），返回一个 DataFrame列出了所有包含性别为男的数据，然后第二个是赋值操作（set）对gender列数据赋值为0，而这个赋值操作是在筛选出所有性别为男的这个新的 DataFrame 数据集上运行的，而压根没有在原始 DataFrame 数据集上运行。

这个时候pandas就分不清到底你是对数据做深拷贝还是浅拷贝操作，也就是修改的数据是本身还是独立的，更不知道最后得出来的结果是不是你想要的，简单理解，就是pandas无法对通过两个方括号选取的数据赋值。最后就产生模棱两可的状态于是就抛出了SettingWithCopyWarning警告。

可以看到，当赋值操作后再次输出data的数据后，发现gender数据里的男数据并没有被修改为0，说明pandas分不清到底你是对数据做深拷贝还是浅拷贝操作，代码执行的结果没有按预期运行。这也就是之前提到为什么遇到SettingWithCopyWarning警告时不能视而不见的原因了。

解决方案（一）

这个解决方案很简单：使用 loc 将链式操作组合到一个 方括号[ ] 操作中，这样用检索所条件作为loc的行列参数，将两次索引变成一次，以便 Pandas 可以确保 set 操作是在原始 DataFrame上执行。Pandas 会始终确保下面这样的非链式 set 操作起作用。其实这个方案也正是SettingWithCopyWarning里建议的操作方式。

解决方案（二）

另一种解决方案就是使用隐式赋值的方式。从原始数据中先提取自己想要的数据后并赋值保存给另一个变量，然后再对这个新变量对象进行操作。

可以从上面的结果看到，代码成功修改了数据值。但是依然出现了SettingWithCopyWarning警告信息，原因是通过隐式赋值变量的方式传递数据pandas不能保证修改的是本身数据。

所以，要解决这个警告信息的办法就是明确使用copy()方法，告诉它我创建的是一个独立副本。

这样就不会出现任何烦人的警告信息了。但是要注意这种办法并不是在原始的数据集上做修改，而是把需要修改的数据完全单独剥离了出来，所以最终输出的数据也仅仅是被筛选过的。

最后提一个中肯的建议，即使SettingWithCopyWarning警告只在 set 操作时才会发生，但在进行 get 操作时，最好也避免使用链式索引。不仅因为链式操作较慢，而且因为链式索引始终是一个潜在的问题，只要你稍后进行赋值操作，就可能不会影响到原始对象，这可能不是你的预期操作。所以，请不惜一切代价避免使用链式索引。

参考资料：https://www.dataquest.io/blog/settingwithcopywarning/

在numpy与pandas的使用中，有个常见的参数axis，根据对axis的设定值不同就会得到截然不同的结果。对于如何正确设置axis的参数值，如果有人与我曾经一样有似懂非懂的经历，那一定是在某方面没有正确的理解作者对这个参数的定义。

为了彻底说清楚axis到底是什么，我会用手稿图结合多个详细的例子来总结。

先借用网上常见描绘的axis图：

其实这个图里的方向指示很容易误导人，包括一些网上的解释也比较模糊，只要你站在不同角度去理解就会造成理解偏差，我曾经就误解绕了进去，虽然用自己理解的方法得到的结果看似是对的，但是放到其它场景时又发现是错的，才导致似懂非懂。

首先，axis定义为轴。在一维的数组中它只有一个轴就是0轴，指的是遍历每一行（所对应的是index），而把每一行拿出来遍历后，它所呈现的操作方向是垂直向下的。当在二维数组中它就有两个轴，除了前面提到的0轴外另一个轴就是1轴，指的遍历每一列（所对应的是columns），而把每一列拿出来遍历后，它所呈现的操作方向是横向水平的。

所以简单只需记住，设定axis是为了确定要删的标签是属于index还是column：

• axis = 0 时跨行操作，一行行，方向垂直向下。所以操作的是行标签，也就是index，遍历的操作方向是垂直向下。如果是以聚合操作（如sum、mean、apply等），指的是按照遍历的操作方向竖着，跨行操作。如果从直观上理解就是，常规形式下axis = 0是对整行操作，聚合形式下是对整列操作计算。
• axis = 1时跨列操作，一列列，方向水平横向。所以操作的是列标签，也就是columns。如果是以聚合操作（如sum、mean、apply等），指的是按照遍历的操作方向横着，跨列操作。如果从直观上理解就是，常规形式下axis = 1是对整列操作，聚合形式下是对整行操作计算。

在用例子说明前，先定义一个初始DataFrame数据类型，也就是前面提到的二维数据结构，所以它具有0轴和1轴。

1、当axis=0时移除缺失值

本例子的目的是删除含有存在缺失值的标签索引行，这里使用dropna()方法移除缺失值，由于设置了axis=0所以对每一行进行遍历（但是它遍历的方向是垂直向下的，这就是为什么说axis=0时方向向下），如果发现任何一行中包含缺失值，那么删除该行。

操作步骤图解如下：

经过上图的图解过程后，会发现0，1，2行都包含缺失值所以都被删除。最后只剩下索引为3的这一行，也就是之前代码运行得到的正确结果。

2、当axis=1时移除缺失值

本例子的目的是删除含有存在缺失值的标签列，这里使用dropna()方法移除缺失值，由于设置了axis=1所以对每一列进行遍历（但是它遍历的方向是水平横向的，这就是为什么说axis=1时方向水平），如果发现任何一列中包含缺失值，那么删除该列。

操作步骤图解如下：

经过上图的图解过程后，会发现col3，col4两列都包含缺失值所以都被删除。最后只剩下col1，col2的这两列，也就是之前代码运行得到的正确结果。

3、使用thresh参数搭配axis参数

当遍历的每一行或列，满足大于等于thresh参数所指定数量的非缺失值（non-NA），便显示这一行或列。例：如果thresh=4，当遍历的行或列中存在大于等于4个非缺失值，即可显示，否则移除。

如下图解，当thresh=4，axis=0时，使用dropna()方法，代表遍历每一行寻找是否满足4个非缺失值，不满足就移除。最后发现只剩下索引为3的这一行满足，所以只显示该行数据。

同理，当thresh=4，axis=1时，使用dropna()方法，代表遍历每一列寻找是否满足4个非缺失值，不满足就移除。最后发现只剩下标签列为col1和col2的这两列满足，所以只显示该两列数据。

4、使用聚合操作搭配axis参数

当使用聚合操作时和平时会有些不一样，比如当axis=0时是向下跨行，axis=1时是水平跨列。

怎么理解跨行和跨例呢，看下下面的例子：

首先，计算mean、sum等批量计算的函数操作时，可以理解为属于聚合操作。如下例子，当要计算当axis=0时的sum总和操作时，从图解中可以看到，它并不是按照之前理解的那样计算每一整行的总和，而是指的是按照遍历的操作方向向下跨行计算，输出每列的总和。

要计算当axis=1时的sum总和操作时，从图解中可以看到，它并不是按照之前理解的那样计算每一整列的总和，而是指的是按照遍历的操作方向水平跨列计算，输出每行的总和。

所以，使用聚合操作的时候，axis = 0 就是对整列操作，axis = 1就是对整行操作，正好和常规情况相反。

为了更好的以不同例子说明这几个方法的作用与区别，如下先声明一个初始数据集。

一、set_index() 的使用

set_index() 主要可以将数据表中指定的某列设置为索引或复合索引，如下是常涉及使用的几个参数：

• keys：列标签或列标签/数组列表，需要设置为索引的列。
• drop：默认为True，删除用作新索引的列，也就是当把某列设置为索引后，原来的列会移除。
• append：是否将列附加到现有索引，默认为False。
• inplace：输入布尔值，表示当前操作是否对原数据生效，默认为False。

如下代码说明，当keys参数设置为height字段时，height字段将变成索引列。同时因为drop=True 代表同时删除原数据中的height列，默认情况本身为True。

可以发现，当drop=False时，保留了原数据中的height列。

如下代码与之前的代码基本相似，也是将height字段将变成索引列。同时drop=False，代表保留原数据中的该列。但有明显不同之处是设置了append=True，这就代表将height列是以添加的额外附加的形式添加到现有索引中，但并不会把原来最初的索引移除。

可以看到正因为之前使用了append参数以附加索引的方式添加，如下当查看index索引时会发现得到的索引是复合索引。

二、reset_index() 的使用

reset_index() 主要可以将数据表中的索引还原为普通列并重新变为默认的整型索引。如下是常涉及使用的几个参数：

• level：数值类型可以为：int、str、tuple或list，默认无，仅从索引中删除给定级别。默认情况下移除所有级别。控制了具体要还原的那个等级的索引 。
• drop：当指定drop=False（默认为False）时，则索引列会被还原为普通列；否则，如设置为True，原索引列被会丢弃。
• inplace：输入布尔值，表示当前操作是否对原数据生效，默认为False。

如下代码我先通过set_index() 方法将字段age列数据设置为索引，并赋值为新变量df_new，之后会使用reset_index() 方法来说明该方法的作用。

如下，现在针对df_new对象变量使用reset_index()方法，同时设置drop=False（默认为False）代表原来的索引列age会保留并被还原为普通列，同时索引列变为默认的整型索引。

可以发现当drop=True时，除了通过reset_index()重新添加了整型索引后，age索引列已经被移除。

• reset_index() 案例场景说明

理解了方法固然也要明白什么场景下使用，这样才能便于更好的记住。通常在做数据预处理的时候，会对数据表的缺失值或者其它不符合预期的值做移除，这个时候很容易会遇到索引值不在连贯，那么就需要重新设置索引排序。

如下例子我先模拟生成了一个5行5列的数据表，它的默认索引为整型索引。因为一些需求通过drop() 方法移除了索引行为2、4的数据。可以从下图看到处理过后的数据只剩下0，1，3索引行数据。

因为经过处理后索引值的排序已经不连贯。为了使整形索引排序连贯，可以回忆下之前的reset_index() 方法可以重新还原整型索引。所以如下代码使用reset_index() 方法很顺利的添加了连贯的整型索引，也正因为默认参数drop=False，所以原来的索引值[0, 1, 3]被还原成了普通列，同时被系统冠名为index字段。

但是通常我们会选择移除该列，所以就像下图那样可以显示设置参数drop = True 进行移除。

但是千万别忘了，任何移除操作并不会真正修改数据对象本身，如需要一定要设置inplace参数。

三、reindex() 的使用

reindex() 顾名思义它的作用是用来重定义索引的，如果定义的索引没有匹配的数据，默认将已缺失值填充。而索引可以分 “行” 索引与 “列” 索引，所以reindex自然对于两者的修改都可以胜任。它在Series和DataFrame中都非常有用：

对于DataFrame，reindex() 可以修改行、列索引或者两个都修改。

对于Series，reindex() 会创建一个适应新索引的新对象，如果某个索引值当前不存在，就会引入缺失值。

另外，对于以上两个数据类型都可以通过fill_value参数填充默认值，也可以通过method参数设置填充方法。而method包含几个参数可以选择：

• None (默认): 不做任何填充
• pad / ffill: 用上一行的有效数据来填充。
• backfill / bfill: 用下一行的有效数据来填充。
• nearest: 用临近行的有效数据来填充。

Series中的reindex() 使用

如下例子，先讲解下在Series中reindex() 方法的使用。同样，为了更好的说明，先使用Series初始创建数据，该数据默认的索引为5个字母a, b, c, d, e。

接着，使用reindex() 方法对数据的索引进行重定义为a, b, c, d, e, f。因为f之前并不存在，是我们新定义的，所以也没有匹配的数据值，那么默认就以缺失值NaN填充。当然也可以像下图最后一段代码那样使用 fill_value 参数设定填充数值，这里我设置为0.0。

也可以使用设置method参数来填充，这里使用ffill参数值设定使用上一行的数据来作为填充数据。

DataFrame中的reindex() 使用

同样，为了更好的说明，先使用DataFrame初始创建数据，该数据默认的索引为5个字母a, b, c, d, e。

之前提到过在DataFrame中reindex() 可以修改行、列索引。所以如下代码分别对index行索引和columns列索引进行了设置。

当然也可以同时指定index和columns来重定义索引，因为定义的索引既没有匹配的数据也没设置默认填充，所以就以缺失值填充。

当然不仅可以通过重定义索引来增加新索引，反过来也可以做一些移除索引达到drop的效果。

可以看到如在Series中操作reindex() 一样，因为重定义的索引没有找到匹配的对象数据，所以默认填充为缺失值。而在DataFrame中依然可以使用fill_value与method两个参数方法来设定默认填充值。

method参数只适用于index是单调递增或者单调递减的情形。同时引用一张网络图很好的说明了method参数的使用逻辑。

蒙特卡络不是一个人名，而是一个地名，因摩纳哥著名的赌场而得名，而该方法的提出者是大名鼎鼎的数学家冯·诺伊曼（现代计算机之父）。

蒙特卡洛(Monte Carlo)方法，又称为随机抽样或统计试验方法，是以概率和统计理论方法为基础的一种计算方法，本质是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。它将所求解的问题同一定的概率模型相联系，以获得问题的近似解。

这里要说明的一点是蒙特卡络方法是一种基于概率方法的统称，包含蒙特卡洛算法、模拟、过程、搜索树等。而且使用可能得到的数值并非是最终最精确的那个，会有一定的误差，而误差的大小与模拟的样本大小直接相关，模拟样本越大，误差越小，它是一种概率数值的逼近。对于简单问题来说，蒙特卡洛是个“笨”办法。但对许多问题来说，它往往是个有效，有时甚至是唯一可行的方法。

用蒙特卡洛方法求圆周率π值

考虑到π和圆形是紧密相关的，所以我首先构造一个边长为1的正方形和半径为1的圆形。为什么选正方形呢，因为正方形的边长相等，又和圆形半径类似，而且正方形在某个角度可以看做是平面坐标系，方便后面计算。

然后使两个图案相交，重叠部分是一个$\frac{1}{4}$个单位圆部分（将这个重叠部分简称为P(A)）。

接着，在边长为1的正方形区域内随机落点，统计落在重叠部分P(A)区域内点的所占比例。这个比例其实很好计算，因为随着随机落点数量的增加，每个随机点的落点概率会越来越接近四分之一圆面积与正方形面积之比（概率比值即为$\frac{π}{4}$），而这个概率就是随机落点到重叠部分P(A)的比例。

π值求解的推导公式如下：

$P(A)=\frac{四分之一圆面积}{正方形面积}=\frac{\frac{1}{4}πr^2}{r^2}=\frac{π}{4}$

$P(A) = \frac{π}{4}$

$π = P(A)\times4$

考虑到用文字难以形象的表达清楚计算π的过程，下面用我的手绘稿来讲解下原理。

注：由于在我的手稿中以平面坐标系来看，将$\frac{1}{4}$圆形放在第二象限，所以x横坐标取值就为负数。（当然无论实际设定在哪个象限，计算出来的结果都是一样的。）

为了更好的呈现这个落点的过程，我做成了动态可视化，就像下面这样。

可以从上面的动态模拟随机落点看到，使用计算机不停的模拟随机落点，由于圆的半径是1，所以小于等于1单位的随机点都标记为红色，其余则是绿色。最终红色部分很明显现形成了之前提到的$\frac{1}{4}$个单位圆部分。

代码部分如下：

import numpy as np  
import pandas as pd  
import math  
import random  
import matplotlib.pyplot as plt

random.seed(123)

# 总随机数量
total = 50000

# 区域内的数量
in_count = 0  
x, y = [], []

# 红点x，y坐标
x_red, y_red = [], []  
# 绿点x，y坐标
x_green, y_green = [], []


for i in range(total):  
    # random()方法默认为0，1边界。另外，由于演示的可视化为第二象限，所以x坐标取值为负数。
    x.append(-random.random())
    y.append(random.random())

    # 也可以直接计算(x**2 + y**2)**0.5，因为开根号就等于0.5次方
    distance = math.sqrt((x[i]**2 + y[i]**2))

    if distance <= 1:
        in_count += 1
        x_red.append(x[i])
        y_red.append(y[i])
    else:
        x_green.append(x[i])
        y_green.append(y[i])

fig, ax = plt.subplots(figsize=(6, 6))  
plt.scatter(x_red, y_red, c='red', edgecolor='black')  
plt.scatter(x_green, y_green, c='green', edgecolor='black')  

print(f'π: {4*in_count/total}')

上面的代码会输出最终的可视化图形，并且输出π值为: 3.13232，可见这个数值并非很精准，因为只设定了最大随机落点50000次，如果随着数量的增大就越会逼近最优解。

最后，要说的是蒙特卡洛方法并没有什么高深的理论支撑，如果一定要说有理论也就只有概率论或统计学中的大数定律了。蒙特卡洛的基本原理简单描述是先大量模拟，然后计算一个事件发生的次数，再通过这个发生次数除以总模拟次数，得到想要的结果。蒙特卡洛方法当然也可以运用在很多领域，如金融，工程，物理，生物医学等等。

描述是指对现有数据的总结和提炼，原始数据是杂乱无章的，所以将原始数据通过某种形式浓缩成一个有意义的统计量，比如通过图表形式对所收集的数据进行加工处理和显示；或将一系列复杂的数据序列减少为几个能够起到描述作用的数字（比如一套多难度复杂的体育动作浓缩为9.8分）。但是任何一种简化都会面临被滥用的危险。

在对任何数据做分析时有个前提，那必然是在拿到数据后，结合业务对数据做一个充分了解，我所拥有的数据是来自淘宝的购买商品 与 婴儿信息 两个数据集。

接下来要对淘宝婴儿的商品数据集做探索性数据分析（EDA），EDA对于数据分析是十分重要的环节，很多新手在拿到数据的第一刻就急冲冲的开始套用各种分析方法紧接着生成一个好看的统计图表，可最终到分析报告时也没有得出个值得的分析结果。而仅仅是对现有粗糙数据的一个“表现”统计而已。

既然是EDA那么对数据分析的过程中先尽量不考虑任何理论先验假设，而是做一些初步的数据探索，比如是否存在缺失值，异常值，查找数据结构和规律等，这样在探索的过程中会随着不断的深入对数据理解更加深刻，也方便后续进一步开始更深入数据分析。

这里我使用Excel工具来对数据进行EDA，这里使用的数据依然是前篇提到的淘宝用户购买婴儿用品的数据集，来源于：Baby Goods Info Data。

数据探索

数据列重命名

考虑到原始数据集都是英文字段也没有对应的信息，为了方便后续的分析展现和理解，这里做了一份副本复制，并对数据集里的数据列标题进行中文重命名。

检查是否有缺失值与重复值

重复值检查：从实际业务角度考虑，由于这是用户购买的数据记录，那么不同时间段相同用户存在相同的购物记录也是合情合理的，所以这里就不做重复值的检查了。

缺失值检查：我发现在商品属性里缺少了144条数值，均为空白值，但是考虑到拿到的数据集目前商品属性没有对应的信息说明，暂时也用不到我先将他隐藏起来。

统一日期格式

由于两个数据集里都包含有日期字段，而且类型为 常规 类型。这里使用excel中 【数据】-->【分列】-->【日期】 功能将原始常规类型转换为统一的日期格式。

多表关联

表2婴儿信息表中含有 生日日期 与 性别 字段，可以使用vlookup函数把这两个字段关联到表1购买商品信息表中，这样让数据集产生关联性就可以得到更多的数据信息。

由于性别里有些是空值或者属于未知性别（数值为2），所以将它们进行筛选剔除，对应的生日日期也显示正常了。

现在，在表1中关联合并中都包含了生日日期与购买时间，思考下就想到了，利用这两个字段的数据可以计算在用户购买商品时他们的婴儿年龄是多少，只需要用【购买时间】减去【生日日期】就可以得到了。所以我添加 【年龄】字段，使用DATEDIF(startdate,enddate,unit)函数来计算两个日期的时间，也就是（生日日期，购买时间，"M"），这里的M代表以月龄来计算，因为考虑到婴儿年龄一般都不大，可能会出现很多0岁的情况。但是可能会遇到一种情况就是用户在购买婴儿商品的时候，婴儿也存在还没出生的情况。所以当我使用如上时间计算月龄的时候就遇到了部分错误值，这里还需要使用IFERROR函数对错误值做一个替换，将错误值统一解释标记为”未出生“，使用方法是输入函数【=IFERROR(DATEDIF(H2,F2,"m"),"未出生")】。

分析数据结构和构建模型

基本的数据清洗都完成了，接下来可以使用数据透视表来进行一进步的统计分析观察。

哪个季度是销售旺季？

不同婴儿性别的家庭，用户购买的销量是否有显著差别？

在整理数据后，通过数据和数据之间的关系得到了更多新的信息，接下来也会从可视化角度去展现这些问题的答案。

淘宝婴儿商品数据之图表可视化

一、不同婴儿性别的家庭中，用户购买的二级商品最受欢迎的前10商品？

从下图可以发现，从上往下看商品 50018831 至 50011993 是最总计排名最高的最受欢迎十项商品，并且可以发现相比较之下拥有女性婴儿的家庭购买这些物品的比例最高。

二、购买这些商品的用户的婴儿性别占比

可以很直观的发现拥有女性婴儿家庭的用户占了购买商品的大部分，可能某些商品更受女性婴儿家庭欢迎。

三、四个季节的销售情况

可以非常直观的看到第四季度销售情况最好，该季自然也就是销售旺季。

四、一年中每个月的销售分布情况

从如下折线图中可以发现一年中明显销售火爆的月份是5月份与11月份。

五、TOP5购物最多的用户排行 这五位用户359601689、259538915、917524288、299196791、486110123是所有数据中购物最多的前5名。

数据分析之业务指标

对于互联网产品数据分析师来说，搭建指标体系可以很好的梳理业务关系，提高问题分析效率。本篇是对业务数据指标以及各指标能解决的相关问题做一个梳理和案例分析。6

常见的指标有哪些？分别有什么用？

在理解指标有哪些时，先要明白什么是指标，指标就是用某个统一标准去衡量业务，这个统一标准就是指标，通常是数值或者比率。

在不同的业务场景中有不同的数据指标和对指标的表达，但是经过高度抽象化，有几大类数据指标是所有类型应用共同分析的需要和参考。如下我用思维导图罗列几个常见的数据指标为：用户指标、行为指标、产品指标。

案例1：淘宝婴儿商品数据的 业务指标分类 以及 对应指标能解决的问题

这里使用上面提到的 淘宝婴儿商品数据集 进行业务指标分析，该数据包一共有2个数据表，我们将对该数据集进行指标分类和对不同数据指标解决的问题做一个说明。

一、淘宝婴儿商品数据 业务指标分类

按照之前罗列的几种常见业务指标对该数据集每列字段进行大分类，如下：

用户数据：用户id、出生日期、婴儿性别 行为数据：购买数量、购买时间 产品数据：商品id、商品属性、商品一级分类、产品二级分类

二、对应的业务指标能解决哪些问题

如上已经按照 淘宝婴儿商品 数据表中的字段大致分了3大类型数据指标，分别为：用户数据、行为数据、产品数据，这几大类数据指标中分别有不同的指标可以解决一些问题，如下：

• 用户数据：用户id、出生日期、婴儿性别

1.「新增用户」指标，可以知道网站每日新增长用户。
2.「留存率」指标，可以计算某时间段内的一批用户是否依然还保持访问网站（即使是每日签到打卡），以此观察网站对用户的粘贴性强弱。
3.「活跃用户率」指标，可以用于统计不同时间段内，访问网站的人数净值（去重复）占总用户数的比率，对于在网站的任何操作、停留时长等都可以统计为活跃用户范围，该比率比「留存率」更凸显网站的健康度和质量。
因为针对的婴儿商品，可以统计来自不同婴儿性别的用户比率比重，以此增加对应性别比重的产品，增大销量。

• 行为数据：购买数量、购买时间

1.通过用户购买商品的数量和购买时间排行可以观察出哪些时间段是销售旺季，以及用户青睐哪些商品。
产品数据：商品id、商品属性、商品一级分类、产品二级分类 「成交数量」指标，可以分析哪些一级或二级商品种类是被用户购买最多的。

案例2：喜马拉雅app是如何根据业务来选择指标，进行数据分析的？

对于互联网产品数据分析师来说，搭建指标体系可以很好的梳理业务关系，提高问题分析效率。确立关键指标（北极星指标），是统一各团队的努力方向。

如下是对喜马拉雅app的一个业务指标体系罗列。

生日悖论（Birthday paradox）是指假设一个班级有50个人，如果说在这个班级里概率大到可以肯定的说至少有2个人的生日相同（当然这里还不包括双胞胎，不包括闰年2月29日的情况），你信吗？

一般情况下，我们的直觉会认为班级里至少有两个人生日相同的概率会比较低，毕竟每个人的生日有365种选择，而班级只有50人，但是实际上计算得到在50个人的班级里出现同生日的概率甚至达到了惊人的97%！

正是因为理性的计算与日常的直觉经验产生了如此明显的矛盾，该问题才被称为生日悖论。

生日的概率计算

我们选择要用理性的方式来计算生日相同的概率是多少，但是对于计算有多少人生日相同，这个理解起来可能有点复杂和一时摸不着头脑。其实可以逆向思考下，反过来计算下50个人的班级里每个人生日不在同一天的概率是多少，再用总概率1减去生日不同的概率，最后不就能得到该班级找那个生日相同的人的概率是多少了吗。

先随便抽了1位同学，他的生日肯定可以是365天的任意一天，那么用数字计算表示就是：$\frac{365}{365}$，因为他有365中可能性。

注意：或许，很多人会认为应该是$\frac{1}{365}$，理由是每个人的生日只占一年中的一天。这里要理清一点，我们不是求每位同学日期的占比率，如果是占比率那任何一个人的生日都只可能是$\frac{1}{365}$。我们需要计算的是第一位同学的生日在一年中有多少种可选性，这是在你不知道的情况下计算的，所以第一位同学的生日可以有365种可能，（事件可能数量 / 总事件数）所以就是$\frac{365}{365}$。

接下来再抽第2位同学，因为要计算的是每个人生日不同的概率，那么前一位同学无论他的生日日期最终是哪一天，肯定会占用365天中的一个日期名额。所以第2位同学的生日可选择的数量只剩下364种了，就是$\frac{364}{365}$。而第3位同学就是$\frac{363}{365}$，以此类推到第50位同学时可选择的范围前面已占用49个名额，剩下就是$\frac{365-49}{365}=\frac{316}{365}$。

由于每个人都是随机独立事件，所以我们将该班级每个人的生日可能概率想乘，得到的即是全班50个人生日各自不同情况下的概率。

$\hat{p}$为假设班级里50位同学各自生日都不同的概率：$\hat{p}=\frac{365}{365}\times\frac{364}{365}\times\frac{363}{365}...\times\frac{316}{365}$

$p$为我们最终想知道的至少两人生日相同的概率，既然前面已经求得所有人生日不同的概率$\hat{p}$，那么拿总数1减一下就可以得出至少有两人生日相同的概率：$p=1-\hat{p}=1-\frac{365}{365}\times\frac{364}{365}\times\frac{363}{365}...\times\frac{316}{365}$

现在，计算包含 n 位同学时至少有两人生日相同的概率$p$为多少，转换成通用公式如下：

$\hat{p}=\frac{365\times364\times363...\times(365-n+1)}{365^{n}}=∏\limits_{i=1}^{n}\frac{365-i+1}{365}$

$p=1-\hat{p} = 1-∏\limits_{i=1}^{n}\frac{365-i+1}{365}$

同生日的概率分布可视化

总结出通用公式后就可以计算在出基于不同人数时至少两人生日相同的概率。

我用Python可视化了不同人数时的概率分布情况，你可以很直观的看到下图中红色横向的虚线是50%概率的位置，当x轴人数达到23左右的时对应的y轴概率基本已经在50%，也就是有50%的情况会出现至少有两个人生日相同，而当人数在60左右的时候蓝色曲线几乎水平并可以说接近100%了。

实现代码如下：

import math  
import matplotlib.pyplot as plt 

def birthday_paradox(size):  
    if size > 0:
        x_size = [i+1 for i in range(size)]
        y_prob = []
        log_x = 0

        for i in range(size):
            # 考虑到乘积数字太大，使用log的方式求解。
            # log(a)N = x 表示以a=10为底，真数N为365/365..365-i/365求值，该得的值为指数最后被使用。
            log_N = (365-i)/365
            log_x += math.log(log_N,10)
            p_hat = 10 ** log_x
            p = 1 - p_hat
            y_prob.append(p)

    plt.figure(figsize=(10,5))
    plt.plot(x_size, y_prob, linewidth=2.5, label='prob(n)', color='blue')

    # 设置水平线
    plt.axhline(0.5, linestyle='--', color='red', label='50% prob')
    # 绘制图例
    plt.legend()
    # 设置x轴范围
    plt.xlim(0, size)
    # 设置网格线型
    plt.grid(linestyle='-.', alpha=0.5)
    # 设置窗口标题
    plt.title('Visualization of Birthday Paradox', fontsize=16);
    # 设置坐标轴标签
    plt.xlabel('Peolple', fontsize=12)
    plt.ylabel('Probability', fontsize=12)


birthday_paradox(365)  

模式验证

通过之前对生日悖论的概率计算及可视化，我们明确知道了在不同人数群体时获得至少两个人生日相同的概率分布情况。但是如果实际真的出现指定数量的人数时，是否真的能接近我们计算的概率呢？

我做了一个模拟实验，实验的构思大致为先生成一年365天真实的日期范围（不用在意年份），测试人数为1至100人。因为我们要观察不同人数群体生日相同的概率，所以也就会出现100个群体样本。如，第1个样本只有1人，第2个样本有2人，以此类推，第100个样本就包含100个人。而每个样本的人群都会随机分配生日日期，接着针对这100个样本再分别做1次，10次，100次，1000次的实验，观察生日是否出现相同的实验。直到每个群体的反复实验完成后，用 相同生日的累积次数 / 总的循环实验次数 就可以获得每个特定人群出现生日相同的概率值。理论上来说，每个人群数量得出的实验概率值会与生日悖论的概率相同，那么两条概率曲线应当是吻合的。

如下是可视化的结果，可以发现只做1次实验的时候，1至50人时的紫色概率曲线波动十分明显，结果不是0.0就是1.0（100%）。回忆起之前在生日悖论概率计算时，当群体达到20人时存在两人生日相同的概率应当在0.41左右，但我们的实验结果却是1.0。这其实很好解释，正因为概率低所以实验的次数越少实验结果才越不稳定，只执行1次实验的话结果可想而知，要么存在要么不存在，类似于一个阶跃函数。要知道我们得到的生日悖论概率其实是一个期望概率，只有经过多次实验也就是根据大数定律才会趋向于最终概率值。可以看到当我们运行的次数越来越多直到1000次的时候，紫色的实验概率曲线最终与蓝色曲线基本吻合在一起，如果继续增加实验次数肯定会重叠，也就说明模式验证的概率是符合的。注：另外要说明的是，模拟测试的生日概率是平均分布的（现实生活中，出生机率不是平均分布的）。

代码部分如下：

import pandas as pd  
import numpy as np  
import datetime  
import time  
from collections import Counter

# 验证生日悖论概率
x = []  
y = []  
date_range = pd.date_range('2019/1/1', periods=365).strftime("%Y-%m-%d").to_list()  
def birth_paradox_test(size, times):  
    global x
    global y
    x = [s+1 for s in range(size)]

    for i in range(1, size+1):
        same_date_sum = 0
        for _ in range(times):
            value = Counter(np.random.choice(date_range, size=i, replace=True)).most_common(1)[0][1]
            if value > 1:
                same_date_sum += 1
        prob = same_date_sum/times
        y.append(prob)


# 生日悖论概率分布        
x_size = []  
y_prob = []  
def birthday_paradox(size):  
    global x_size
    global y_prob
    if size > 0:
        x_size = [i+1 for i in range(size)]
        y_prob = []
        log_x = 0

        for i in range(size):
            # 考虑到乘积数字太大，使用log的方式求解。
            # log(a)N = x 表示以a=10为底，真数N为365/365..365-i/365求值，该得的值为指数最后被使用。
            log_N = (365-i)/365
            log_x += math.log(log_N,10)
            p_hat = 10 ** log_x
            p = 1 - p_hat
            y_prob.append(p)



plt.figure(figsize=(20,12))  
# 循环4次，每次以1, 10, 100, 1000循环次数。
times_list = [1, 10, 100, 1000]  
map_no = 1  
# 设置人数最大为100人
size = 100

for time in times_list:  
    birthday_paradox(size)
    birth_paradox_test(size, times=time)

    plt.subplot(2,2, map_no)
    map_no += 1

    plt.plot(x_size, y_prob, linewidth=2.5, label='prob(n)', color='blue')
    plt.plot(x, y, linewidth=2.5, label='prob_test(n)', color='purple')

    # 设置水平线
    plt.axhline(0.5, linestyle='--', color='red', label='50% prob')
    plt.legend()
    plt.xlim(0, size)
    plt.grid(linestyle='-.', alpha=0.5)
    plt.title(f'Visualization of Birthday Paradox Test in {time} times', fontsize=16);
    plt.xlabel('Peolple', fontsize=12)
    plt.ylabel('Probability', fontsize=12)

    x = [] 
    y = []
    x_size = [] 
    y_prob = []

扩展思考

既然我们已经通过理性的计算和了解了一个班级存在至少两个人生日相同的概率如此之高，但是你有没有发觉，回想自己的学生生涯，在自己经历各个年级阶段时，却很少发现自己与班级里的同学生日相同，这是为什么呢？似乎又产生矛盾了。这其实也可以通过理性计算来解惑。

要计算与我（或者说班级里特定某个人）的生日相同的概率是多少，我们依然使用之前的逆向思维来思考，可以计算与我生日不同的人的概率有多少，最后相减一下就可以得出有多少概率可以使生日相同。

假设还是在50个人的班级里，首先我的生日（也可以说特定某个人的生日）可以是一年中的任意一天（实际情况肯定知道自己的生日是哪一天，但是从理性计算角度而言可以允许是一年里的任一天），那么我自己的生日可以有365种可能性，概率就是$\frac{365}{365}$。别忘了，现在计算的是生日不相同的概率，所以第二个人的生日不能与我的生日相同，那么他的生日只能有364种可选性，概率就是$\frac{364}{365}$。重点来了，第三个人他的生日依然有364种可能，为什么呢？因为生日的日期参照值是我的生日，只要与我生日的日期不同，哪怕第三个人和别人生日存在相同都没关系，所以依然还有364种可选性，那么概率依然是$\frac{364}{365}$，以此类推，直到班级里第50个人，他与我的生日不同的概率依然是$\frac{364}{365}$。

所以在50个人的班级里与我生日不同的概率计算如下：

$\hat{p}=(\frac{364}{365})^{49} = \frac{364}{365}\times\frac{364}{365}...\times\frac{364}{365}$

我们转换成通用公式如下：

$\hat{p} = (\frac{364}{365})^{n-1}$

最后用总数1减去$\hat{p}$就可以得到和我的生日相同的概率：$p = 1-\hat{p}=1-(\frac{364}{365})^{n-1}$。

插播一个新问题：如果问题的问法是，在这50个人的班级里是否有人的生日是8月8日（注意，这里并不是问是否和班级里特定某个人生日一样，而是问生日日期是否和指定的一个日期相同）？

如果这样问，那么生日不是指定日期的概率$\hat{p}=(\frac{364}{365})^{50}$，应该是50次方而不是49次方，因为我们计算的$\hat{p}$是逆向计算生日不同于8月8日这个日子的人数，而这个班级中完全有可能50人的生日都不是8月8日。

回归正题，老样子我们同样用可视化描绘出不同人数时与自己生日相同的概率分布情况。

我们可以发现，大概到250人的时候，才有50%的几率会有至少1个人与我的生日相同。所以如果我们仅仅关注的是一群人里有至少两个人的生日相同，那只需要满足有23个人时就可以达到50%，而如果以自己或者某一个人的生日为参照中心的时候这个人数的要求就很高了。

实现代码如下：

import math  
import matplotlib.pyplot as plt 

def birthday_paradox(size):  
    if size > 0:
        x_size = [i+1 for i in range(size)]
        y_prob = []
        y_prob_2 = []
        log_x = 0

        # 便于理解，统一序号以1开始，range(1, n)只显示1,2,3..n-1。所以最终
        for i in range(size):
            # 考虑到乘积数字太大，使用log的方式求解。
            # log(a)N = x 表示以a=10为底，真数N为365/365..365-i/365求值，该得的值为指数最后被使用。
            # 注：因为这里i是从0开始，所以直接365-i即可，不需要365-i+1。
            log_N = (365-i)/365
            log_x += math.log(log_N,10)
            p_hat = 10 ** log_x
            p = 1 - p_hat
            y_prob.append(p)

            # 注：需要循环累积乘以n-1次，因为这里i是从0开始，所以不需要i-1。
            p_hat_2 = math.pow(364/365, i)
            p_2 = 1 - p_hat_2
            y_prob_2.append(p_2) 

    plt.figure(figsize=(10,5))
    plt.plot(x_size, y_prob, linewidth=2.5, label='prob(n)', color='blue')
    plt.plot(x_size, y_prob_2, linewidth=2.5, label='prob2(n)', color='green')

    # 设置水平线
    plt.axhline(0.5, linestyle='--', color='red', label='50% prob')
    # 绘制图例
    plt.legend()
    # 设置x轴范围
    plt.xlim(0, size)
    # 设置网格线型
    plt.grid(linestyle='-.', alpha=0.5)
    # 设置窗口标题
    plt.title('Visualization of Birthday Paradox', fontsize=16);
    # 设置坐标轴标签
    plt.xlabel('Peolple', fontsize=12)
    plt.ylabel('Probability', fontsize=12)


birthday_paradox(365)  

所以，综合上述来说，我们的潜在直觉并没有错，错的是我们没有从正确的角度去理解问题。因此，当我们拨开直觉的谎言去理解问题，才会觉得如此不可以思议。

整个数据文件一共有两个分别为data_ncov.xlsx和chinadata.json。前一个文件是我们的疫情数据集文件，后一个则是后续创建地图时用到的全国各地省市的地理数据信息。

# COVID-19 数据分析

import pandas as pd  
import numpy as np  
import matplotlib.pyplot as plt  
import warnings

# 设置不弹出警告
warnings.filterwarnings('ignore')

# 中文乱码设置
plt.rcParams['font.sans-serif']=['SimHei']  
plt.rcParams['axes.unicode_minus']=False

# retian屏幕显示设置
%config InlineBackend.figure_format = 'retina'

# 读取数据
url = '/Users/Anders/Documents/Jupyter/DataSet/COVID-19/data_ncov.xlsx'  
# 数据集文件格式为xlsx，故使用read_excel方法读取文件
df = pd.read_excel(url)  
df.head()  

先读取数据集文件并输出前5条数据内容格式如下：

在进行基本的数据清洗前，可以通过info方法先了解数据的基本信息，如字段类型、数据一致性、是否存在空值等。

df.info()  

输出如下：

<class 'pandas.core.frame.DataFrame'>  
RangeIndex: 782 entries, 0 to 781  
Data columns (total 7 columns):  
 #   Column      Non-Null Count  Dtype 
---  ------      --------------  ----- 
 0   Unnamed: 0  782 non-null    int64 
 1   区域编码        782 non-null    int64 
 2   省市          782 non-null    object
 3   疑似          782 non-null    int64 
 4   确诊          782 non-null    int64 
 5   死亡          782 non-null    int64 
 6   date        782 non-null    int64 
dtypes: int64(6), object(1)  
memory usage: 42.9+ KB  

从如上输出可以发现date字段为int64类型，并不是日期类型，考虑到后续需要与日期字段打交道，为了更方便比较，这里选择先将该字段转化为字符串类型更合适。

# 当前数据日期信息的类型为整型，所以先将字段转化为字符串
df['date'] = df['date'].astype('str')  
df['date'] = pd.to_datetime(df['date'])  

一、各省市疫情排行

为了查看各省市的疫情数据排行，编写一个函数可以方便查看（除了湖北省以外）10大省市排行最高的确诊类型数据，该函数参数支持自定义，如指定时间、类型、几大排行。

def fig_top_type(time, ty, top_num):  
    # 根据传入的时间参数'time'筛选指定时间日期。
    df_data = df[df['date']==time]

    # 根据传入的指定类型'ty'排序，按照从大到小顺序排序，并修改原数据。
    df_data.sort_values(by=ty, ascending=False, inplace=True)

    # 因为考虑到当前数据湖北省远远领先，只取湖北省外的数据所以切片从1开始
    df_data.iloc[1:top_num+1].plot(x = '省市', y=ty, 
                                kind='bar', figsize=(15,5), 
                                grid=True, rot=45,
                                title='湖北省外{0}日，{1}病例最多的{2}省市排行'.format(time, ty, top_num))

# 输出一个除了湖北省外2020年2月1日关于确诊和疑似的柱状图
fig_top_type('20200201',['确诊','疑似'], 10);  

这里很清楚的可以看到2020年2月1日排名前三的分别是浙江省、广东省、河南省。

二、(全国) 每日病例与增长率数据展示

之前我们以省市为单位展示了数据排行的情况，但是如果想知道全国每天的新增数据情况呢？

由于目前提供的数据都是每个省市每一天的累积更新，想获得每天新增的各类病例，就需要计算得到。

既然以每天的日期为单位，就可以使用groupby针对日期进行分组，把所有城市每天的疫情数据汇总起来，按照 "疑似"、"确诊"、"死亡"三类进行分类分组并汇总求和。

# 以日期分组，获取全国每天的 疑似、确诊、死亡 数据。
data_china = df.groupby('date')[['疑似','确诊','死亡']].sum()  
data_china.head()  

输出如下：

因为通过观察得知每天的原始数据是累积的数据，所以要得到每天各类的新增数据就必须拿后一天的数据减去前一天数据，那么使用shift方法可以使数据进行前移或者后移，其中shift(1)代表向下平移1位，shift(-1)代表向上平移一位。同时我们添加【新增**】列字段保存每日新增的数据。

# 为了计算后一日对前一日的数据差计算，使用shift方法可以用于同行单纯的前移或后移操作
data_china['新增疑似'] = data_china['疑似'] - data_china['疑似'].shift(1)  
data_china['新增确诊'] = data_china['确诊'] - data_china['确诊'].shift(1)  
data_china['新增死亡'] = data_china['死亡'] - data_china['死亡'].shift(1)  
data_china.head()  

有了各项新增数据后，也可以添加更直观的增长率（百分比）。比如确诊的增长率就等于用当天的新增确诊数据 除以 前一日的确诊数据就代表当天的确诊增长率。由于直接使用新增确诊除以确诊在这里除以的是相同行内当天的确诊数据，所以我们必须使用shift使数据向下移动一位。

# 计算新增确诊的增长率就等于用新增的数据除以新增前的数据。
data_china['确诊_增长率'] = data_china['新增确诊'] / data_china['确诊'].shift(1)  
data_china['确诊_增长率'] = data_china['确诊_增长率'].round(3)  
data_china['疑似_增长率'] = data_china['新增疑似'] / data_china['疑似'].shift(1)  
data_china['疑似_增长率'] = data_china['疑似_增长率'].round(3)  
data_china.head()  

基本的数据清理完毕后，我们打算开始进行可视化，这里选择使用百度的Echarts做图形可视化。

# https://pyecharts.org/#/zh-cn/intro
# Echarts 是一个由百度开源的数据可视化，而pyecharts是封装了Echarts的python库，
# 为了使用python能更方便的使用Echarts图形。注：新版本V1以后版本要使用新的调用方法。
from pyecharts.charts import Bar, Grid, Line  
from pyecharts import options as opts  
from pyecharts.globals import ThemeType

# x轴为日期为单位分类。考虑到第一行数据本身不存在新增对比，数据存在Nan值，所以我们舍弃第1行，当然如果需要也可以做合理的填充处理。
x = data_china.iloc[1:].index.astype('str')  
y1 = data_china[['确诊','疑似']].iloc[1:]  
bar = (  
    Bar(init_opts=opts.InitOpts(theme=ThemeType.INFOGRAPHIC, width="700px", height="350px"))
    .add_xaxis(list(x))
    .add_yaxis('确诊',list(y1['确诊']))
    .add_yaxis('疑似',list(y1['疑似']), gap='10%')
    .extend_axis(yaxis=opts.AxisOpts())
    .set_series_opts(label_opts=opts.LabelOpts(is_show=False))

    .set_global_opts(title_opts=opts.TitleOpts(title="全国累计确诊病例 柱状图", title_textstyle_opts=opts.TextStyleOpts(font_size=15)),
                     datazoom_opts=opts.DataZoomOpts(),
                     tooltip_opts = opts.TooltipOpts(axis_pointer_type='cross'))
)

y2 = data_china['确诊_增长率'].iloc[1:]  
y3 = data_china['疑似_增长率'].iloc[1:]  
line = (  
    Line()
    .add_xaxis(x)
    .add_yaxis('确诊增长率', y2, yaxis_index=1, is_smooth=True)
    .add_yaxis('疑似增长率', y3, yaxis_index=1, is_smooth=True)
    .set_series_opts(label_opts=opts.LabelOpts(is_show=False))
)

bar.overlap(line)  
bar.render_notebook()  

使用Echarts可视化后输出的图形可以很方便的进行数据的交互。

三、中国疫情地图

用地图展现疫情数据可以更直观的从地域上来了解疫情的严重程度的分布情况。当然地图的可视化方法有很多，这里我们分别使用GeoPandas库和Echarts包内自带的Map库制作。

• 使用GeoPandas库实现可视化

# GeoPandas是一个开源项目，它的目的是使得在Python下更方便的处理地理空间数据。
# GeoPandas扩展了pandas的数据类型，允许其在几何类型上进行空间操作。

import geopandas as gpd

geo_url = '/Users/Anders/Documents/Jupyter/DataSet/COVID-19/chinadata.json'

# 通过geopandas.GeoDataFrame.from_file()方法来读取数据，对于空间数据的格式常见的为*.json与.shapefile(GIS类软件常用)
china_spatial = gpd.GeoDataFrame.from_file(geo_url)  
china_spatial.info()  

读取chinadata.json文件后使用info方法查看信息输出如下。

<class 'geopandas.geodataframe.GeoDataFrame'>  
RangeIndex: 34 entries, 0 to 33  
Data columns (total 4 columns):  
 #   Column     Non-Null Count  Dtype   
---  ------     --------------  -----   
 0   name       34 non-null     object  
 1   centerlng  34 non-null     float64 
 2   centerlat  34 non-null     float64 
 3   geometry   34 non-null     geometry
dtypes: float64(2), geometry(1), object(1)  
memory usage: 1.2+ KB  

可以看到geometry字段数据类型为geopandas.geodataframe.GeoDataFrame，并不是pandas读取后的pandas.DataFrame，但由于geopandas结合了pandas和shapely的功能，所以对于表格的处理方法基本和pandas一致。

可以这样理解，GeoDataFrame是向DataFrame增加了地理数据支持的功能。我们也能发现新的数据字段类型：geometry，该类型字段代表的是空间信息，比如这里的数据其实是我国省级行政区划的面数据，面数据的空间属性为：面中每个节点的经纬度坐标构成。

china_spatial.head()  

通过head方法查看提取的chinadata.json前几行内容如下。

为了让地图里每个省市有自己的对应数据，我们就需要把疫情数据和对应的省市匹配起来。假设我们需要用地图的方式显示2020年2月1日的疫情数据，那么我们

data_0201 = df[df['date'] == '2020-2-1']  
data_0201.head()  

接下来开始绘制地图，逻辑分为：先把地理数据与疫情数据合并，接着创建画布，然后使用合并后带有地理坐标的疫情数据集去结合plot方法绘制，最后是借助plt.text文本方法给地图上的每个地理位置标注省市的名称。

# 把地图数据与疫情数据匹配起来
data_china_0201 = pd.merge(china_spatial, data_0201, left_on = 'name', right_on = '省市', how = 'left')  
data_china_0201.drop(['name'], axis=1, inplace=True)  
data_china_0201.head()

# 创建画布
plt.figure(figsize=(12,9))  
plt.title('2020-2-1 全国确诊病例', fontsize = 20)

# 绘制疫情地图，这里我们只显示确诊信息，所以设置参数column='确诊'
data_china_0201.plot(ax=plt.subplot(1,1,1), alpha=1, edgecolor='k', linewidth = 0.5,  
                  legend=True, scheme = 'FisherJenks', column='确诊', cmap = 'Reds')
# 设置网格线
plt.grid(True,alpha=0.5)

# 添加省市信息
lst = data_china_0201[['省市','centerlng','centerlat','确诊']].to_dict(orient = 'record')

for i in lst:  
    plt.text(i['centerlng'], i['centerlat'], i['省市'] +':' + str(i['确诊']))

现在就可以很方便的从地图上看到哪个省市到2月1日为止的疫情最为严重，颜色深度越深代表疫情越严重。

除了如上以不同颜色来展示外，还可以在地图中绘制以气泡图的方式来展现。绘图逻辑和上一关大体相同，不过地理空间图将会作为底图，也就不需要指定column参数。这里需要用到matplotlib的scatter方法将xy轴参数去匹配地理坐标的经纬度数据，而其中s参数就是疫情的确诊数据。

# 把地图数据与疫情数据匹配起来
data_china_0201 = pd.merge(china_spatial, data_0201, left_on = 'name', right_on = '省市', how = 'left')  
data_china_0201.drop(['name'], axis=1, inplace=True)

plt.figure(figsize=(14,18))  
plt.title('2020-2-1 全国确诊病例', fontsize = 23)

# 绘制底图地图时不需要显示一列数据，所以不需要指定column参数，只需要显示一个底图颜色。
data_china_0201.plot(ax=plt.subplot(1,1,1),  
                  edgecolor='k', linewidth = 0.5, color = 'gray', alpha = 0.1)

# 添加气泡图
plt.scatter(data_china_0201['centerlng'], data_china_0201['centerlat'],  
            s = data_china_0201['确诊'], edgecolors='k', alpha = 0.8)

# 设置网格线
plt.grid(True,alpha=0.5)

# 添加省市信息
lst = data_china_0201[['省市','centerlng','centerlat','确诊']].to_dict(orient = 'record')  
for i in lst:  
    plt.text(i['centerlng'], i['centerlat'], i['省市'] +':' + str(i['确诊']))

如下图，就可以很直观的以气泡的方式发现哪个省市的疫情气泡大就说明疫情更严重。

甚至还可以包装成一个函数，批量创建气泡地图，考虑到批量创建，地图显示局促，所以我们取消了地理名称的文本显示。

def create_map(time, tp, x, y, n):  
    # 按照日期筛选数据
    datai = df[df['date'] == time] 
    # 匹配数据
    data_chinai = pd.merge(china_spatial, datai, left_on = 'name', right_on = '省市', how = 'left')
    del data_chinai['name']  # 删除多余字段
    # 绘制底图
    data_chinai.plot(ax=plt.subplot(x,y,n),
                      edgecolor='k', linewidth = 0.5,
                      color = 'gray', alpha = 0.1)
    # 添加气泡图
    plt.scatter(data_chinai['centerlng'],data_chinai['centerlat'], 
                s = data_chinai[tp], edgecolors='k', alpha = 0.8)
    # 设置标题及网格线
    plt.title('%s 全国%s病例' % (time, tp), fontsize = 20)
    plt.grid(True,alpha=0.5)


# 构建for循环批量出图
# 设置日期列表
datelst = ['2020-1-22','2020-1-30','2020-2-1','2020-2-5']

# 创建绘图对象
plt.figure(figsize=(20,18))

# 批量出图
m = 1  
for i in datelst:  
    create_map(i,'确诊', 3,2,m)
    m += 1

下图批量显示了'2020-1-22','2020-1-30','2020-2-1','2020-2-5'4个时间的疫情气泡大小变化。

• 使用pyecharts库实现可视化

pyecharts库中负责地理坐标系的模块是Geo，负责地图的模块是Map，负责百度地图的模块是BMap。

其中Geo实现了一个地理坐标系，地图上的点可以与经纬度进行转换(即可以利用经纬度向地图中插入点，也可以获取地图上某一点的经纬度)，实现地图上的打点功能主要依靠Geo类来进行，而Map功能类似于Geo，但只有地图，没有坐标系，即地图上的点无法与经纬度进行转换。而负责图表配置的模块是options。在 pyecharts 中，图表的一切皆通过 options来修饰调整。

from pyecharts.faker import Faker  
from pyecharts import options as opts  
from pyecharts.charts import Map

# 整理原先数据集里的地理名称，为了使与pyecharts内的map包内的地理名称所匹配
dict_special_str = {'广西壮族自治区':'广西', '内蒙古自治区':'内蒙古',  
 '宁夏回族自治区':'宁夏', '西藏自治区':'西藏', '新疆维吾尔族自治区':'新疆'}
special_str_func1 = lambda x: x.replace('市','').replace('省','').replace('特别行政区','')  
provinces = data_china_0201['省市'].map(special_str_func1)

def special_str_func2(s):  
    if dict_special_str.setdefault(s):
        return dict_special_str[s]
    else:
        return s
provinces = list(map(special_str_func2, provinces))


value = list(data_china_0201['确诊'])


def map_base():  
    c = (
        Map()
        .add("确诊", [list(z) for z in zip(provinces, value)], "china")
        .set_global_opts(title_opts=opts.TitleOpts(title="中国2月1日疫情地图"),
                         visualmap_opts=opts.VisualMapOpts(is_piecewise=True, max_=max(value), min_=0))
    )
    return c


city_map = map_base()  
city_map.render_notebook()  

网上的动态数据排行在我了解之后主要发现是用javascript写出来的，但是基于对javascript没有那么深入了解，我找到了其它可替代方案，就是使用matplotlib的animation方法来绘制动图。

效果展示

最终的效果还不错，可以看下GIF效果动图。

设计逻辑

对于整个代码的设计逻辑主要就是两方面，数据 和 动态展示。在实际编写前首先需要思考下如何获取数据和如何使数据动起来达到自己的设想。

对于获取数据而言，本身打算自己去用爬虫获取数据，但是鉴于各个信息源的数据都比较杂乱，我搜索发现已经有人做了相关的数据汇总并开放使用，可以说已经有了半成品的疫情数据集，那就可以直接拿来用，只需要对数据做基于自己场景下的数据清洗就可以了。

数据问题已经解决，接下来就是考虑如何达到最终目的使数据动态化。按照文章开头提到的，只需要使用matplotlib的animation方法就可以绘制动图。

代码分解

开始的部分主要是做数据清洗的部分，尽管我们获取的公开数据集已经很好的汇总了国内以及全球的疫情每日信息，但是由于制作动图的所需数据格式不一样，所以还是需要做一定的数据清洗和整理。

# 完整代码
import os  
import pandas as pd  
import matplotlib as mpl  
import matplotlib.pyplot as plt  
import matplotlib.ticker as ticker  
import matplotlib.animation as animation  
from IPython.display import HTML

# 数据源地址
url='https://raw.githubusercontent.com/canghailan/Wuhan-2019-nCoV/master/Wuhan-2019-nCoV.csv'  
df_csv = pd.read_csv(url, parse_dates=['date'], low_memory=False)  
df_csv.sample(10)  

完成导入数据后，随机抽取10行数据查看内容如下。

由于数据里既包含了国家又包含了国内各省市的数据，而我们只需要罗列各国的数据，那就需要做数据分类处理。经过观察，可以看到上面的内容里country字段包含了各个国家的信息，而只要是以国家为统计单位的话province这一列的字段是Nan值，所以为了方便处理这里做数据填充，统一把Nan空值改为0。

# 由于原始数据包含的外国数据是不存在省份信息，所以统一设置为0
df_csv['province'].fillna(0, inplace=True)

# 只要省份信息为0，说明它属于一个国家。同时country字段筛选出'中国'
china_data = df_csv[(df_csv['country']=='中国') & (df_csv['province']==0)]  
non_china_data = df_csv[df_csv['country']!='中国']

为了验证数据是否正确，使用head()方法查看前5行数据的内容如下。

china_data.head()  

non_china_data.head()  

我们发现上面chinadata和nonchina_data两个变量中储存的中国和非中国的疫情数据集完全符合预期的格式结果。接下来需要去掉一些不相关的字段列，并且将这两大数据集合并起来。

# 为了后续的数据叠加，只需要各个国家以及确诊，死亡等相关信息，其它删除。
china_data.drop(['countryCode', 'cityCode', 'province', 'provinceCode', 'city'], axis=1, inplace=True)  
non_china_data.drop(['countryCode', 'cityCode', 'province', 'provinceCode', 'city'], axis=1, inplace=True)

# 将中国和其它外国数据叠加
all_countries_data = china_data.append(non_china_data)  
all_countries_data.reset_index(drop=True, inplace=True)  

数据叠加后，养成好习惯再次查验下内容是否是我们想要的。

all_countries_data.sample(10)  

格式精简过的数据看起来十分干净并且很正确。

接着，开始进入可视化环节。为了更好的区分各个国家，所以用不同颜色来表示，世界上的国家不少，这里只罗列出疫情数据相对排名突出的一些国家。

# 疫情中主要出现的一些国家的相应颜色
colors = dict(zip(  
    ['美国', '德国', '伊朗', '瑞士', '荷兰', '韩国','比利时', '奥地利', '土耳其',
     '葡萄牙', '加拿大', '挪威', '澳大利亚', '巴西', '以色列', '瑞典',
     '捷克', '丹麦', '马来西亚', '爱尔兰', '日本','巴基斯坦', '俄罗斯', '泰国', 
     '沙特阿拉伯', '南非', '芬兰', '印度尼西亚', '菲律宾', '希腊',
     '冰岛', '印度', '新加坡', '巴拿马', '墨西哥', '意大利', '钻石公主号邮轮', '阿根廷',
     '法国','塞尔维亚', '伊拉克', '巴林','新西兰', '黎巴嫩', '阿尔及利亚', '阿联酋', 
     '乌克兰', '拉脱维亚', '西班牙', '越南', '中国', '英国'],
    ['#FFB7DD', '#FF88C2','#FF44AA','#FF0088','#C10066','#A20055','#8C0044',
    '#FFCCCC','#FF8888','#FF3333','#FF0000','#CC0000','#AA0000','#880000',
    '#FFC8B4','#FFA488','#FF7744','#FF5511','#E63F00','#C63300','#A42D00',
    '#FFDDAA','#FFBB66','#FFAA33','#FF8800','#EE7700','#CC6600','#BB5500',
    '#FFEE99','#FFDD55','#FFCC22','#FFBB00','#DDAA00','#AA7700','#886600',
    '#FFFFBB','#FFFF77','#FFFF33','#FFFF00','#EEEE00','#BBBB00','#888800']))

如下代码主要创建了一个构建图例的函数，这个函数我们会在之后通过matplotlib下的animation函数不停的循环传入日期参数以此达到动态效果。

fig, ax = plt.subplots(figsize=(15, 8))

# 定义画横向柱状图函数
def draw_barchart(current_date):  
# 只筛选排名前10的国家
    df = all_countries_data[all_countries_data['date']==pd.to_datetime(current_date)].sort_values(by='confirmed', ascending=False).head(10)

  # 因为横向显示，所以数据多的国家考前排名 
    df = df[::-1]

    ax.clear()

    # 开始画柱状图，并根据之前颜色字典分配指定国家颜色，由于使用了setdefault方法那么字典中不存在就使用默认颜色。
    ax.barh(df['country'], df['confirmed'], color=[colors.setdefault(x,'#90d595') for x in df['country']])

    dx = df['confirmed'].max() / 200
    for i, (country, confirmed) in enumerate(zip(df['country'], df['confirmed'])):
        ax.text(x=confirmed-dx+300, y=i+0.02, s=country, size=12, weight=600, ha='right', va='center')  # Tokyo: 名字
        ax.text(x=confirmed+dx, y=i, s=f'{confirmed:,.0f}', size=16, ha='left', va='center')

    # 对样式的修改
    ax.text(0.8, 0.15, current_date, transform=ax.transAxes, color='#777777', size=30, weight=800)
    ax.text(0, 1.05, '确认感染病例', transform=ax.transAxes, size=12, color='#777777')
    ax.text(0.97, 0.1, 'code by Anders', transform=ax.transAxes, size=12, color='#777777', ha='right')
    ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:,.0f}'))
    ax.xaxis.set_ticks_position('top')

    # 将x轴确诊病例的最大范围设置为所有国家最高的确诊病例的1.1倍
    plt.xlim(0, all_countries_data['confirmed'].max()*1.1)
    ax.tick_params(axis='x', colors='#777777', labelsize=12)
    ax.set_yticks([])
    ax.margins(0, 0.01)
    ax.grid(which='major', axis='x', linestyle='-')
    ax.set_axisbelow(True)
    ax.text(0, 1.1, '世界COVID-19确诊病例排行榜',
            transform=ax.transAxes, size=24, weight=600, ha='left')
    plt.box(False)

接下来我们要用到animation.FuncAnimation，而函数FuncAnimation(fig,func,frames,init_func,interval,blit)是绘制动图的主要函数，其参数如下：

• fig 绘制动图的画布名称
• func自定义动画函数，即下边程序定义的函数update
• frames动画长度，一次循环包含的帧数，在函数运行时，其值会传递给函数update(n)的形参“n”。
• init_func自定义开始帧，即传入刚定义的函数init,初始化函数
• interval更新频率，以ms计。
• blit选择更新所有点，还是仅更新产生变化的点。应选择True，但mac用户请选择False，否则无法显。

最后我们只需要设置一个时间列表，并传入时间列表给FuncAnimation方法，即可以html的形式显示我们的动态可视化数据排行。

# 调用图例表生成函数
date_time = ['2020-03-25', '2020-03-26', '2020-03-27', '2020-03-28',  
             '2020-03-29', '2020-03-30', '2020-03-31', '2020-04-01']


animator = animation.FuncAnimation(fig, draw_barchart, frames=date_time)  
HTML(animator.to_jshtml())  

完整代码

# 完整代码
import os  
import pandas as pd  
import matplotlib as mpl  
import matplotlib.pyplot as plt  
import matplotlib.ticker as ticker  
import matplotlib.animation as animation  
from IPython.display import HTML

# 数据源地址
url='https://raw.githubusercontent.com/canghailan/Wuhan-2019-nCoV/master/Wuhan-2019-nCoV.csv'  
df_csv = pd.read_csv(url, parse_dates=['date'], low_memory=False)

# 由于原始数据包含的外国数据是不存在省份信息，所以统一设置为0
df_csv['province'].fillna(0, inplace=True)

# 只要省份信息为0，说明它属于一个国家。同时country字段筛选出'中国'
china_data = df_csv[(df_csv['country']=='中国') & (df_csv['province']==0)]  
non_china_data = df_csv[df_csv['country']!='中国']

# 为了后续的数据叠加，只需要各个国家以及确诊，死亡等相关信息，其它删除。
china_data.drop(['countryCode', 'cityCode', 'province', 'provinceCode', 'city'], axis=1, inplace=True)  
non_china_data.drop(['countryCode', 'cityCode', 'province', 'provinceCode', 'city'], axis=1, inplace=True)

# 将中国和其它外国数据叠加
all_countries_data = china_data.append(non_china_data)  
all_countries_data.reset_index(drop=True, inplace=True)

# 疫情中主要出现的一些国家的相应颜色
colors = dict(zip(  
    ['美国', '德国', '伊朗', '瑞士', '荷兰', '韩国','比利时', '奥地利', '土耳其',
     '葡萄牙', '加拿大', '挪威', '澳大利亚', '巴西', '以色列', '瑞典',
     '捷克', '丹麦', '马来西亚', '爱尔兰', '日本','巴基斯坦', '俄罗斯', '泰国', 
     '沙特阿拉伯', '南非', '芬兰', '印度尼西亚', '菲律宾', '希腊',
     '冰岛', '印度', '新加坡', '巴拿马', '墨西哥', '意大利', '钻石公主号邮轮', '阿根廷',
     '法国','塞尔维亚', '伊拉克', '巴林','新西兰', '黎巴嫩', '阿尔及利亚', '阿联酋', 
     '乌克兰', '拉脱维亚', '西班牙', '越南', '中国', '英国'],
    ['#FFB7DD', '#FF88C2','#FF44AA','#FF0088','#C10066','#A20055','#8C0044',
    '#FFCCCC','#FF8888','#FF3333','#FF0000','#CC0000','#AA0000','#880000',
    '#FFC8B4','#FFA488','#FF7744','#FF5511','#E63F00','#C63300','#A42D00',
    '#FFDDAA','#FFBB66','#FFAA33','#FF8800','#EE7700','#CC6600','#BB5500',
    '#FFEE99','#FFDD55','#FFCC22','#FFBB00','#DDAA00','#AA7700','#886600',
    '#FFFFBB','#FFFF77','#FFFF33','#FFFF00','#EEEE00','#BBBB00','#888800',
    '#EE9A00','#EE9572','#EE82EE','#EE8262','#EE7AE9','#EE799F','#EE7942',
      '#EE7621','#EE7600','#EE6AA7']))

fig, ax = plt.subplots(figsize=(15, 8))

# 定义画横向柱状图函数
def draw_barchart(current_date):

    df = all_countries_data[all_countries_data['date']==pd.to_datetime(current_date)].sort_values(by='confirmed', ascending=False).head(10)
    df = df[::-1]

    ax.clear()

    # 开始画柱状图，并根据之前颜色字典分配指定国家颜色，由于使用了setdefault方法那么字典中不存在就使用默认颜色。
    ax.barh(df['country'], df['confirmed'], color=[colors.setdefault(x,'#90d595') for x in df['country']])

    dx = df['confirmed'].max() / 200
    for i, (country, confirmed) in enumerate(zip(df['country'], df['confirmed'])):
        # 设置国家名字
        ax.text(x=confirmed-dx+300, y=i+0.02, s=country, size=12, weight=600, ha='right', va='center') 
        # 设置确诊数值
        ax.text(x=confirmed+dx, y=i, s=f'{confirmed:,.0f}', size=16, ha='left', va='center')

    # 对样式的修改
    ax.text(0.8, 0.15, current_date, transform=ax.transAxes, color='#777777', size=30, weight=800)
    ax.text(0, 1.05, '确认感染病例', transform=ax.transAxes, size=12, color='#777777')
    ax.text(0.97, 0.1, 'code by Anders', transform=ax.transAxes, size=12, color='#777777', ha='right')
    ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:,.0f}'))
    ax.xaxis.set_ticks_position('top')

    # 将x轴确诊病例的最大范围设置为所有国家最高的确诊病例的1.1倍
    plt.xlim(0, all_countries_data['confirmed'].max()*1.1)
    ax.tick_params(axis='x', colors='#777777', labelsize=12)
    ax.set_yticks([])
    ax.margins(0, 0.01)
    ax.grid(which='major', axis='x', linestyle='-')
    ax.set_axisbelow(True)
    ax.text(0, 1.1, '世界COVID-19确诊病例排行榜',
            transform=ax.transAxes, size=24, weight=600, ha='left')
#     ax.text(1, 0, ',by QIML',, transform=ax.transAxes, ha=',right',,
#             color=','#777777',, bbox=dict(facecolor=',white',, alpha=0.8, edgecolor=',white',))
    plt.box(False)


date_time = [  '2020-03-25', '2020-03-26', '2020-03-27', '2020-03-28',  
               '2020-03-29', '2020-03-30', '2020-03-31', '2020-04-01']


animator = animation.FuncAnimation(fig, draw_barchart, frames=date_time)  
# 在jupyter上显示出动态效果
HTML(animator.to_jshtml())  

最后可以使用FFmpeg视频解码器，将我们的动态效果一帧帧转为指定的mp4视频文件。

# 输出mp4动画，由于我们使用了ffmpeg来解码，所以需要指定系统中ffmpeg执行文件的位置
ffmpegpath = os.path.abspath("/Users/Anders/Documents/Jupyter/ffmpeg")  
plt.rcParams["animation.ffmpeg_path"] = ffmpegpath  
writer = animation.FFMpegWriter(extra_args=['-vcodec', 'libx264'])  
animator.save('covid19_animation.mp4', writer=writer, dpi=180)  

I. 问题的定义

• 项目概述
• 问题陈述
• 评价指标

II. 分析

• 数据的探索
• 探索性可视化
• 算法和技术
• 基准和模型

III. 方法

• 数据预处理
• 执行过程
• 完善

IV. 结果

• 模型的评价与验证
• 合理性分析

V. 项目结论

• 结果可视化
• 对项目的思考
• 需要作出的改进

I. 问题的定义

项目概述

无论是飞速发展的互联网电商行业还是传统的零售行业，销售预测在企业的整个运营体系中都是必不可少的环节。所谓销售预测，是在对影响市场供求变化的众多因素上进行系统地调查和研究，并在此基础上运用科学的方法对未来市场产品的供需发展趋势以及有关的各种因素的变化，进行分析、预见、估计和判断。传统的销售预测方法往往只考虑了一部分影响销售的因素，其建立的模型也相对简单。而数据挖掘技术是一种科学有效的数据处理方式，它为应对信息爆炸，海量信息的处理提供了科学有效的手段。计算机数据挖掘技术顺应了时代和社会的发展，也逐渐成为社会关注的焦点。为了优化企业商品销售决策方案，提高商品销售预测的准确率，应用数据挖掘方法对销售数据库进行分析，提高销售预测的准确率无疑是十分有意义的。麻省理工学院(MIT)和宾夕法尼亚大学(University of Pennsylvania)的经济学家进行的一项研究发现，在数据驱动的决策规模上，一个标准差高的公司生产率高5%，利润高6%，市值高50%。

本项目是Kaggle¹上的一个数据科学竞赛项目，项目的目的是通过由一家名为Rossmann的公司所提供的往年销售数据，运用机器学习的手段进行数据挖掘并构建出一个高效和尽可能精准的数据预测模型，通过这个模型以后只需要输入相关特征信息就可以预测未来一段时间内的销售额。相对于传统低效而又不精准的人工销售预测方式，对于Rossmann这样一家拥有3500家门店、47500名员工、2015年销售收入79亿欧元的公司来说，通过这个数据模型来预测未来一定时间段内的销售额，相信这对公司的市场战略把握和对门店方面的管理来说是非常有价值的。²

问题陈述

Rossmann是德国最大的日化用品超市，它在欧洲7个国家拥有超过3500家的连锁商店。当前，为了更加有效地管理采购和库存，Rossmann商店的门店经理有一项重要的任务就是预测未来6周时间内每家分店的每日销售额。而商店的销售额会受到很多因素的影响，其中包括促销活动、竞争对手、节假日和一些季节性和局部的因素。

我们需要通过Rossmann公司所提供的1115家门店的历史销售数据来训练和测试一个最佳销售预测模型，这个模型最终可以用来预测未来一段时间内接近真实销售额的数据。

预测销售额这在机器学习中属于一个回归性问题，其本质就是构建一个最佳映射函数 $\hat{y}=f(x)$，既给定任意新的输入变量来预测输出变量值。在这里，x 为输入变量，变量经过函数 $f$ 后预测输出最终的销售额 $\hat{y}$。预测建模主要关注的是如何最小化模型的误差，这里的误差则通过真实的销售额 $y$ 值与预测的销售额 $\hat{y}$ 之间的差值来评估。

评价指标

模型所采纳的评估指标为Kaggle在竞赛中所推荐的 Root Mean Square Percentage Error (RMSPE) 指标。

$RMSPE = \sqrt{\frac{1}{n}\sum\limits_{i=1}^n\left(\frac{y_i-\hat{y}_i}{y_i}\right)^2} = \sqrt{\frac{1}{n}\sum\limits_{i=1}^n\left(\frac{\hat{y}_i}{{y}_i}-1\right)^2}$

其中 $y_i$ 代表门店当天的真实销售额，$\hat{y}_i$ 代表相对应的预测销售额，$n$ 代表样本的数量。如果有任何一天的销售额为0，那么将会被忽略。最后计算得到的这个RMSPE值越小代表误差就越小，相应就会获得更高的评分。

II. 分析

数据的探索

Rossmann所提供的数据可以到Kaggle网站上下载³，下载的数据包中一共包含四个文件，分别是train.csv, test.csv, store.csv, sample_submission.csv，这几个文件的相关描述如下。

train.csv 数据集文件共有 1,017,209 个数据样本（由 9列 x 1017209行 数据组成），是用于训练模型的数据集。其中包含了1115家商店从 2013年1月1日 至 2015年7月31日 范围内每家商店每天销售的相关信息，例如有当天是否营业、当天的顾客数量、是否参加促销、当天的销售额，等信息。其中有180家商店缺少从2014年7月1日至 2014年12月31日约半年的数据。

要说明一点是，其中存在的Sales字段就是我们用于训练的预测数据标签，其余的字段信息将归类于特征数据，但最终需要经过一定特征工程处理后用来作为模型的训练输入数据。

store.csv 数据集文件共有 1115 个数据样本（由 10列 x 1115行 数据组成），其中包含了对应1115家门店的额外补充信息，例如有商店等级、商店类型、与附近的竞争对手商店，等相关信息。

提示：由于原始字段的长度不易于排版展示，所以对某些字段里的部分单词进行了缩写展示，如：CompetitionOpenSince缩写为COS、Promo2Since缩写为P2S。

test.csv 数据集文件共有 41,088 个数据样本（由 8列 x 41088行 数据组成），其中包含了1115家商店从 2015年8月1日 至 2015年9月17日 范围内每家商店每天销售的相关信息，包含的信息与 train.csv 数据集类似，仅仅区别是其中剔除了Sales销售额字段，因为这是一个独立的测试数据集，用来测试我们的模型最终得分。由于真实的销售额并没有在数据集中公开，所以在进行模型最终测试后预测出来的结果（销售额）最终需要上传到Kaggle网站上才能看到实际得分排名。

samplesubmission.csv 作用是提供一个上传到Kaggle测试评分的标准格式文件，也就是基于之前提到的预测模型与test.csv 测试集文件所测试出来的预测数据，并把这个预测数据以按照 samplesubmission.csv 的格式上传到Kaggle上。

考虑到数据预处理在进行机器学习训练数据前是非常重要的一个步骤，经过技术手段初步筛查、统计后发现：在train.csv数据集中一共有172,871个样本销售额为0，但其中有54个销售额为0的样本当天属于营业状态。在store.csv数据集中字段名为 CompetitionDistance 缺失值有 3 个，CompetitionOpenSinceMonth和CompetitionOpenSinceYear 缺失值 354 个，Promo2SinceWeek和Promo2SinceYear和PromoInterval 缺失值 544 个。在test.csv数据集中Open字段中存在 11 个缺失值。在后面数据预处理环节，我们会对这部分缺失的数据和异常值的进行预处理。

在所有这些数据集中，大部分字段的含义都是一目了然的。另外，仍然有一些需要说明的数据字段，包括如下：

探索性可视化

在项目的早期阶段，我们已经大致弄清楚如何构建数据，确定如何操纵可用的数据以获得所需的答案，比如包括获取对数据的直觉、比较变量的分布、对数据进行检查、发现数据中的缺失值和异常值等。

而通过可视化去探索数据也是数据分析工作中的重要组成部分。在庞大的数据面前，通过图形的方式去呈现数据分析的结论，可以帮助我们更加直观的理解数据。

• 所有数据样本的销售额分布情况

如上，图1 左小图展示的是在训练集中，所有数据样本的销售额的分布情况。从图中可以观测到，数据主要是左偏态分布的情况，销售额为0的样本量非常高，接近竖轴刻度标记的17500数值。除此以外，从局部看剩下的大部分数据形成了一个小正态分布，从蓝色柱状图比例看大致集中的范围在0到15000。为了便于观测，我将销售额集中的区域进行了选取放大。如上图中，右小图是基于左小图进行了缩放后的结果，可以发现销售额在5000左右的样本量频率是最高的，也就是说大约有125000个样本量数据它们所产生的消费额在5000左右。

为了进一步探究为什么之前显示销售额为0的样本量非常高，我将所有销售为0的数据筛选出来，进行了查询分类并做了可视化，如下图。

从 图2 左小图中可以发现，可以说几乎占比100%的商店由于没有营业属于休业状态才导致了销售为0。由于在之前的数据探索部分，我已经通过技术手段筛查出所有销售额为0的数据，并从数据中发现了在超过17万个样本中仅仅只有54个样本显示销售额为0但是商店是属于营业状态的。初步认为，这种极小的情况是合理的，因为不管当天有没有顾客访问商店，当天没有产生实际销售额的情况也是有可能的。同时，为了可视化能更清晰的查看，右小图通过缩放范围来展示。

• 促销对销售额的影响

图3 使用小提琴图形的方式展示了每周几中某一天，进行促销与不促销对销售额产生的影响。可以很清楚的发现进行了促销活动后的商店销售额（橙色）明显高于没有进行促销的商店。我们还可以发现一个特别的现象，在周六和周日，商店均没有开展过促销活动。且周日的销售额非常低，这也是符合逻辑的，因为在德国，绝大多数情况下 (除了节假日前后)，绝大部分商店都是不营业的，因此总的销售额会远低于工作日以及周六。

另外，从图上可以发现小提琴图的顶端非常高，说明样本存在很高的异常值区间，我们会在后面数据处理环节进行处理。

• 销售与竞争对手商店距离的影响

考虑到store.csv数据集是对商店相关信息的额外补充，为了探索潜在的数据信息，所以我将store.csv数据集和train.csv数据集进行了合并，并对 商店的销售额 与 附近竞争对手商店的距离（也就是数据集里的CompetitionDistance字段）进行了一个可视化输出，从 图4 展示可以发现并没有因为竞争对手离Rossmann商店的距离越远，Rossmann商店的消费额就变多，所以没有产生什么影响（其中要说明的是，由于总样本数据超过百万数量，为了方便展示所以我抽取了1%的数据点在不影响分布的情况下进行展示。以及原点的颜色越深代表该样本量的销售额越高，原点的大小越大代表该样本的顾客量越大）。相反离竞争对手越近的一些商店它的销售额都很高。

• 销售额与顾客数量之间的关系

图5 第一二两张子图分别展示了 不同商店的销售额分布 与 不同商店的顾客数量分布。从两张图的比对结果可以发现，当第一张子图中每家商店对应的销售额增高时，第二张子图中相同区域的顾客数量也出现一定峰值，所以可以考虑给为数据集增加一个 人均消费 (PerCustomerSale字段) 。第三张子图展现的是不同商店的人均消费分布，我们可以通过观察大致发现在人均消费分布上并没有之前两张子图显示的波动明显，相对适中。

最后从第一二张子图与第三张子图对应的人均消费来看，峰值区域没有形成相互对应，说明在人均消费额差距不大的情况下，几家商店销售额的突出是以跑量的方式换来的。

• 销售额与商店类型之间的关系

• 假期对销售额的影响

图7 显示了国定假日与学校放假时对销售额的影响情况，可以观察发现，国定假日期间的销售额有一定的提升说明假日对销售额是有一定影响的。而学校放假与不放假对销售额几乎没有多少影响。

• 2013~2015年中每月的销售额分布情况

因为数据集中包含了不同时间年份的数据，为了根据不同日期维度来统计信息，我将日期（Date字段）进行了分割处理，其中由于训练集的数据只包含时间跨度为 2013年1月1日 至 2015年7月31日 的范围，所以2015年8月份以后的数据并没有显示。

算法和技术

最初，考虑的是决策树回归作为基准模型，但是观察看到许多商店的销售模式明显偏离了平均行为，一个单一的模型将很难处理这些商店提供的所有特殊情况，以及考虑到在Kaggle上得分比较高的模型算法都使用了XGBoost集成算法。所以我认为集成大量模型的方法最适合在当前的情况下进行预测，在这里将XGBoost作为目标模型。

决策树回归使用的是 CART 算法，既可用于分类也可用于回归，如果待预测结果是连续型数据，则 CART 生成回归决 策树。区别于 ID3 和 C4.5，CART 假设决策树是二叉树，回归树选取 Gainσ为评 价分裂属性的指标。选择具有最小 Gainσ的属性及其属性值作为最优分裂属性 以及最优分裂属性值。Gain_σ值越小，说明二分之后的子样本的“差异性”越小，说明选择该值作为分裂值的效果越好⁷。

XGBoost是集成学习中 Boosting 流派中的其中一种算法，其本质上还是一个 GBDT，XGBoost与GDBT的区别主要有¹²：

1. XGBoost生成CART树考虑了树的复杂度，GDBT未考虑，GDBT在树的剪枝步骤中考虑了树的复杂度。

2. XGBoost是拟合上一轮损失函数的二阶泰勒展开，GDBT是拟合上一轮损失函数的一阶泰勒展开，因此，XGBoost的准确性更高，且满足相同的训练效果，需要的迭代次数更少。

3. XGBoost与GDBT都是逐次迭代来提高模型性能，但是XGBoost在选取最佳切分点时可以开启多线程进行，大大提高了运行速度。

XGBoost核心算法原理基本为：

（1）不断地添加树，不断地进行特征分裂来生长一棵树，每次添加一个树，其实是学习一个新函数，去拟合上次预测的残差。

（2）当我们训练完成得到k棵树，我们要预测一个样本的分数，其实就是根据这个样本的特征，在每棵树中会落到对应的一个叶子节点，每个叶子节点就对应一个分数。

（3）最后只需要将每棵树对应的分数加起来就是该样本的预测值。⁸。

基准模型

预测销售额是一个定量的数值，所以使用回归模型而不是分类模型。由于这个项目与波士顿房价（Boston Housing）项目比较类似都是属于回归预测的问题，而且波士顿房价项目作为入门项目，其中所用到的模型简单易于理解，所以这里决定选用决策树作为基准模型。可以通过引入 sklearn 库中的 DecisionTreeRegressor模块来创建一个回归决策树模型，然后用网格搜索训练法指定不同的参数组合，把所有可能的组合生成网格，然后用交叉验证进行评估，最后遍历寻找出最优的模型参数，使用这个最优参数得到的 RMSPE 得分即基准模型的分数。

III. 方法

数据预处理

在训练一个模型之前需要做数据的预处理，数据预处理的目的是为了保证数据的质量，以便能够更好的为后续的分析、建模工作服务，因为模型的最终效果决定于数据的质量和数据中蕴含的有用信息的数量。通常在拿到数据以后，首先要判断此数据是否可为我们所用，也就是我们根据需求目标所拿到的数据的质量是否过关。³

在之前的数据的探索环节，已经提到过在正式训练模型前有一些前提准备工作要做。如，解决数据的不一致，用合理的值替换缺失值，删除不必要的数据，将分类变量转换为数值，检查异常值并纠正等。下面我将对这些处理的步骤进行描述。

• 数据集合并

在之前的数据的探索环节已经提及到为了更好的探索数据，考虑到store.csv数据集是对商店相关信息的额外补充，所以已经使用pandas模块库将store.csv数据集和train.csv数据集进行了成功合并。

在合并时考虑到store.csv数据集包含了1115家商铺信息正好与train.csv里的1115家商店编号标签吻合，所以我使用pandas.DataFrame.merge函数并设置on参数来合并，on参数值就是用于指定连接两个数据表的相同列名，这里设置为(on="Store")。

• 数据不一致与缺失值处理

在之前的（图2）探索性可视化部分，从数据中探索出了超过17万个数据样本因为门店是休业状态所以销售额为0，但是其中依然有45家门店为营业状态而销售额也是0。据Kaggle推测，这在某些情况下是由于商店进行了改造，或进行了试营业，实际上并没有向公众开放。所以我将这45家销售额为0的门店的Open状态从营业更改为休业状态（Open = 1改为Open = 0），这样它们就不会对预测模型构建产生负面影响。⁴

接着是对缺失值的处理，通过对缺失值筛查后发现在store.csv数据集中字段名为 CompetitionDistance 缺失值有 3 个，CompetitionOpenSinceMonth和CompetitionOpenSinceYear 缺失值 354 个，Promo2SinceWeek和Promo2SinceYear和PromoInterval 缺失值 544 个。在test.csv数据集中Open字段中存在 11 个缺失值。处理缺失值的方式有很多，必须要以分析环境和模型的需求来考虑。通常在数据量很多的情况下可以考虑删除法来删除整行数据，也可以考虑数据填补补齐的方式使数据表更加完备。因为已经提前将store.csv与train.csv两个数据表进行了整合，再次查看显示的缺失值信息与单独查看store数据集的空值会多很多，在这里我选择填补的方式进行缺失值处理。我分别对CompetitionOpenSinceMonth、CompetitionOpenSinceYear、Promo2SinceWeek、Promo2SinceYear进行了中位数的填充。

• 检测并处理异常值

异常值会大幅度地改变数据分析和统计建模的结果。检测异常值的方法有许多，这里我使用了箱线图的四分位距（IQR）对异常值进行检测⁵，也叫Tukey's Test，下图（图9）就是箱线图去除异常值的原理。

箱线图中间是一个箱体，也就是粉红色部分，箱体左边，中间，右边分别有一条线，左边是下四分位数（Q1），右边是上四分位数（Q3），中间是中位数（Median），上下四分位数之差是四分位距（称IQR），IQR也可以说指的是两个四分值之间的范围大小。用 Q1-1.5*IQR 得到下边缘（最小值），Q3+1.5*IQR 得到上边缘（最大值），在上边缘之外的数据就是极大异常值，在下边缘之外的数据极小异常值，总之在上下边缘之外的数据就是异常值。另外要说明的是图中的数值 1.5 是一个可变的系数，表示的是中度异常；对于重度异常的情况，可以视情况而定将系数提高至 3。

我在处理异常值时预先使用了常规的中度异常系数1.5去实现箱线图最小值与最大值的划分。经过异常值处理后保留了990,515个样本数据，正常数据保留占比约97.37%，去除了约2.63%的异常数据，所以对整体的数据质量不会产生多少影响。

• 时间格式字段分离处理

在数据集中Date字段默认显示的时间格式为年-月-日，这个格式的数据是不能在机器学习中直接被处理的，所以需要进行一个时间格式的字段分离处理，这里直接单独分离为年（Year）、月（Month）、日（Day）三个新字段变量。同时去除原始的Date日期字段变量数据。

• 数据类型转换

对于一些特殊字段它们的值是a、b、c、d，为了使模型能够顺利对它们进行处理，分别修改为对应的1、2、3、4。

执行过程

在本项目中，如之前的算法和技术环节提到，我使用了两个模型进行销售额预测分别为DecisionTreeRegressor作为基准模型以及XGboost作为目标模型。

• 数据探索

在训练模型前将所有原始数据进行数据探索与可视化探索，从而可以进一步发现各类数据变量之间的关联和潜在信息。

• 数据清洗与特征工程

这个项目存在多个数据集文件，所以在进一步数据探索后进行数据集文件的合并，由于数据集的合并产生新的特征信息，结合特供工程化处理对特征变量进行了构建补充与删除。

• 定义评价函数

基于本项目推荐使用的RMSPE评价指标，所以定义了RMSPE函数来作为模型的预测结果评分。

• 构建数据模型并预测

在构建模型前先需要划分训练模型和验证模型的必备数据，而测试数据集test.csv已经额外提供给我们，所以可以将train.csv中的数据全部都用来训练，换句话说，只需要将train.csv中的数据分为训练集和验证集。考虑到这个项目属于一个有关时间序列的预测类问题，因为这是和时间序列有关的数据集，所以我将train.csv训练集数据进行了细分，按照时间顺序划分验证集和训练集，其中取最后6周时间作为验证集，剩下的为训练集。

接着是对模型的参数定义，DecisionTreeRegressor的参数定义是我鉴于在曾经用机器学习预测波士顿房价（Boston Housing）项目时的经验，使用了网格化搜索指定不同的参数组合的方式来选择最优参数。

对于XGBoost我之前并不了解，XGBoost模型有许多参数可以设置，我参考了Kaggle和其它资料进行了参数的指定⁶。因为在构建模型过程中定义了RMSPE评估函数，所以将参考最优的train eval rmspe的得分来决定最终模型。为了得到最低的RMSPE得分，我对参数进行了一定调整其中调整最多的是eta参数，这个参数可以理解为是学习速率，我将eta从0.03调整为0.01后发现迭代的次数也逐步增加，在eta为0.03时2900次就终止了迭代，当设置为0.01后迭代次数增加，直到6000次才停止，花费了更多的训练迭代次数，也味着要花很长时间才会收敛。最后我根据最优的RMSPE得分设定eta为0.02训练了出最终模型。最后，用这个模型进行对项目提供的test.csv测试集进行预测，预测的值同时会乘以权重因子以得到更精确的预测值（权重因子参考Kaggle¹⁰），然后提交的到Kaggle进行最终评测。

要说明的一点是，关于模型的调参其实可以选择GridSearch方式来选择一系列最佳参数，在本项目中考虑到当前的模型得分已经满足条件，所以受限于时间和运算能力的情况下我并没有最终使用GridSearch来进行参数优化。

完善

在XGBoost模型测试的过程中，我在参考了网上以及Kaggle上一些资料后，曾尝试调整增加了一些看似有用的特征，但发现在上传到Kaggle后最终模型的预测得分却反而非常不理想。然后在参数选择方面，一开始是使用默认的参数，发现训练结果一直存在不少的过拟合现象，本地的RMSPE得分非常不错，但是上传到Kaggle上最终预测时十分不理想。然后通过查阅资料⁹在XGBoost中调整参数（Complete Guide to Parameter Tuning in XGBoost with codes in Python），在参数调整后尽管得分有了一定的上升，但是我发觉经过多次的迭代后上升的幅度很小。于是我开始尝试在特征组合选取方面进行改动，发现特征的改动极大的使得分出现了改变。 所以，如果要进一步提高预测的精度，我认为未来应该需要重点在特征工程和 Xgboost 的调参上下更多功夫。

IV. 结果

模型的评价与验证

DecisionTreeRegressor模型使用网格化搜索对参数maxdepth, minsamplesleaf参数进行了搜索选择了，maxdepth = 20，minsamplesleaf = 6 为最优参数，使用这些参数建立的模型的最终训练集RMSPE得分是 0.21064，Kaggle最终得为0.63305。XGBoost的本地测试集得分为0.084679，Kaggle最终得分为0.11170，存在一定过拟合现象，但XGBoost依然优于DecisionTreeRegressor。

合理性分析

因为这是 Kaggle 上的一个竞赛，而且项目中提供的test.csv文件是一个独立的测试数据集，在一开始的test.csv作用描述中提到过它是用来测试模型的最终得分。由于真实的销售额并没有在数据集中公开，所以在进行模型最终测试后预测出来的结果（销售额）最终需要上传到Kaggle网站上才能看到实际得分排名。

对于DecisionTreeRegressor模型我多次测试下来本地的R2分数在0.2左右，通过上传到Kaggle后最好得分一直徘徊在0.63上下。而使用XGBoost的算法模型一开始并没有得到理想的得分，通过参考Kaggle上的资料¹⁰以及对特征变量的持续优化后，我得到了符合预期的RMSPE得分并在Kaggle上的Private Score得分最好的为0.11170，显然XGBoost作为目标模型是优于DecisionTreeRegressor。同时也已经满足项目要求的达到Top 10%。

V. 项目结论

结果可视化

图10展现了预测模型构建好之后的特征重要性，从图中看出 有关于时间的Day、PromoOpen 与 Store 这几特征的分数远高于其他特征。在构建模型的过程中我尝试调整了许多特征，但这几个特征的重要性一直都比较高，可以从这些特征中了解到销售额的高低与不同商店的类型和时间以及附近竞争对手开业的时间有极大的关系。

对项目的思考

尽管之前做过类似几个价格预测的项目，但是更多的是带有引导性去完成。通过这个数据挖掘项目，从头到尾拿到原始数据去构建整个分析和建模的过程，使我对整个数据挖掘的过程有了更清晰的认识。如：了解行业领域背景、数据探索、数据预处理、数据特征工程、模型的选择、模型的评价、最终验证。

基于时间问题，很多探索性测试没有完全去尝试，可以考虑把原来删除的异常值，用反预测的方式，去预测那些应该被删除的异常值数据他们原始的正常值应该是多少，这样在保持数据更齐全的情况下再去构建预测模型可能效果更好。

另外，从Kaggle上提交的分数排名看大家的分数之间差距并不是太大，也就是说大家在方法使用上，或者数据理解上，基本上差异性很小。而且基本上都使用了最好的XGBoost模型，可能在具体的尝试过程中，由于数据的准备不同才带来了一些细微差异。我们也并没有看到出现那种差距极大的队伍出现，也许提供的数据本身就代表了这个项目得分的局限性。因为在这个项目中可能存在更多潜在未知的信息可以影响模型的预测精准性，如当天天气变化情况，当地国情等。

需要作出的改进

在这个项目中依然还有特别多需要改进和完善的地方，比如数据探索的可视化部分，如果先从行业背景角度去理解然后再分析数据的关联性，结合之前的行业背景信息逐步推理相信可以进一步挖掘潜在的信息。而且这也为后续的特征工程构建和删除特征时可以提供更进一步的合理支撑。

在模型方面也可以尝试比其他优秀的Boosting类模型，如LightGBM，LightGBM是个快速的、分布式的、高性能的基于决策树算法的梯度提升框架，在速度上比XGBoost更快。¹¹

在特征选取方面，我认为特征的选择对模型的影响也是巨大的，在测试模型的过程中，我尝试了多种特征组合得到了非常大的分数反差，所以特征选择这块依然有很大的改进空间。考虑到Kaggle在讨论区里也有提供额外相关数据，例如各个药店所在的州统计的天气和降雨量数据等，通过增加更多的特征维度提升模型的精准度。

对于算法参数的改进，XGBoost是在数据分析中比较流行的集成算法，对于这个算法的很多参数我还没有很好的去理解，所以由于我的参数选择可能不是很合理导致模型最终得分并不十分理想，调参方面可以考虑使用GridSearchCV或者RandomizedSearchCV来进行合理参数的选择。

参考文献

检测到了异常值，我们需要对其进行一定的处理。而一般异常值的处理方法可大致分为以下几种：

• 直接将含有异常值的记录删除。

• 视为缺失值：将异常值视为缺失值，利用缺失值处理的方法进行处理。

• 平均值修正：可用前后两个观测值的平均值修正该异常值。

• 不处理：直接在具有异常值的数据集上进行数据挖掘。

具体如何处理异常值我们不在这里涉及，而是主要讲如何检测异常值，检测的分析方法多种多样，一般异常值的检测方法有基于统计的方法，基于聚类的方法，以及一些专门检测异常值的方法等。异常值会大幅度地改变数据分析和统计建模的结果。

由于检测分析的方法有许多，这里主要使用两种检测异常值的方法分别为：3σ准则（又称为 拉依达准则）与

箱型图 两种方式对异常值进行检测分析。

一、3σ准则

3σ准则 又称为 拉依达准则。所谓3σ，当数据被定义为在一组测定值中与平均值的 σ（标准偏差）超过3倍时，这个值就被认为是异常值，而这个异常值的概率通常小于0.3%，用公式可以理解为 $p(|x-\mu|>3 \sigma) \leqslant 0.003$ 。

也可以从如下的正态分布图中更好的理解（其中，$\mu$代表均值，$x=\mu$即为图像的对称轴），可以发现大部分数值都在 $(\mu-3σ，\mu+3σ)$ 范围内概率占到了99.73%，而超出这个范围区间的可能性仅占不到0.3%，这个区间范围内的数值就是异常值。

另外要注意的是，这种判别处理原理及方法仅局限于对正态或近似正态分布的样本数据处理，以及建议当数据量大的时候使用。

3σ准则 可视化部分

接下来，看看如何用可视化的方式查看异常值。之前提到过使用3σ准则有个前提就是数据必须服从正太分布，所以为了便于实验，这里使用numpy.random.randn()函数随机产生基于正太分布的数据值。

第一张图使用了KDE画出密度曲线后，根据3倍标准偏差画出对应红色虚线位置，这个红色虚线以内的范围占99.73%，而以外的范围占不到0.3%（为异常值范围）。

第二张图使用了散点图的方式将所有数据值以数据点的方式画出来，其中蓝色点为正常值，4个红色点就是异常值。而判断异常值的方式就是我们之前提到的公式：$|x-\mu|>3 \sigma$

总结3σ准则的分析大致为如下几个步骤：

1. 首先需要保证需要检验的数据列大致上服从正态分布。
2. 然后计算需要检验的数据列的标准差。
3. 最后比较数据列的每个值，是否大于标准差的3倍。
4. 选择用可视化的方式呈现出来。

3σ准则 代码实现

如下是代码的实现部分：

import numpy as np  
import matplotlib.pyplot as plt  
import pandas as pd  
from scipy import stats  
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  
plt.rcParams['axes.unicode_minus'] = False

# 创建虚拟数据
data = pd.Series(np.random.randn(1000))  
data.head()


# 用K-S检测正态性
p = stats.kstest(data, 'norm', (mu, std))  
if p.pvalue > 0.05:  
    print(f'p值为 {p.pvalue}是正态分布')
else:  
    print(f'p值为 {p.pvalue}不是正态分布')

# 构建画布
fig = plt.figure(figsize=(10, 8))  
ax1 = fig.add_subplot(211)

### 绘制kde密度图 ###
data.plot(kind='kde', grid=True, style='-k', title='密度曲线')

std = data.std()  
plt.axvline(3*std, color='red', linestyle='--')  
plt.axvline(-3*std, color='red', linestyle='--')



### 绘制散点图 ###
mu = data.mean()  
# 筛选出异常值|x-u|>3σ
outlier = data[np.abs(data - mu) > 3*std]

# 筛选出正常值
data_correct = data[np.abs(data - mu) <= 3*std]

ax2 = fig.add_subplot(212)  
ax2.scatter(data_correct.index, data_correct.values, edgecolor = 'black', alpha = 0.5)  
ax2.scatter(outlier.index, outlier.values, color = 'red', edgecolor = 'black')  

二、箱型图

有些数据不一定符合正太分布，如果用几倍的σ去检测异常值就不合适，那么我们如何来判断这些变化是否在合理的范围呢？可以考虑使用箱型图，箱型图的四分位距（IQR）对异常值进行检测，也叫Tukey's Test。我们可以看下如下箱型图的样式：

箱型图中间是一个箱体，也就是粉红色部分，箱体左边，中间，右边分别有一条线，左边是下四分位数（Q1），右边是上四分位数（Q3），中间是中位数（Median），上下四分位数之差是四分位距（称IQR）。IQR也可以说指的是两个四分值之间的范围大小。用 Q1-1.5*IQR 得到下边缘（最小值），Q3+1.5*IQR 得到上边缘（最大值），在上边缘之外的数据就是极大异常值，在下边缘之外的数据极小异常值，总之在上下边缘之外的数据就是异常值。另外要说明的是图中的数值 1.5 是一个可变的系数，表示的是中度异常；对于重度异常的情况，可以视情况而定将系数提高至 3。

箱型图 可视化

如下是基于我用python模拟的数据并以（横向）箱型图的方式呈现的可视化结果，从下图中可以很清楚的发现，左侧也就是下边缘（最小值）外有5个红色的数据点，右侧也就是上边缘（最大值）外有3个红色的数据点，总计一共有8个异常值。

然后我又通过散点图的方式来可视化异常值的分布情况，如下我做了一些定制，为每个异常值标注数值信息。最终我们可以发现一共有个8个异常值，因为这8个异常值属于比最大值还要大，所以在上限上方，显然这与之前的箱型图出现的异常值数量相吻合。

箱型图 代码实现

代码部分如下 （其中data数据沿用最早的data模拟数据），或许你会问为什么同样的数据得出的异常值个数不同呢，通常来说通过箱型图的以分位数做判断的方式得到的异常值比标准偏差更精确细腻，所以异常值得到的更多。毕竟很少会有大于3倍标准偏差的数据，那么通过3σ准则分析得到的异常值自然就少的多。

### 构建一个箱型图 ###
fig = plt.figure(figsize = (12, 3))  
ax = fig.add_subplot(111)  
# vert = True 时 箱型图为垂直方式
red_circle = dict(markerfacecolor='r', marker='o')  
data.plot.box(vert = False, grid = True, ax = ax,  
              title='箱型图', flierprops=red_circle)



### 构建一个散点图 ###
# 通过describe()函数获得数据的统计量描述信息
data_des = data.describe()  
q1 = data_des['25%']  
q3 = data_des['75%']  
iqr = q3 - q1  
minimum = q1 - 1.5 * iqr  
maximum = q3 + 1.5 * iqr  
outlier = data[(data < minimum) | (data > maximum)]  
correct = data[(data >= minimum) | (data <= maximum)]

print('共有{}个异常值'.format(len(outlier)))

fig = plt.figure(figsize = (12, 4))

plt.scatter(correct.index, correct.values, s = 50,color = 'green', edgecolor = 'black', alpha = 0.5, label='正常值')  
plt.scatter(outlier.index, outlier.values, color = 'red', s = 50, edgecolor = 'black', label='异常值')

# 为异常值附上文字说明
for x, y in zip(outlier.index, outlier.values):  
    plt.text(x = x + 10, y = y - 0.15, s = '{:0.2f}'.format(y), color = 'black', fontsize = 12)

# 添加x横轴标线
plt.axhline(y = maximum, color = 'blue', linestyle = '--', label = '上限')  
plt.axhline(y = minimum, color = 'red', linestyle = '--', label = '下限')

plt.ylim(-5, 5)  
plt.legend(loc = 8, ncol = 4);  

Series是一种类似于一维数组的对象，它由一组数据和一组与之相关的数据标签索引组成。DataFrame是一个表格型的数据结构，它既有行索引也有列索引，它可以被看做有Series组成的字典。

一、将Series转换为DataFrame数据结构

而将Series转换为DataFrame也是会经常遇到的，看下常用的几种方式：

• 直接使用字典方式

  使用字典的方式，创建好对应的列名与字典值就可以了。

import numpy as np  
import pandas as pd

# 将Series转换为DataFrame
data = pd.Series(np.random.randn(10)*500+1000,  
                 index=['A37','A50','R7S','Note5',
                        'G7','R9_Plus','5C','X5_Pro','MX3','M5'])

df = pd.DataFrame({'Product_Name':data.index, 'Price':data.values})  
print(df)

# 输出内容如下：
#   Product_Name        Price
# 0          A37  1000.464248
# 1          A50   657.992057
# 2          R7S   288.075879
# 3        Note5  2094.297636
# 4           G7  1886.582215
# 5      R9_Plus  1843.629256
# 6           5C   931.553668
# 7       X5_Pro   880.684009
# 8          MX3   932.247291
# 9           M5   650.789302

• 使用reset_index()方法转换

import numpy as np  
import pandas as pd

# 将Series转换为DataFrame
data = pd.Series(np.random.randn(10)*500+1000,  
                 index=['A37','A50','R7S','Note5',
                        'G7','R9_Plus','5C','X5_Pro','MX3','M5'])

# 通过reset_index()转换为DataFrame
df = data.reset_index(name='Price')

print(type(df))  
print(df)  

从如下输出可以发现，经过reset_index()方法转换后输出的df已经由一个Series结构类型转换为DataFrame，并且原来的Series标签索引列转换成了DataFrame下名为index的一列。

注：和之前使用字典方式转换不同的是我们还需要后续额外修改列名。

<class 'pandas.core.frame.DataFrame'>

     index        Price
0      A37   629.039268  
1      A50    49.832131  
2      R7S  1693.139415  
3    Note5   445.894181  
4       G7  1133.240339  
5  R9_Plus   886.009704  
6       5C  1115.016608  
7   X5_Pro   894.055677  
8      MX3  1014.297011  
9       M5   520.398571  

二、修改DataFrame的列名

由于之前我们使用reset_index()方法将Series转换为DataFrame后，但是转换后的列名还需要修改，这里可以通过使用DataFrame.rename(mapper=None, index=None, columns=None, axis=None, copy=True, inplace=False, level=None)方法修改。除了修改columns还可以修改index，十分方便。

import numpy as np  
import pandas as pd

# 将Series转换为DataFrame
data = pd.Series(np.random.randn(10)*500+1000,  
                 index=['A37','A50','R7S','Note5',
                        'G7','R9_Plus','5C','X5_Pro','MX3','M5'])

df = data.reset_index(name='Price')  
print(df)

# 输出内容如下：
#          index        Price
# 0          A37  1000.464248
# 1          A50   657.992057
# 2          R7S   288.075879
# 3        Note5  2094.297636
# 4           G7  1886.582215
# 5      R9_Plus  1843.629256
# 6           5C   931.553668
# 7       X5_Pro   880.684009
# 8          MX3   932.247291
# 9           M5   650.789302

df.rename(columns={'index':'Product_Name'}, inplace=True)  
print(df)

# 输出内容如下：
#   Product_Name        Price
# 0          A37  1000.464248
# 1          A50   657.992057
# 2          R7S   288.075879
# 3        Note5  2094.297636
# 4           G7  1886.582215
# 5      R9_Plus  1843.629256
# 6           5C   931.553668
# 7       X5_Pro   880.684009
# 8          MX3   932.247291
# 9           M5   650.789302

先来说说什么是帕累托法则，其原型是19世纪意大利经济学家帕累托所创的库存理论。帕累托运用大量的统计资料分析当时的一些社会现象，概括出一种关键的少数和次要的多数的理论，并根据统计数字画成排列图，后人把它称为帕累托曲线图。简单的说，帕累托法则其实就是我们常说的二八法则，在经济学定律中说的是80%的财富掌握在20%的人手中，而在运营中说的则是80%的贡献度来自于20%的用户。

而基于帕累托法则的帕累托分析法（Pareto Analysis）是制定决策的统计方法，用于从众多任务中选择有限数量的任务以取得显著的整体效果。

下图是我基于Pandas随机模拟出的一组产品数据，并结合Matplotlib来展现这组数据的相关帕累托曲线图。

如图可见，蓝色的条形柱状图代表了每款产品的销售金额，每一个橘色百分点标记代表了从左往右的销售累计占比。其中红色的垂直虚线则是一条数据占>=80%的分界线，也就代表了这条分界线之前的7款产品（即，'G7', 'MX3', 'R9Plus' , 'A37', 'Note5', 'X5Pro', 'A50'）占到了总销售额中83.91%的销售份额，之后其余的3款产品占比只有16.09%。假如放到现实场景中，就可以根据市场情况制定某些决策，如是否要撤销某些产品的研发。

接下来为了更好的理解代码是如何实现的，我将分解代码成若干段来说说明。首先我用Pandas模拟10组虚拟产品数据。

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

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

# 创建模拟数据，10个品类产品的销售额     
data = pd.Series(np.random.randn(10)*500+1000,  
                 index=['A37','A50','R7S','Note5',
                        'G7','R9_Plus','5C','X5_Pro',
                        'MX3','M5'])
print(data)  

通过输出data变量值查看如下，左侧的index是我们虚拟的10组产品名称，而右侧对应的数值则是通过random.randn随机产生的10组基于标准正态分布的伪数据。

A37        1784.729567  
A50         198.971079  
R7S         985.258417  
Note5      1348.619119  
G7         1330.237773  
R9_Plus     188.307306  
5C         2105.157302  
X5_Pro     1225.447782  
MX3         957.038300  
M5         1346.797999  
dtype: float64  

接着，我们将所有的数据从大到小降序排列并以此构建柱状图。

其中要注意的一点是，为什么需要降序排列呢？因为考虑到我们要制作的是累计占比，同样是计算累计相加值满足>=80%，从小到大升序排列 与 从大到小降序排列 相比 降序排列能以更少的数据个数预先满足累计相加值>=80%，也就能体现帕累托法则的二八法则，少数掌握多数的原理。

# 将数值按照从大到小排序
data.sort_values(ascending=False, inplace=True)  
# 构建画布
plt.figure(figsize=(10,6))

# 先构建一个从大到小排列的，标准柱状图。
data.plot(kind='bar', colormap='GnBu_r', width=0.6, edgecolor='black', rot=0)

plt.xlim(-1,11)  
plt.ylim(0,2000)

plt.xlabel('产品名')  
plt.ylabel('销售金额')  
plt.title('销售情况')  
plt.legend(['销售金额'], loc='upper right')  

输出图如下：

之后，需要绘制另一条累计占比的百分比图。首先需要通过cumsum函数得到每个产品时的当前累计值，并用这个累计值去除以总销售金额得到的就是该产品当前的累计占比的百分比数值。得到每个产品当前的累计占比的百分比数值后就可以通过布尔索引筛选出>=80%的产品名称。有了这个产品名称就可以获得该产品的索引位置编号，而这个索引位置十分重要，因为这有助于我们之后用来标注那条红色分界线的具体位置。

# 累加值 除以 总数值 得到每个阶段累加值的占比值。
proportion = data.cumsum() / data.sum()

# 用布尔索引列出大于等于0.8的数值列表，而列表里第一行就是大于80%的关键点，我们取得它的索引名称。
key = proportion[proportion >= 0.8].index[0]

# 根据之前得到的索引名称，就可以从data数据中得到该名称的索引位置。
key_number = data.index.tolist().index(key)

proportion.plot(style='--ko', secondary_y=True, color='#FFA500')

plt.axvline(key_number, linestyle='--', color='red')

for x, y, s in zip(range(len(data)), proportion.values, proportion):  
    if x == key_number:
        plt.text(key_num+0.2, y-0.035, '累计占比为 {:.2%}'.format(s), color='red', fontsize=12)
    else:
        plt.text(x+0.3, y-0.035, "{:.2%}".format(s), fontsize=12, horizontalalignment='center', color='black')

plt.ylabel('销售（比例）')  
plt.legend(['累积占比'], loc='upper center')  

之后将两段代码合并，也就是两张图合并为一张图就形成了我们一开始提到的帕累托曲线图。

最后放上完整代码如下：

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

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

# 创建模拟数据，10个品类产品的销售额     
data = pd.Series(np.random.randn(10)*500+1000,  
                 index=['A37','A50','R7S','Note5',
                        'G7','R9_Plus','5C','X5_Pro',
                        'MX3','M5'])


# 将数值按照从大到小排序
data.sort_values(ascending=False, inplace=True)  
# 构建画布
plt.figure(figsize=(10,6))

# 先构建一个从大到小排列的，标准柱状图。
data.plot(kind='bar', colormap='GnBu_r', width=0.6, edgecolor='black', rot=0)  
plt.xlim(-1,11)  
plt.ylim(0,2000)  
plt.xlabel('产品名')  
plt.ylabel('销售金额')  
plt.title('销售情况')  
plt.legend(['销售金额'], loc='upper right')

# 累加值 除以 总数值 得到每个阶段累加值的占比值。
proportion = data.cumsum() / data.sum()

# 用布尔索引列出大于等于0.8的数值列表，而列表里第一行就是大于80%的关键点，我们取得它的索引名称。
key = proportion[proportion >= 0.8].index[0]

# 根据之前得到的索引名称，就可以从data数据中得到该名称的索引位置。
key_number = data.index.tolist().index(key)

proportion.plot(style='--ko', secondary_y=True, color='#FFA500')

plt.axvline(key_number, linestyle='--', color='red')

for x, y, s in zip(range(len(data)), proportion.values, proportion):  
    if x == key_number:
        plt.text(key_num+0.2, y-0.035, '累计占比为 {:.2%}'.format(s), color='red', fontsize=12)
    else:
        plt.text(x+0.3, y-0.035, "{:.2%}".format(s), fontsize=12, horizontalalignment='center', color='black')

plt.ylabel('销售（比例）')  
plt.legend(['累积占比'], loc='upper center')  

总结

画帕累托图最重要的几点是要将数据从大到小排序，这主要是为了能体现帕累托法则的二八法则，少数掌握多数的原理。当数据以降序的方式排列后，计算每个列数据累计占比百分比数值。有了百分比数值，我们就可以将数据可视化，可视化后可以很清楚的从图中来区分哪些对应的产品满足了>=80%的贡献。

• 什么是生成器？

生成器其实是一种特殊的迭代器(iterator)，但是不需要像迭代器一样自己去实现__iter__()和__next__()方法，简单的说生成器是通过一个或多个yield表达式构成的函数，生成器是为迭代器产生数据的。如果一个函数包含yield关键字，这个函数就会变为一个生成器。生成器并不会一次返回所有结果，而是每次遇到yield关键字后返回相应结果，并保留函数当前的运行状态，等待下一次的调用。

由于 生成器(generator) 自动实现了迭代器协议，而迭代器协议对很多人来说，也是一个较为抽象的概念。所以，为了更好的理解生成器，我们需要简单的梳理一下迭代器协议的概念。

迭代器协议是指：对象需要提供next方法，它要么返回迭代中的下一项，要么就引起一个StopIteration异常，以终止迭代。

可迭代对象就是：实现了迭代器协议的对象。

协议是一种约定，可迭代对象实现迭代器协议，Python的内置工具 (如for循环，sum，min，max函数等) 使用迭代器协议访问对象。

• 两种不同的方式提供创建生成器
  • 生成器函数：使用yield关键字语句定义的常规函数就被认为是一个生成器，但而不是return语句返回结果。因为其中的区别是yield语句一次返回一个结果，在每个结果中间，挂起函数的状态，以便下次重它离开的地方继续执行。而return语句一旦退出就真的退出函数体了。
  • 生成器表达式：类似于列表推导，但是生成器返回按需产生结果的一个对象，而不是一次构建一个结果列表。

我们主要说的就是如何使用yeild关键字创建的生成器。如下例子是一个关于计算斐波那契数列的生成器。其中 fibonacci 函数中我们没有用 return 关键字。你也可以看到当运行 fib = fibonacci(20) 的时候，我们打印出它的类型信息显示它返回的是一个生成器对象。如果你直接调用这个（生成器的）实例对象fib，并不会运行 fibonacci 函数中的代码。前面我们提到过，因为生成器其实是一种特殊的迭代器(iterator)，所以我们可以看到在代码最后部分只有当循环调用 next() 的时候才会真正运行其中的代码。而且用这种方式，我们可以不用担心它会使用大量的内存资源。

def fibonacci(n):  
    x, y = 0, 1
    for _ in range(n):
        yield x
        x, y = y, x + y

fib = fibonacci(20)  
print(type(fib))  
# 输出如下：
# <class 'generator'>

for _ in fib:  
    print(next(fib))
# 输出如下：
# 1
# 2
# 5
# 13
# 34
# 89
# 233
# 610
# 1597
# 4181

相反，如果不用生成器，我们用常规函数定义的话就会像下面这样写。如果我们传入的参数n值增大，返回的列表占用的空间将会显著提升，这显然是我们不希望看到的。

def fibonacci(n):  
   i, a, b = 1, 0, 1 
   L = [] 
   while i < n: 
       L.append(b) 
       a, b = b, a + b 
       i = i + 1 
   return L

fib = fibonacci(20)  
print(type(fib))

for f in fib:  
    print(f)

• 判断函数是否是生成器

我们可以用inspect类里的isgeneratorfunction类方法判断是否是一个生成器函数，以及使用 isgenerator类方法判断是否是一个生成器。

from inspect import isgeneratorfunction, isgenerator

print(f'fibonacci is a generator function: {isgeneratorfunction(fibonacci)}')  
print(f'fib is a generator: {isgenerator(fib)}')

# 输出如下：
# fibonacci is a generator function: True
# fib is a generator: True

• 应用生成器的场景与好处
  • 生成器可用于产生数据流，而且并不立刻产生返回值，而是等到被需要的时候才会产生返回值，相当于一个主动拉取的过程(pull)，比如现在有一个日志文件，每行产生一条记录，对于每一条记录，不同部门的人可能处理方式不同，但是我们可以提供一个公用的、按需生成的数据流。
  • 还有做爬虫的时候，爬取大量数据的时候如果使用生成器每次需要的时候执行输出也可以大大降低资源的消耗。

使用生成器的好处当然不仅限于此，让我们来看一下下面的例子，我们打算读取小说《三国演义》的所有文字内容，如果直接对文件对象调用 read() 方法，会导致不可预测的内存占用。好的方法是利用固定长度的缓冲区来不断读取文件内容。而且同时通过 yield来执行每次输出，就可以轻松实现文件读取。

from pathlib import Path

file = Path('三国演义.txt')

def read_file(fpath):  
   BLOCK_SIZE = 1024 
   with file.open(encoding='GB18030') as f: 
       while True: 
           block_content = f.read(BLOCK_SIZE) 
           if block_content: 
               yield block_content 
           else: 
               return

for c in read_file(file):  
    print(c)

我们先来直观的区分这两者有什么不同。

可迭代 (iterable)：如果一个对象具备有__iter__() 或者 __getitem__()其中任何一个魔术方法的话，这个对象就可以称为是可迭代的。其中，__iter__()的作用是可以让for循环遍历，而__getitem__()方法可以让实例对象通过[index]索引的方式去访问实例中的元素。所以，列表List、元组Tuple、字典Dictionary、字符串String等数据类型都是可迭代的。

迭代器 (iterator): 如果一个对象同时有__iter__()和__next__()魔术方法的话，这个对象就可以称为是迭代器。__iter__()的作用前面我们也提到过，是可以让for循环遍历。而__next__()方法是让对象可以通过 next(实例对象) 的方式访问下一个元素。列表List、元组Tuple、字典Dictionary、字符串String等数据类型虽然是可迭代的，但都不是迭代器，因为他们都没有next( )方法。

• 如何判断可迭代(iterable) & 迭代器(iterator)

我们可以借助Python中的isinstance(object, classinfo)函数来判断一个对象是否是一个已知类型。如下例子中，通过isinstance( )函数分别判断列表、元组、字典、字符串是不是可迭代或迭代器。

from collections import Iterable  
from collections import Iterator

print(f"List is 'Iterable': {isinstance([], Iterable)}")  
print(f"Tuple is 'Iterable': {isinstance((), Iterable)}")  
print(f"Dict is 'Iterable': {isinstance({}, Iterable)}")  
print(f"String is 'Iterable': {isinstance('', Iterable)}")

print("="*25)

print(f"List is 'Iterator': {isinstance([], Iterator)}")  
print(f"Tuple is 'Iterator': {isinstance((), Iterator)}")  
print(f"Dict is 'Iterator': {isinstance({}, Iterator)}")  
print(f"String is 'Iterator': {isinstance('', Iterator)}")

# 输出如下：
# List is 'Iterable': True
# Tuple is 'Iterable': True
# Dict is 'Iterable': True
# String is 'Iterable': True
# =========================
# List is 'Iterator': False
# Tuple is 'Iterator': False
# Dict is 'Iterator': False
# String is 'Iterator': False

通过对定义的分析和比较我们得知：迭代器都是可迭代的（因为迭代器都包含__iter__()函数），但可迭代的不一定是迭代器（因为未必每个可迭代就包含__next__()方法）。

• 创建一个迭代器

得益于Python的鸭子类型特性，只要我们实现类似具备某一特征的方法，就可以认为它就是什么。所以我们定义了一个类并在类中实现__iter__()和__next__()方法，那么这个类就可以当做是一个迭代器了。

from collections import Iterator

class Data:  
    def __init__(self, x):
        self.x = x

    def __iter__(self):
        return self

    def __next__(self):
        if self.x >= 10:
            raise StopIteration
        else:
            self.x += 1
            return self.x

data = Data(0)

print(f"data is 'Iterator': {isinstance(data, Iterator)}")

# 输出如下：
# data is 'Iterator': True

如上例子中我们可以看到，最后我们用isinstance()函数判断得到结果为True，证明我们定义的实例对象是一个真正的迭代器了。因为是迭代器，我们就可以用for循环来验证试试。

from collections import Iterator

class Data:  
    def __init__(self, x):
        self.x = x

    def __iter__(self):
        return self

    def __next__(self):
        if self.x >= 10:
            raise StopIteration
        else:
            self.x += 1
            return self.x

data = Data(0)

for d in data:  
    print(d)

# 输出如下：
# 1
# 2
# 3
# 4
# 5
# 6
# 7
# 8
# 9
# 10

上述例子中，我们定义的类对象内部，把x的值显示在大于等于10以内，否则就会抛出StopIteration异常错误（当然实际使用时并不会出错，只是中断继续执行。）我们先创建了一个初始值为0的实例对象，最后顺利的用for循环遍历了从1到10的数字，因为内部对大于等于10的限制，所以输出到10的时候就停止了。特别要注意的是，如果你再次单独去执行for循环的话不会有任何输出，因为迭代器默认只运行一次。

除了自己定义__iter__()和__next__()魔术方法的外，我们还可以使用Python内置的iter()函数来返回一个迭代器，像下面这样。

list_a = [1,2,3,4,5,6]  
my_iterator = iter(list_a)

print(f"my_iterator is 'Iterator': {isinstance(my_iterator, Iterator)}")

# 输出如下：
# my_iterator is 'Iterator': True

我们知道，迭代器必须具备两个基本方法__iter__()和__next__()，而__next__()方法是让对象可以通过 *next(实例对象) *的方式访问下一个元素。所以让我们验证下用next()的方式去访问这个我们转换过的迭代器是否能正常运行。

next(my_iterator)

# 输出如下：
# 1
# 2
# 3
# 4
# 5
# 6

通过 next(实例对象)的方式可以访问出每个元素。但是这里要特别说明的一点是，next()函数只能每次执行一次才输出一次结果，如上从1输出到6，是我们手动执行了6次分别执行出来的，如果在输出6后再次执行的话，它就会报一个StopIteration异常错误。

最后，我们还可以使用Python内置的dir()函数来看看传入参数的属性，方法等信息，比如我们用它来看看之前从list转换成的my_iterator迭代器。

dir(my_iterator)

# 输出如下：
# ['__class__',
#  '__delattr__',
#  '__dir__',
#  '__doc__',
#  '__eq__',
#  '__format__',
#  '__ge__',
#  '__getattribute__',
#  '__gt__',
#  '__hash__',
#  '__init__',
#  '__init_subclass__',
#  '__iter__',
#  '__le__',
#  '__length_hint__',
#  '__lt__',
#  '__ne__',
#  '__new__',
#  '__next__',
#  '__reduce__',
#  '__reduce_ex__',
#  '__repr__',
#  '__setattr__',
#  '__setstate__',
#  '__sizeof__',
#  '__str__',
#  '__subclasshook__']

可以发现，也是意料之中的，my_iterator迭代器包含了两个基本方法__iter__()和__next__()方法。