学习一下神经网络。其在当下之热，应用之广就不必赘述。

就群体和进化遗传学来说，通过神经网络识别自然选择信号也是这两年的一个研究方向。通过计算机模拟，计算很多传统识别自然选择信号的统计量，比如Tajima's D 之类的，然后构建一个二维向量矩阵，然后使用卷积神经网络（CNN）来训练，识别，判断一个基因组区域是否受到自然选择。

本文初试一下神经网络，相当于一个学习笔记。

导入iris数据集

我们从sklearn中导入iris数据集，然后我们简单查看以下该数据集的简单情况

from sklearn import datasets

iris = datasets.load_iris()
iris.keys()

dict_keys(['data', 'target', 'frame', 'target_names', 'DESCR', 'feature_names', 'filename', 'data_module'])

iris.data.shape

(150, 4)

iris.data[:5,]

array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2]])

该数据集是一个150x4的格式，即有150行观测，有4列变量，这4列变量分别描述花萼长度，宽度，花瓣长度，宽度，如下：

iris.feature_names

['sepal length (cm)',
 'sepal width (cm)',
 'petal length (cm)',
 'petal width (cm)']

我们在看以下这150行观测数据对应的因变量（Y）花的种类，如下：

iris.target

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

其中0 1 2 对应的名字如下：

iris.target_names

array(['setosa', 'versicolor', 'virginica'], dtype='<U10')

我们的目的就是通过花萼和花瓣的四个变量也预测是哪一种花。但是四个变量的尺度不同，这种情况下，我们一般需要对数据进行标准化，各个变量的分布均在0-1之间。我们可以通过下面的标准化公式，对每一个变量进行标准化。

当然，也可以使用sklearn里面的功能对整个数据集进行标准化，如下：

from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler()
iris.data = scaler.fit_transform(iris.data)

对于因变量（Y）我们同样需要进行处理，此处变量类型为无需（三）多分类，我们可以将其转换为哑变量。keras.utils中的功能可以很方便地将多分类变量转很为哑变量的形式。

from keras.utils import to_categorical
y_dummy = to_categorical(iris.target)
y_dummy[(1,90,140),]

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]], dtype=float32)

在建模之前，将数据集分为训练集train和测试集test。可以通过随机抽样，把整个iris数据集分开。不过，sklearn已经为我们准备好了更方便的划分函数train_test_split。此处，我们用80% （150 x 0.8 = 120）数据训练,剩下的20% （30）用来测试。

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(iris.data, y_dummy, test_size=0.2)

print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)

(120, 4) (30, 4) (120, 3) (30, 3)

建立神经网络

因为我们有4个变量，所以输入数据维度为4，然后构建256个神经元组成的一个层，当然，也可以是其他数量的神经元，看自己爱好。activation在中间隐藏层一般用relu函数。随后的输出层有3个神经元，对应我们的预测变量。我们的预测变量为二分类变量，这儿activation使用sigmoid函数。

optimizer有很多种，可以自行尝试其他类型的优化。因为我们的预测变量为二分类，所以损失函数使用分类交叉熵函数。

callbacks是为了保留整个训练优化的过程，如果你对训练过程不感兴趣，可以省略掉。

import tensorflow as tf
from tensorflow import keras

model = keras.Sequential([
    keras.layers.Dense(256, input_shape=(4,), activation='relu'),
    keras.layers.Dense(3, activation='sigmoid')
])

model.compile(optimizer='adam',
             loss='categorical_crossentropy',
             metrics=['accuracy'])

tb_callback = tf.keras.callbacks.TensorBoard(log_dir='tf_logs/', histogram_freq=1)

with tf.device('GPU:0'):
    model.fit(X_train, y_train, epochs=50, callbacks=[tb_callback])

Epoch 1/50
4/4 [==============================] - 0s 8ms/step - loss: 1.1937 - accuracy: 0.5583
Epoch 2/50
4/4 [==============================] - 0s 9ms/step - loss: 1.0127 - accuracy: 0.5083
Epoch 3/50
4/4 [==============================] - 0s 8ms/step - loss: 0.9612 - accuracy: 0.3417

....

Epoch 48/50
4/4 [==============================] - 0s 7ms/step - loss: 0.2488 - accuracy: 0.9667
Epoch 49/50
4/4 [==============================] - 0s 8ms/step - loss: 0.2439 - accuracy: 0.9583
Epoch 50/50
4/4 [==============================] - 0s 7ms/step - loss: 0.2384 - accuracy: 0.9500

另外，此处还用到了GPU加速（本例子数据量小，没有必要用）。可以通过下面的例子看看是否可以使用GPU。如果此处没有显示出‘GPU’，但是你自己的电脑有GPU硬件，那么你需要安装cuda库和cuDNN库，同时TensorFlow必须是支持GPU版本的。

一般而言，GPU能够大大加快模型训练速度，可能有十几上百倍，不过本文测试只有提升了数倍的速度。可能是本例子训练量小？或者设备GPU比较弱？

此外，callbacks的结果我们保存在了tf_logs的文件夹里，在终端输入tensorboard.exe --logdir .\tf_logs (Windows)，你会得到一个类似 http://localhost:6006/的地址，将该地址粘贴到浏览器中，即可以查看训练的详细过程。此处不再展示。

tf.config.experimental.list_physical_devices()

[PhysicalDevice(name='/physical_device:CPU:0', device_type='CPU'),
 PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]

上面的神经网络模型进行了50次训练，最终预测的准确度达到了95%。

模型的评价和预测

下面在我们的测试集上验证一下我们模型的预测能力。

model.evaluate(X_test, y_test)

1/1 [==============================] - 0s 87ms/step - loss: 0.2453 - accuracy: 1.0000


[0.24527209997177124, 1.0]

此处，模型的预测准确度为100%，这是一个非常理想的结果，可能测试数据量比较小。这也和本身iris数据集有关，如果你做主成分分析，会看到主成分能够将三类花很好的区分开来。下面是我们具体的预测结果

import numpy as np
y_pred = np.argmax(model.predict(X_test),axis=1)
iris.target_names[y_pred]

1/1 [==============================] - 0s 30ms/step


array(['setosa', 'setosa', 'virginica', 'versicolor', 'virginica',
       'virginica', 'versicolor', 'virginica', 'setosa', 'setosa',
       'versicolor', 'versicolor', 'versicolor', 'versicolor',
       'versicolor', 'virginica', 'setosa', 'setosa', 'versicolor',
       'versicolor', 'virginica', 'virginica', 'setosa', 'virginica',
       'versicolor', 'versicolor', 'setosa', 'setosa', 'versicolor',
       'virginica'], dtype='<U10')

对比上面的模型预测结果，我们看一下实际的花的种类（如下），可以看到二者完全一致。

y_test_raw = np.argmax(y_test, axis=1)
iris.target_names[y_test_raw]

array(['setosa', 'setosa', 'virginica', 'versicolor', 'virginica', 'virginica', 'versicolor', 'virginica', 'setosa', 'setosa', 'versicolor', 'versicolor', 'versicolor', 'versicolor', 'versicolor', 'virginica', 'setosa', 'setosa', 'versicolor', 'versicolor', 'virginica', 'virginica', 'setosa', 'virginica', 'versicolor', 'versicolor', 'setosa', 'setosa', 'versicolor', 'virginica'], dtype='<U10')

下面是对模型预测的总结：

from sklearn.metrics import confusion_matrix, classification_report
print(classification_report(y_test, y_pred))

              precision    recall  f1-score   support

           0       1.00      1.00      1.00         8
           1       1.00      1.00      1.00        10
           2       1.00      1.00      1.00        12

    accuracy                           1.00        30
   macro avg       1.00      1.00      1.00        30
weighted avg       1.00      1.00      1.00        30

在深度学习中，混淆矩阵confusion matrix是我们常用来查看预测准确性的图。下图即是我们模型预测的混淆矩阵。

import seaborn  as sn
import matplotlib.pyplot  as plt
cm = tf.math.confusion_matrix(labels=y_test_raw, predictions=y_pred)
sn.heatmap(cm, annot=True, fmt='d')
plt.xlabel('Prediction')
plt.ylabel('Truth')
plt.xticks(np.arange(3)+0.5, iris.target_names)
plt.yticks(np.arange(3)+0.5, iris.target_names)

([<matplotlib.axis.YTick at 0x2a5a7bdbb50>,
  <matplotlib.axis.YTick at 0x2a5a7bdb3d0>,
  <matplotlib.axis.YTick at 0x2a5a8918340>],
 [Text(0, 0.5, 'setosa'),
  Text(0, 1.5, 'versicolor'),
  Text(0, 2.5, 'virginica')])

从上图可以看到，9个setosa, 12个versicolor和9个virginica全部准确预测。